MATLAB音频信号的分析与合成程序.ppt

1,第二讲 MATLAB的,21 脚本文件和函数文件 22 函数调用和参数传递 23 MATLAB的程序结构和控制流 24 M文件的调试 25 程序设计实例 251 音频信号的分析与合成 252 音频频谱分析仪 253 幅度调制,程序设计,2,21 脚本文件和函数文件,211 M脚本文件,对于一些比较简单的问题 ，在指令窗中直接输入指令计算 。,对于复杂计算，采用脚本文件（Script file）最为合适 。,MATLAB只是按文件所写的指令执行 。,M脚本文件的特点是：,脚本文件的构成比较简单，只是一串按用户意图排列而成的（包括控制流向指令在内的）MATLAB指令集合。,脚本文件运行后 ，所产生的所有变量都驻留在 MATLAB基本工作空间（Base workspace）中。只要用户不使用清除指令（clear）， MATLAB指令窗不关闭，这些变量将一直保存在基本工作空间中。,M文件有两种形式 ：脚本文件（Script File）和函数文件（Function File ）。这两种文件的扩展名，均为“ . m” 。,3,21 脚本文件和函数文件（续1）,212 M函数文件,与脚本文件不同 ，函数文件犹如一个“黑箱”，把一些数据送进并经加工处理，再把结果送出来。,MATLAB提供的函数指令大部分都是由函数文件定义的。,M函数文件的特点是：,从形式上看 ，与脚本文件不同 ，函数文件的笫一行总是以 “function”引导的“函数申明行”。,从运行上看 ，与脚本文件运行不同 ，每当函数文件运行， MATLAB就会专门为它开辟一个临时工作空间，称为函数工作空间（ Function workspace） 。当执行文件最后一条指令时 ，就结束该函数文件的运行，同时该临时函数空间及其所有的中间变量就立即被清除。,MATLAB允许使用比 “标称数目 ”较少的输入输出宗量，实现对函数的调用 。,4,21 脚本文件和函数文件（续2）,213 M文件的一般结构,由于从结构上看 ，脚本文件只是比函数文件少一个“函数申明行”，所以只须描述清楚函数文件的结构 。,典型 M函数文件的结构如下 ：,函数申明行：位于函数文件的首行，以关键字 functio 开头，函数名以及函数的输入输出宗量都在这一行被定义。,笫一注释行：紧随函数申明行之后以%开头笫一注释行。该行供lookfor关键词查询和 help在线帮助使用 。,在线帮助文本区 ：笫一注释行及其之后的连续以%开头的所有注释行构成整个在线帮助文本。,编写和修改记录：与在线帮助文本区相隔一个“空”行，也以%开头，标志编写及修改该M文件的作者和日期等 。,函数体：为清晰起见，它与前面的注释以“空”行相隔。,例 2_2_1：M函数文件示例。,5,22 函数调用和参数传递,221 局部变量和全局变量,局部（Local）变量：它存在于函数空间内部的中间变量，产生于该函数的运行过程中，其影响范围也仅限于该函数本身 。,全局（Global）变量：通过 global 指令，MATLAB也允许几个不同的函数空间以及基本工作空间共享同一个变量，这种被共享的变量称为全局变量。,222 函数调用,在MATLAB中，调用函数的常用形式是：,输出参数1,输出参数2, = 函数名(输入参数1,输入参数2, ),函数调用可以嵌套，一个函数可以调用别的函数，甚至调用它自己 （递归调用）。,6,22 函数调用和参数传递（续）,223 参数传递,MATLAB在函数调用上有一个与众不同之处 ：函数所传递的参数具有可调性 。,传递参数数目的可调性来源于如下两个MATLAB永久变量：,函数体内的 nargin 给出调用该函数时的输入参数数目。,函数体内的 nargout 给出调用该函数时的输出参数数目。,只要在函数文件中包括这两个变量，就可以知道该函数文件调用时的输入参数和输出参数数目。,值得注意：nargin、 nargout 本身都是函数，不是变量，所以用户不能赋值，也不能显示。,“变长度”输入输出宗量：varargin 、 varrgout。具有接受 “任意多输入” 、返回“任意多输出”的能力 。,跨空间变量传递：evalin。,（参考 ： circle .m，am1.m）,7,23 MATLAB的程序结构和控制流,231 程序结构,循环结构： MATLAB提供两种循环方式。,顺序结构,分支结构： ifelseend 。,forend 循环和while-end循环。,232 程序流控制,常用指令：return，echo，input，pause，keyboard，break。,switch-case 结构。,try-catch 结构。,警示指令：error， warning 。,8,23 MATLAB的程序结构和控制流（续）,233 图形用户界面（GUI）编程,现代的主流应用程序已经从命令行的交互方式转变为以图形界面为主的交互方式，这主要是由于它给用户带来了操作和控制的方便与灵活性。（面向对象编程）,MATLAB能够以比较简单的方式实现一系列的图形界面功能。通过对控件、菜单属性的设置和 Callback 的编写，就能够满足大多数用户的需求。,控件的 Callback 属性： Callback 属性的取值是字符串，可以是某个M文件名或一小段MATLAB语句。当用户激活控件对象（例如 ：在控件对象图标上单击鼠标左键 ）时，应用程序就运行 Callback 属性定义的子程序。,菜单的 Callback 属性： Callback 属性的取值是字符串，可以是某个M文件名或一小段MATLAB语句。当用户激活菜单对象时，若没有子菜单就运行 Callback 属性定义的子程序。若有，先运行 Callback 属性定义的子程序，再显示子菜单。,9,24 M文件的调试,编写 M文件时，错误（Bug）在所难免。错误有两种：语法（Syntax）错误和运行（Run-time）错误。,语法错误是指变量名、函数名的误写，标点符号的缺、漏等。对于这类错误，通常能在运行时发现，终止执行，并给出相应的错误原因以及所在行号。,运行错误是算法本身引起的，发生在运行过程中。相对语法错误而言，运行错误较难处理 。尤其是M函数文件，它一旦运行停止，其中间变量被删除一空，错误很难查找。,有两种调试方法：直接调试法和工具调试法。,10,24 M文件的调试（续1）,直接调试法：可以用下面方法发现某些运行错误。,在M文件中，将某些语句后面的分号去掉， 迫使M文件输出一些中间计算结果，以便发现可能的错误。,在适当的位置，添加显示某些关键变量值的语句（包括使用 disp 在内）。,利用 echo 指令，使运行时在屏幕上逐行显示文件内容。echo on 能显示M脚本文件；echo FunNsme on 能显示名为FunNsme 的M函数文件。,在原M脚本或函数文件的适当位置，增添指令 keyboard 。 keyboard 语句可以设置程序的断点 。,通过将原M函数文件的函数申明行注释掉，可使一个中间变量难于观察的M函数文件变为一个所有变量都保留在基本工作空间中的M脚本文件。,11,24 M文件的调试（续2）,GUI 界面调试法：,MATLAB 5.x 版提供了一个基于GUI界面的调试。使用它，可以对函数进行调试。,Debug菜单的使用：,Continue：恢复程序运行至结束或另一个断点 。,Single Step：单步执行函数。,Step In：深入下层局部工作区 。,Quit Debugging：退出调试状态。,Set/Clear Breakpoint：设置/清除光标处的断点 。,Clear All Breakpoints：清除程序中的所有断点 。,Stop if Error：运行至出错或结束。,Stop if Warning：运行至警告消息或结束。,Stop if NaN of Inf：运行至运算结果出现 NaN 或 Inf。,12,25 程序设计实例,251 音频信号的分析与合成,（韩利竹P262-sigspec.m，stereospec.m）,采用MATLAB分析WAV文件。,步骤一：选择一个WAV文件作为分析的对象 。（ding.wav）,步骤二：读WAV文件数据并画时域图形 。,步骤三：进行 FFT变换并画频域图形 。,步骤四：进行该声波主要频谱的分析。,步骤五：根据该声音的频谱，反演时域图形 。（有失真）,常用函数：,x，fs，bits=waveread(filename),d= FFT (w，l),sound(w，fs，bits),步骤六：进行付立叶逆变换IFFT并画频域图形 。,13,252 音频频谱分析仪,（车晴P199-audspec.m）,M函数文件 audspec.m 用于对音频频谱分析仪进行仿真，它通过鼠标来选择一个WAV文件，它可以播放该声音文件、显示音频输入信号的波形和频谱。,WAV文件的输入方式采用 uigetfile 函数，该函数打开一个标准的对话框，然后可以选择声音文件的路径和文件的名称。,信号波形生成使用绘图函数就可以完成，而信号的频谱是通过快速付立叶变换而得到的。,常用函数：,strcmp(action,initialized);,strcmp(action,Play);,strcmp(action,Zoom);,strcmp(action,Spectrum);,14,253 幅度调制,（车晴P203- am.m）,幅度调制，即载波幅度随调制信号变化的调制（AM）。幅度调制有多种实现方式，包括：标准幅度调制（SAM），双边带幅度调制（DSB-AM），单边带幅度调制（SSB），残留边带调制（VSB）和平衡正交幅度调制（QAM）等。,假定调制信号为m(t)，其幅度为1v ，载波信号为 ，各种调幅方式的已调波的表达式如下 :,SAM:,是调制度。,DSB-AM：,SSB：,是调制信号 的希尔伯特变换。,QAM：,15,M脚本文件入门,返回,通过M脚本文件，画出下列分段函数所表示的曲面。,（s_file.m）,16,M脚本文件入门（续）,返回,%s_file.m This is my first example. a=2;b=2; clf; x=-a:0.2:a;y=-b:0.2:b; for i=1:length(y) for j=1:length(x) if x(j)+y(i)1 z(i,j)=0.5457*exp(-0.75*y(i)2-3.75*x(j)2-1.5*x(j); elseif x(j)+y(i)=-1 z(i,j)=0.5457*exp(-0.75*y(i)2-3.75*x(j)2+1.5*x(j); else z(i,j)=0.7575*exp(-y(i)2-6.*x(j)2); end end end axis(-a,a,-b,b,min(min(z),max(max(z); colormap(flipud(winter);surf(x,y,z);,17,例 2_2_1：M函数文件示例。,返回,sa = circle(r,s) % CIRCLE plot a circle of radii r in the line specified by s. % r 指定半径的数值 % s 指定线色的字符串 % sa 圆面积 % circle(r) 利用蓝实线画半径为 r 的圆周线. % circle(r,s) 利用串 s 指定的线色画半径为 r 的圆周线. % sa=circle(r) 计算圆面积，并画半径为 r 的蓝色圆面. % sa=circle(r,s) 计算圆面积，并画半径为 r 的 s 色圆面. % 编写于1999年4月7日，修改于1999年8月27日。张志涌EX8.4.4.-1（P299） 程序文本。,18,控件的 Callback 属性,返回,M2_ex2_3_1.m（张志涌P426）,对于传递函数为 的归一化二阶系统，制作一 个能绘制该系统单位阶跃响应的图形用户界面。,在编辑框输入 0.3，可以看到 一根曲线 。 在编辑框输入 0.1:0.1:1，可以看到十根曲线 。,本例设置断点演示：（A）图形界面的大致生成过程；,（B）静态文本和编辑框的生成；,（C）坐标方格控制键的形成；,（D）如何使用该界面。,运行程序 sub_of_M2_ex2_3_1.m， 验证 Callback 属性定义的子程序。,19,菜单的 Callback 属性,返回,本例的目标是：在图形窗上自制一个名为【Test】的“顶层菜单项”；当用鼠标点动该菜单项时，将产生一个带分格的封闭坐标轴。通过本例说明：,（A）回调属性的运作机理；,（B）用户顶层菜单项的制作；,（C）uimenu属性的设置方法；,（D）复杂字符串的构成方法和注意事项。,M2_ex2_3_2.m （张志涌P432）,20,音频信号的分析与合成程序,（韩利竹P262-sigspec.m，stereospec.m）,返回,y=sigspec(action); % SIGNAL SPECTRUM ANALYSER THEE 2001/10 % y=sigspec(action); % Usage:y=sigspec; %= clear; if nargin1;action=initialized;end; fname,pname=uigetfile(*.wav,Open Wave File); file=pname,fname; x,fs,bits=wavread(file); % 读入声音文件（*.wav） sound(x,fs,bits); % 数据通过声卡转换为声音 %=选择声道= %y=x(:); % 双声道信号波形数据 y=x(:,1); % 左声道信号波形数据 %y=x(:,2); % 右声道信号波形数据 %=,21,音频信号的分析与合成程序（续1）,%= disp(按任意键显示左声道声音波形！); disp(-); pause plot(y); set(gca,FontName,Arial,FontSize,10); title(左声道声音波形); xlabel(采样点数); disp(按任意键听左声道信号三次！); disp(-); pause sound(y,fs,bits); % 数据通过声卡转换为声音 pause sound(y,fs,bits); pause sound(y,fs,bits); k=fft(y,32768); disp(按任意键显示左声道声音频域的幅值！); disp(-);,22,音频信号的分析与合成程序（续2）,disp(-); pause plot(abs(k); set(gca,FontName,Arial,FontSize,10); title(左声道声音频域的幅值); xlabel(采样点数); m1,i1=max(abs(k); % 找出频域最大值 F1=i1/32768*fs; F1t=num2str(F1) Hz; text(1500,450,F1t); m2,i2=max(abs(k(4000:5000); % 找出频域次大值 F2=(4000+i2)/32768*fs; F2t=num2str(F2) Hz; text(3000,40,F2t); m3,i3=max(abs(k(10000:12000); % 找出频域次次大值 F3=(10000+i3)/32768*fs; F3t=num2str(F3) Hz; text(8000,20,F3t); %=,23,音频信号的分析与合成程序（续3）,%= t=(0:0.0001:1); y1=(sin(2*pi*788*t)+sin(2*pi*3174*t)*(22.65/490)+sin(2*pi*6936*t)*(2.996/490)*0.18; disp(按任意键显示左声道简单合成声音波形！); disp(-); pause plot(t,y1); set(gca,FontName,Arial,FontSize,10); title(左声道简单合成声音波形); disp(按任意键听简单合成左声道信号三次！); disp(-); pause sound(y1,fs,bits); % 数据通过声卡转换为声音 pause sound(y1,fs,bits); pause sound(y1,fs,bits); k1=fft(y1,32768); disp(按任意键显示左声道简单合成声音的频域幅值！); disp(-);,24,音频信号的分析与合成程序（续4）,disp(-); pause plot(abs(k1); set(gca,FontName,Arial,FontSize,10); title(左声道简单合成声音的频域幅值); xlabel(采样点数); m1_r,i1_r=max(abs(k1); % 找出简单合成声音频域最大值 F1_r=i1_r/32768*10000; F1_rt=num2str(F1_r) Hz; text(3000,850,F1_rt); m2_r,i2_r=max(abs(k1(10000:12000); F2_r=(10000+i2_r)/32768*10000; % 找出简单合成声音频域次大值 F2_rt=num2str(F2_r) Hz; text(9000,60,F2_rt); end; %=,25,音频频谱分析仪程序,（车晴P199-audspec.m）,返回,y=audspec(action); % AUDIO SPECTRUM ANALYSER THEE 2001/10 %= global x t dt fs c1 c2 H_osc H2 fn y ystr Ktype if nargin1;action=initialized;end; if strcmp(action,initialized); % 示波器 clc; % 清除工作区 v=version; % MATLAB版本 fname,pname=uigetfile(*.wav,Open Wave File);% 打开对话框输入文件 file=pname,fname; % 文件名称和路径 if strcmp(v(1),5); % 读入声音文件（*.wav） x,fs,nbits=wavread(file); % x是采样数据，fs是采样速率（Hz） m,n=size(x); % nbits是每次采样比特数 file=file fs= int2str(fs) Hz, int2str(nbits) bits; if n=2;file=file ,stero; x=x(:); % 若是双声道，将数据合并为一列 %x=x(:,1); % 取左声道数据 %x=x(:,2); % 取右声道数据 end; elseif strcmp(v(1),4); x,fs=wavread(file); file=file fs= int2str(fs) Hz,; end;,26,音频频谱分析仪程序（续1）,% 增加菜单项（Zoom，Play，Spectrum），生成信号波形，回应（Callback）分别调用 % 函数audspec(Zoom)，audspec(Play)和 audspec(Spectrum)。 H_osc=figure(Name,OSCILLOSCOPE THEE 2001/8 blanks(6) file,. NumberTitle,off,Position,10 78 640 480); uimenu(Label,27,音频频谱分析仪程序（续2）,elseif strcmp(action,Play); % 播放声音文件 v=axis; % 返回当前轴坐标的设置 del=t(2)-t(1);del=del*.5; % v(1)和v(2)是X轴坐标的设置 n1=find(abs(t-v(1)del);n2=find(abs(t-v(2)del); if isempty(n2);n2=length(x);end; xx=x(n1:n2);sound(xx,fs); % 以采样速率fs输出声音文件数据 elseif strcmp(action,Zoom); % 放大 zoom(H_osc,xon); % 只能X轴放大 elseif strcmp(action,Spectrum); % 频谱分析仪 v=axis; % 返回当前轴坐标的设置 del=t(2)-t(1);del=del*.5; % v(1)和v(2)是X轴坐标的设置 n1=find(abs(t-v(1)del);n2=find(abs(t-v(2)del); if isempty(n2);n2=length(x);end;,28,音频频谱分析仪程序（续3）,% 增加菜单项（Log），生成信号频谱分析图形，回应（Callback）调用函数audspec(Log)。 %if H20;delete(H2);end; H2=figure(Name,SPECTRUM ANALYZER THEE 2001/8,. NumberTitle,off,Position,155 32 640 480); uimenu(Label,29,音频频谱分析仪程序（续4）,fn=f1*(0:N2-1);I=find(fn20);fn=fn(I);y=y(I); plot(fn,y,c2); % 画信号频谱分析图形 ystr=(dB);Ktype=-1; set(gca,Ygrid,on,YLim,-60 0,. % gca是返回当前轴的句柄并设置 FontName,Arial,FontSize,10); xlabel(frequency ( Hz);ylabel(ystr);zoom xon; elseif strcmp(action,Log); % 线性/对数坐标（频率） if Ktype=-1; semilogx(fn,y,c2); % 半对数坐标（X轴）画频谱分析图形 elseif Ktype=1; plot(fn,y,c2); end; Ktype=-Ktype; set(gca,Ygrid,on,YLim,-60 0,. FontName,Arial,FontSize,10); xlabel(frequency ( Hz);ylabel(ystr);zoom xon; end; %=,30,幅度调制程序,（车晴P203- am.m）,返回,function =am(action) % AM Modulate THEE 2001/10 % =am(action); %= global y x t fc fs f tstr Kw Kmod m H1 H2 H3 H42 H52 v=version;if nargin1;action=initialized;end; if strcmp(action,initialized); if strcmp(v(1),5); figure(Num,off,Units,pix,pos,. 5 29 792 530,Color,1 1 1); elseif strcmp(v(1),4); figure(num,off,units,pix,pos,5 29 792 530); whitebg(w); end; uicontrol(Style,Frame,Units,Normal,Position,. .82 0 .2 1,Back,.8 .8 .8); H1=uicontrol(Style,Popup,Units,Normal,Position,. .84 .9 .15 .05,String,str2mat(AM,AM-DSB,. AM-SSB Lower Band,AM-SSB Upper Band),. Back,1 1 1,Callback,am method;);,31,幅度调制程序（续1）,H2=uicontrol(Style,Popup,Units,Normal,Position,. .84 .8 .15 .05,String,str2mat(Sine,. Square,Sawtooth,Diric),. Back,1 1 1,Callback,am wave;); H3=uicontrol(Style,Popup,Units,Normal,Position,. .84 .7 .15 .05,String,str2mat(100%,. 90%,87.5%,70%,60%,50%),. Back,1 1 1,Callback,am mr;); H41=uicontrol(Style,Text,Units,Normal,Position,. .83 .595 .05 .04,String,f(Hz):,Back,.8 .8 .8); H42=uicontrol(Style,Edit,Units,Normal,Position,. .89 .6 .1 .04,Back,1 1 1,String,. 15.625e+03,Call,am f;); H51=uicontrol(Style,Text,Units,Normal,Position,. .822 .495 .06 .04,String,fc(Hz):,Back,.8 .8 .8); H52=uicontrol(Style,Edit,Units,Normal,Position,. .89 .5 .1 .04,Back,1 1 1,String,38e+06,Call,am fc;); H6=uicontrol(Style,Push,Units,Normal,Position,. .84 .4 .15 .05,String,RESTORE,Call,am restore;); H7=uicontrol(Style,Push,Units,Normal,Position,. .84 .3 .15 .05,String,DEMOD,Call,am demod;);,32,幅度调制程序（续2）,f=15.625e+03;fc=38e+06;m=1;Kw=1;Kmod=1; y,x,t,fc,fs,tstr=am1(f,fc,m,Kmod,Kw,1) elseif strcmp(action,method); Kmod=get(gco,Value); y,x,t,fc,fs,tstr=am1(f,fc,m,Kmod,Kw,1); if Kmod1;set(H3,value,1);end; elseif strcmp(action,wave); Kw=get(gco,Value); y,x,t,fc,fs,tstr=am1(f,fc,m,Kmod,Kw,1); elseif strcmp(action,mr); k=get(gco,Value);mt=1 .9 .875 .7 .6 .5;m=mt(k); if Kmod=1; y,x,t,fc,fs,tstr=am1(f,fc,m,Kmod,Kw,1); else; set(H3,value,1); end; elseif strcmp(action,f); f=get(gco,string);f=str2num(f); y,x,t,fc,fs,tstr=am1(f,fc,m,Kmod,Kw,1);,33,幅度调制程序（续3）,elseif strcmp(action,fc); fc=get(gco,string);fc=str2num(fc); y,x,t,fc,fs,tstr=am1(f,fc,m,Kmod,Kw,1); elseif strcmp(action,restore); f=15.625e+03;fc=38e+06;m=1;Kw=1;Kmod=1; y,x,t,fc,fs,tstr=am1(f,fc,m,Kmod,Kw,1); set(H42,string,15.625e+03);set(H52,string,38e+06); set(H1,value,1);set(H2,value,1);set(H3,value,1); elseif strcmp(action,demod); x1=demod(y,fc,fs,am);x1=x1-mean(x1); if strcmp(v(1),4); subplot(1,1,1);set(gca,vis,off); subplot(position,.065 .5811 .7 .3439); plot(t,y,b,t,x,r);zoom xon;xlabel(t( tstr ); end; subplot(position,.065 .11 .7 .3439); i=1:length(t)-8;plot(t(i),x1(i),b); xlabel(t( tstr );title(Demodulated Signal); end; %=,34,幅度调制程序（续4）,y,x,t,fc,fs,tstr=am1(f,fc,m,Kmod,