matlab在信号处理中的应用.ppt

第4章 MATLAB在信号处理中的应用,4.1 信号及其表示4.2 信号的基本运算 4.3 信号的能量和功率 4.4 线性时不变系统4.5 线性时不变系统的响应4.6 线性时不变系统的频率响应 4.7 傅里叶(Fourier)变换4.8 IIR数字滤波器的设计方法4.9 FIR数字滤波器设计,4.1 信号及其表示,连续时间信号的表示 在MATLAB中，对连续时间信号用采样点的数据来表示，当采样点很密时可看成连续信号。连续时间信号：时间变化连续。如y=x(t)离散时间信号(序列)：时间离散，如x(nT)=x(t)|t=nT.当T足够小时，可以认为是连续信号。例如：用MATLAB命令绘出 关于t的曲线，t的范围为030s,并以0.1递增。,t=0:0.1:30;%对时间变量赋值x=exp(-0.707*t).*sin(2/3.*t);%计算变量所对应的函数值plot(t,x);%绘制函数曲线grid;%图形窗口加网格xlabel(time(sec);%标注x轴,y轴ylabel(x(t);axis(-0.05,30,-0.05,0.35);%确定x轴，y轴显示范围,工具箱中的信号产生函数,t=0:0.0001:1;x1=sawtooth(2*pi*50*t,0.8)plot(t,x1)axis(0,0.2,-1,1),t=0:0.0001:0.05;x1=square(2*pi*50*t,50)plot(t,x1),例如：产生周期为0.02的三角波,产生频率为50Hz，占空比为50%的方波,t=0:0.001:2y=chirp(t,0,1,150)figure(1)plot(t,y)axis(0,0.5,0,1),产生线性调频信号：,t=0:0.01:1;y1=tripuls(t);y2=tripuls(t,0.6);subplot(211);plot(t,y1);subplot(212);plot(t,y2);,产生非周期三角波信号：,离散时间信号的表示,在MATLAB中，离散时间信号x(n)的表示：需用一个向量x表示序列幅值，用另一个等长的定位时间变量n，才能完整地表示一个序列。,例4-10 绘制离散时间信号的棒状图。其中x(-1)=-1,x(0)=1,x(1)=2,x(2)=1,x(3)=-1,x(4)=0。MATLAB源程序为：n=-3:5;%定位时间变量x=0,0,-1,1,2,1,-1,0,0;stem(n,x);grid;%绘制棒状图line(-3,5,0,0);%画x轴线xlabel(n);ylabel(xn)运行结果如图4.11所示。,图 4.11 离散时间信号图形,几种常用离散时间信号的表示,1单位脉冲序列,直接实现：x=zeros(1,N);x(1,n0)=1;,2单位阶跃序列,直接实现：n=ns:nf;x=(n-n0)=0;函数实现：,function x,n=stepseq(n0,ns,nf)n=ns:nf;x=(n-n0)=0;,3实指数序列,直接实现：n=ns:nf;x=a.n;,4复指数序列,直接实现：n=ns:nf;x=exp(sigema+jw)*n);,5正(余)弦序列,直接实现：n=ns:nf;x=cos(w*n+sita);,复数求实部：real(x)复数求虚部：imag(x)复数求幅度：abs(x)复数求相角：angel(x),复数运算的常用函数：,function fsxl(delta,omig,n1,n2)e=delta+omig.*j;n=n1:n2;x=exp(e*n);x_real=real(x);%生成实部序列x_imag=imag(x);%生成虚部序列x_magnitude=abs(x);%生成幅度序列x_phase=(180/pi)*angle(x);%生成相位序列,fsxl(0.1,pi/3,-10,10),4.2 信号的基本运算,信号的相加与相乘 y(n)=x1(n)+x2(n)y(n)=x1(n)x2(n)MATLAB实现：y=x1+x2;y=x1.*x2两信号相加和相乘的MATLAB实现：首先将两序列时间变量延拓到相同长度，然后再逐点相加相乘,序列移位与周期延拓运算,序列移位：y(n)=x(n-m)。MATLAB实现：y=x;ny=nx-m序列周期延拓：y(n)=x(n)M，MATLAB实现：ny=nxs:nxf；y=x(mod(ny,M)+1),n1=-5:4;%序列x1(n)的起始及终止位置n1s=-5;n1f=4;x1=2 3 1-1 3 4 2 1-5-3;%序列x1(n)不同时间的幅度n2=0:9;%序列x2(n)的起始及终止位置n2s=0;n2f=9;x2=1 1 1 1 1 1 1 1 1 1;%序列x2(n)不同时间的幅度ns=min(n1s,n2s);nf=max(n1f,n2f);%求新信号的时间起始终止位置n=ns:nf;y1=zeros(1,length(n);%延拓序列初始化y2=zeros(1,length(n);y1(find(n=n1s),N=24;M=8;m=3;n=0:N-1;x1=(0.8).n;%产生指数序列x2=(n=0),4.2.4 两序列的卷积运算,两序列卷积运算：,MATLAB实现：y=conv(x1,x2)。序列x1(n)和x2(n)必须长度有限。,4.2.5 两序列的相关运算,两序列相关运算：,MATLAB实现：y=xcorr(x1,x2)。,4.2.3 序列翻褶与序列累加运算,序列翻褶：y(n)=x(-n)。MATLAB可实现:y=fliplr(x),序列累加的数学描述为：,MATLAB实现：y=cumsum(x),例如：求序列,n=0:10;x=3*exp(-0.2*n);y=fliplr(x);n1=-fliplr(n);n2=-fliplr(n-3);%在指定位置m=3处的时间序列翻褶s=cumsum(x);subplot(221);stem(n,x);xlabel(n);ylabel(x(n);subplot(223);stem(n1,y);xlabel(n);ylabel(y(n)=x(-n);subplot(222);stem(n2,y);xlabel(n);ylabel(y(n)=x(-3+n);subplot(224);stem(n,s);xlabel(n);ylabel(s(n);,4.3 信号的能量和功率,1.信号能量,数字定义：,MATLAB实现:E=sum(x.*conj(x);或 E=sum(abs(x).2);,数字定义：,2.信号功率,MATLAB实现:P=sum(x.*conj(x)/N;或 E=sum(abs(x).2)/N;,非周期三角波信号能量的MATLAB计算：,dt=0.0001;t=0:dt:1;x=tripuls(t);E=sum(abs(x).2*dt),dt=0.001;t=0:dt:2*pi;x=cos(t);P=sum(abs(x).2*dt)./(2*pi),E=0.1667,P=0.5001,余弦信号的平均功率MATLAB计算：,4.4 线性时不变系统,一个信号处理系统就是将输入信号变换成输出信号的运算过程，如图所示。在此过程中，输出的信号称为系统对输入信号作出的响应，输入信号称为系统的激励信号。,当一个系统具有可加性和齐次性，则该系统称为线性系统。如果系统响应与激励加于系统的时刻无关，则称该系统为时不变系统。,4.4 线性时不变系统,4.4.1 系统的描述,1常系数线性微分/差分方程,2系统传递函数,连续系统：,离散系统：,连续系统：,离散系统：,注意：在MATLAB中，传递函数由分子、分母两个多项式的系数来表示，系数为降幂排列。例如：,系统传递函数的MATLAB表示：,可以表示为：num=1,0.2,1;den=1,0.5,1,3零极点增益模型,连续系统：,离散系统：,4极点留数模型,离散系统：,连续系统：,5二次分式模型,连续系统：,离散系统：,6状态空间模型,连续系统：,离散系统：,系统模型的MATLAB表示,在MATLAB中，用sos、ss、tf、zp分别表示二次分式模型、状态空间模型、传递函数模型和零极点增益模型。其中sos表示二次分式，g为比例系数，sos为L6的矩阵，即,4.4.2 系统模型的转换函数,1ss2tf函数格式：num,den=ss2tf(A,B,C,D,iu)功能：将指定输入量iu的线性系统(A,B,C,D)状态空间模型转换为传递函数模型num,den。,2zp2tf函数格式：num,den=zp2tf(z,p,k)功能：将给定系统的零极点增益模型转换为传递函数模型，z、p、k分别为零点列向量、极点列向量和增益系数。,3tf2sos函数格式：sos,g=tf2sos(num,den)功能：将给定系统的传递函数模型num,den转换为二次分式模型sos,g为增益系数。,线性系统模型的变换函数,在命令窗口输入：z=3;p=1,2;k=2;A,B,C,D=zp2ss(z,p,k)屏幕显示结果：A=3.0000-1.4142 1.4142 0B=1 0C=2.0000-4.2426D=0,例4-18 将系统 转换为状态空间模型A,B,C,D,例4-19 求离散时间系统,的零、极点向量和增益系数。在命令窗口输入：num=2,3;den=1,0.4,1;num,den=eqtflength(num,den);%使长度相等 z,p,k=tf2zp(num,den)屏幕显示为z=0-1.5000p=-0.2000+0.9798i-0.2000-0.9798ik=2,4.4.3 系统互联与系统结构,MATLAB实现函数series()格式：A,B,C,D=series(A1,B1,C1,D1,A2,B2,C2,D2)或 num,den=series(num1,den1,num2,den2),将系统1、系统2级联，可得到级联连接的传递函数形式为：,1.系统的级联,MATLAB实现函数parallel()格式：A,B,C,D=parallel(A1,B1,C1,D1,A2,B2,C2,D2)或 num,den=parallel(num1,den1,num2,den2),2.系统的并联,将系统1、系统2并联，可得到并联连接的传递函数形式为：,3.两个系统的反馈连接函数feedback格式：A,B,C,D=feedback(A1,B1,C1,D1,A2,B2,C2,D2,sign)或 num,den=feedback(num1,den1,num2,den2,sign)将系统1和系统2进行反馈连接，sign表示反馈方式(默认值为-1)；当sig=+1时表示正反馈；当sig=-1时表示负反馈。,例4-20 求两个单输入单输出子系统,的级联、并联和反馈后系统的传递函数。MATLAB源程序为：num1=1;den1=1,1;%系统1num2=2;den2=1,2;%系统2nums,dens=series(num1,den1,num2,den2)%实现两个系统级联nump,denp=parallel(num1,den1,num2,den2)%实现两个系统并联 numf,denf=feedback(num1,den1,num2,den2)%实现两个系统反馈程序运行结果为：nums=0 0 2；dens=1 3 2nump=0 3 4；denp=1 3 2numf=0 1 2；denf=1 3 4因此，各系统的传递函数分别为：,在实际应用中，也可以把一个复杂的线性时不变（LTI）系统分解为几个简单系统的组合结构，即直接型结构、级联型结构和并联型结构。Matlab所提供的系统模型变换函数，实质就是给出了这几种系统结构的互换关系。系统的传递函数对应于系统的直接型结构，二次分式（sos）模型对应级联型结构,系统传递函数的部分分式（residue或residuez)形式对应于并联型结构。,例4-20 已知FIR数字滤波器的传递函数为：求出其级联型结构和格型结构,clear;b=2,13/12,5/4,2/3;a=1;%设定参数 fprintf(级联型结构系数：);sos,g=tf2sos(b,a)%直接型到级联型转换fprintf(格型结构系数（反射系数）：);K=tf2latc(b,a)%直接型到格型转换,MATLAB源程序：,级联型结构系数：sos=1.0000 0.5360 0 1.0000 0 0 1.0000 0.0057 0.6219 1.0000 0 0g=2格型结构系数（反射系数）：K=0.2500 0.5000 0.3333,b=1,-3,11,-27,18;a=16,12,2,-4,-1;disp(级联型结构系数：)sos,g=tf2sos(b,a)disp(并联型结构系数：)R,P,K=residuez(b,a),MATLAB源程序：,级联型结构系数：sos=1.0000-3.0000 2.0000 1.0000-0.2500-0.1250 1.0000 0.0000 9.0000 1.0000 1.0000 0.5000g=0.0625并联型结构系数：R=-5.0250-1.0750i-5.0250+1.0750i 0.9250 27.1875,P=-0.5000+0.5000i-0.5000-0.5000i 0.5000-0.2500 K=-18,4.5 线性时不变系统的响应,线性时不变系统的时域响应,1连续LTI系统的响应,2离散LTI系统的响应,用MATLAB中的卷积函数conv()来实现。,用MATLAB中的卷积函数conv()来实现。,dt=input(输入时间间隔dt)x=ones(1,fix(10/dt);h=exp(-0.1*0:fix(10/dt)*dt);y=conv(x,h);t=dt*(1:length(y)-1);plot(t,y);grid;,x=ones(1,10);n=0:14;h=0.5.n;y=conv(x,h);stem(y);xlabel(n);ylabel(yn);,nx=-5:4;x=ones(1,10);nh=0:14;h=0.5.n;y=conv(x,h);n0=nx(1)+nh(1);%求卷积序列y起始时间位置N=length(nx)+length(nh)-2;%求卷积序列y序列长度ny=n0:n0+N;%求卷积序列y的时间向量subplot(221);stem(nx,x);title(xn);xlabel(n);xlabel(xn);subplot(222);stem(nh,h);title(hn);xlabel(n);xlabel(hn);subplot(223);stem(ny,y);title(xn*hn);xlabel(n);xlabel(yn);h=get(gca,position);h(3)=2.5*h(3)set(gca,position,h),改写为：,格式：y,x=lsim(a,b,c,d,u,t)功能：返回连续LTI系统,(2)对任意输入的离散LTI系统响应函数dlsim()格式：y,x=dlsim(a,b,c,d,u)功能：返回离散LTI系统,对任意输入时系统的输出响应y和状态记录x，其中u给出每个输入的时序列，一般情况下u为一个矩阵；t用于指定仿真的时间轴，它应为等间隔。,对输入序列u的响应y和状态记录x。,3时域响应函数(1)对任意输入的连续LTI系统响应函数lsim(),num=2,5,1;den=1,2,3;t=0:0.1:10peiod=4;u=(rem(t,peiod)=peiod/2)lsim(num,den,u,t);title(方波响应);,num=2,-3.4,5.5;den=1,-1.2,0.8；u=randn(1,100);dlsim(num,den,u);title(随机噪声响应),4.5.2 LTI系统的单位冲激响应,1.求连续LTI系统的单位冲激响应函数impulse()格式：Y,T=impulse(sys)或impulse(sys)功能：返回系统的响应Y和时间向量T，自动选择仿真的时间范围。其中sys可为系统传递函数、零极增益模型或状态空间模型。,2.求离散系统的单位冲激响应函数dimpulse()格式：y,x=dimpulse(num,den)功能：返回多项式传递函数,的单位冲激响应y向量和时间状态历史记录x向量。,a=-0.55,-0.78;0.78,0;b=1;0;c=1.96,6.45;d=0;impulse(a,b,c,d);title(LTI系统的冲激响应),num=2,-3.5,1.5;den=1,-1.7,0.3;dimpulse(num,den);title(离散LTI系统的冲激响应),4.5.3 时域响应的其它函数1.求连续LTI系统的零输入响应函数initial()格式：y,t,x=initial(a,b,c,d,x0)功能：计算出连续时间LTI系统由于初始状态x0所引起的零输入响应y。其中x为状态记录，t为仿真所用的采样时间向量。,2.求离散系统的零输入响应函数dinitial()格式：y,x,n=dinitial(a,b,c,d,x0)功能：计算离散时间LTI系统由初始状态x0所引起的零输入响应y和状态响应响应x，取样点数由函数自动选取。n为仿真所用的点数。,3.求连续系统的单位阶跃响应函数step()格式：Y,T=step(sys)功能：返回系统的单位阶跃响应Y和仿真所用的时间向量T，自动选择仿真的时间范围。其中sys可为系统传递函数(TF)、零极增益模型(ZPK)或状态空间模型(SS)。4.求离散系统的单位阶跃响应函数dstep()格式：y,x=dstep(num,den)功能：返回多项式传递函数G(z)=num(z)/den(z)表示的系统单位阶跃响应。,4.6线性时不变系统的频率响应,4.6线性时不变系统的频率响应,1求模拟滤波器Ha(s)的频率响应函数freqs()格式一：Hfreqs(B,A,W)功能：计算由向量W(rad/s)指定的频率点上模拟滤器系统函数Ha(s)的频率响应Ha(j)，结果存于H向量中。格式二：H,wfreqs(B,A,M)功能：计算出M个频率点上的频率响应存于向量H中,M个频率存放在向量w中。freqs自动将这M个频率点设置在适当的频率范围。默认w和M时,freqs自动选取200个频率点计算。,求该模拟滤波器的频率响应。MATLAB源程序如下：B=1;A=1 2.6131 3.4142 2.6131 1;W=0:0.1:2*pi*5;freqs(B,A,W),例4-31 已知某模拟滤波器的系统函数,图4.30 模拟滤波器的频率响应,例4-32 已知某滤波器的系统函数为,2求数字滤波器H(z)的频率响应函数freqz()格式一：H=freqz(B,A,W)功能：计算由向量W(rad)指定的数字频率点上(通常指在H(z)的频率响应H(ejw)。格式二：H,W=freqz(B,A,N)格式三：H,W=freqz(B,A,N,whole),求该滤波器的频率响应。MATLAB源程序为：B=1 0 0 0 0 0 0 0 1;A=1;freqz(B,A),图4.31滤波器幅度和相位曲线,该程序运行所绘出的幅频与相频性曲线如图4.31所示。,例：,已知一系统传递函数为：,求系统的单位冲激响应h(n),以及h(n)的幅频和相频响应。,N=128;x=1,zeros(1,N-1);%产生单位冲激函数num=0.01-0.03 0.06-0.03 0.01;den=1 2.5 3 1.8 0.5;y=dlsim(num,den,x);%计算单位冲激响应figure(1);n=1:N;stem(n,y);grid on;title(单位冲激响应);figure(2)freqz(num,den,N);%绘制幅频和相频响应曲线grid on;,dimpulse(num,den),单位冲激响应函数，不带输出变量时，直接绘出单位冲激响应,单位冲激响应,幅频和相频响应,冲激响应,3滤波函数filter从频域角度，无论是连续时间LTI系统还是离散时间LTI系统，系统对输入信号的响应实质上就是对输入信号频谱进行不同选择处理的过程，这个过程称为滤波。格式：y=filter(B,A,x)功能：对向量x中的数据进行滤波处理，即差分方程求解，产生输出序列向量y。B和A分别为数字滤波器系统函数H(z)的分子和分母多项式系数向量。,例4-33 设系统差分方程为,求该系统对信号,的响应。,MATLAB源程序为：B=1;A=1,-0.8;n=0:31;x=0.8.n;y=filter(B,A,x);subplot(2,1,1);stem(x)subplot(2,1,2);stem(y),图4.32系统对信号的响应,该程序运行所得结果如图4.32所示,傅里叶变换：是建立以时间为自变量的“信号”与以频率为自变量的“频谱函数”之间的某种对应关系。遵循以下原则：时域连续周期离散非周期,非周期,离散,周期,连续,频域,4.7傅里叶(Fourier)变换,4.7傅里叶(Fourier)变换,连续时间、连续频率傅里叶变换,正变换：,逆变换：,这说明求和的问题可以用x(t)行向量乘以 列向量来实现。式中 是 t的增量，在程序中用 dt表示.,例：求如图矩形脉冲信号f(t)在-4040rad/s区间的频谱。,clc;clear all;N=input(取时间分隔的点数N=);dt=10/N;t=1:N*dt;%给出时间分割f=ones(1,N/2),zeros(1,N/2);%给出信号(方波)wf=input(需求的频谱宽度wf=);Nf=input(需求的频谱点数Nf=);w1=linspace(0,wf,Nf);F1=f*exp(-j*t*w1)*dt;%求傅里叶变换w=-fliplr(w1),w1(2:Nf);%补上负频率F=fliplr(F1),F1(2:Nf);%补上负频率对应的频谱subplot(121)plot(t,f,linewidth,1.5)grid on;axis(0,10,0,1.1);subplot(122)plot(w,abs(F),linewidth,1.5);grid on;,取时间分隔的点数N=64需求的频谱宽度wf=40需求的频谱点数Nf=256,取时间分隔的点数N=256需求的频谱宽度wf=40需求的频谱点数Nf=64,采样较密,采样稀，有频率泄漏,4.7.2 连续时间、离散频率傅里叶级数,正变换：,逆变换：,式中,为离散频率相邻两谱线间的角频率间隔，k为谐波序号。,4.7.3 离散时间、离散频率离散傅里叶级数,正变换：,逆变换：,离散时间、离散频率离散傅里叶变换(DFT),正变换：,逆变换：,4.7.4 时间离散、连续频率序列傅里叶变换,正变换：,逆变换：,注意，对于序列傅里叶变换，如果x(n)为无限长，那么就不能用MATLAB直接利用上式计算，只可以用它对表达式 在 频率点上求值，再画出它的幅度和相位。如果x(n)为有限长，就可直接用MATLAB计算,离散时间、离散频率离散傅里叶变换(DFT),在MATLAB中，旋转因子矩阵可以由dftmtx函数实现。,dftmtx函数:调用格式：W=dftmtx(n)功能：函数调用后，返回一个n*n的离散傅里叶变换DFT的旋转因子矩阵W。这样对于给定长度为n 的行向量信号x(n),利用Xk=x*W就可以获得x(n)的傅里叶变换X(k)。其实，旋转因子矩阵也可以自己通过表达式构造。,例如：,编写函数，首先利用MATLAB计算序列x(n)的DFT和IDFT函数.,function xn=idft(xk)N=length(xk);n=0:N-1;k=0:N-1;WN=exp(-j*2*pi/N);WNnk=WN.(-n*k);xn=xk*WNnk/N;,IDFT:,function xk=dft(xn)N=length(xn);WNnk=dftmtx(N);xk=xn*WNnk;,DFT:,WNnk=conj(dftmtx(N),DFT参数对频率分辨率的影响：频率分辨率就是，在频率轴上能分辨出的两个频率点的最小间隔,要能区分两个频率点f1,f2,，有效数据长度N必须满足：,可以对序列进行补零操作，以增加数据的长度，这样做虽然不能提高频率分辨率，但补零后对原来的X(k)起到了插值作用，可以平滑频谱的包络。但是，如果增加有效数据的长度，则可以提高频率分辨率。,利用上面编写的DFT函数，求,的频谱特性，其中取fs=20Hz,N=128;n=0:N-1;fs=20;xn=sin(3.1*2*pi*n/fs)+cos(3*2*pi*n/fs);,xk=dft(xn);k=(0:N/2-1)*fs/Nsubplot(211);plot(n,xn);grid on;title(x(n)(0=n128)subplot(212);plot(k,abs(xk(1:N/2);grid on;,N=128;n=0:N-1;fs=20;xn=sin(3.1*2*pi*n/fs)+cos(3*2*pi*n/fs);N1=256;n1=0:N1-1;xn1=xn,zeros(1,N1-N);xk1=dft(xn1);m=(0:N1/2-1)*fs/N1subplot(211);plot(n1,xn1);grid on;title(x(n)补零后);subplot(212);plot(m,abs(xk1(1:N1/2);grid on;title(x(n)补零后),N=402;n=0:N-1;fs=20;xn=sin(3.1*2*pi*n/fs)+cos(3*2*pi*n/fs);xk=dft(xn);k=(0:N/2-1)*fs/Nsubplot(211);plot(n,xn);grid on;title(x(n)(0=n401)subplot(212);plot(k,abs(xk(1:N/2);grid on;,N=128;n=0:N-1;fs=20;xn=sin(3.1*2*pi*n/fs)+cos(3*2*pi*n/fs);xk=dft(xn);subplot(211)plot(abs(xk);grid on;subplot(212)k=(0:N-1)*fs/N;plot(k,abs(xk);grid on;,1一维快速正傅里叶变换函数fft格式：X=fft(x,N)功能：采用FFT算法计算序列向量x的N点DFT变换，当N缺省时，fft函数自动按x的长度计算DFT。当N为2整数次幂时，fft按基-2算法计算，否则用混合算法。2一维快速逆傅里叶变换函数ifft格式：x=ifft(X,N)功能：采用FFT算法计算序列向量X的N点IDFT变换。3.将零频分量移至频谱中心的函数格式：Y=fftshift(X)功能：用来重新排列X=fft(x)的输出，把X 的左右两半进行交换，从而将零频分量移至频谱中心。,fs=1000;%采样频率t=0:1/fs:0.5;%时间轴取值范围及间隔x=cos(160*pi*t);N=512;%FFT点数subplot(311)；plot(t,x);grid on;y=fft(x,N);%N点FFTx(t)xk=fftshift(abs(y);%求幅度谱并将零频分量移至频谱中心位置 f=(-N/2:N/2-1)*fs/512;%频率范围subplot(312)plot(f,xk);%做频谱图grid on;%用IFFT恢复原始信号z=real(ifft(y,N);ti=0:length(z)-1/fs;subplot(313)plot(ti,z);grid on;,例：求,的FFT序列和IFFT序列,例4-36 用快速傅里叶变换FFT计算下面两个序列的卷积。,图4.35 快速卷积框图,MATLAB程序(部分)：%线性卷积xn=sin(0.4*1:15);%对序列x(n)赋值,M=15hn=0.9.(1:20);%对序列h(n)赋值,N=20ticyn=conv(xn,hn);%直接调用函数conv计算卷积toc%园周卷积L=pow2(nextpow2(M+N-1);tic Xk=fft(xn,L);Hk=fft(hn,L);Yk=Xk.*Hk;yn=ifft(Yk,L);toc,图4.36 x(n),h(n)及其线性卷积波形,Elapsed time is 0.005074 seconds.Elapsed time is 0.000072 seconds.,练 习,求x(t)在有噪声和无噪声条件下的FFT.(提示：噪声由randn函数生成),1.,4.8 IIR数字滤波器的设计方法,1.数字滤波器的频率响应函数,幅度响应：,相位响应：,图4.37 理想低通、高通、带通、带阻数字滤波器幅度特性,2.滤波器的技术指标 幅度响应指标、相位响应指标,图4.38 数字低通滤波器的幅度特性,通带要求：,阻带要求：,通带最大衰减：,阻带最小衰减：,冲激响应不变法,2.MATLAB信号处理工箱中的专用函数impinvar()：格式：BZ,AZ=impinvar(B,A,Fs)功能：把具有B,A模拟滤波器传递函数模型转换成采样频率为Fs(Hz)的数字滤波器的传递函数模型BZ,AZ。采样频率Fs的默认值为Fs=1。,1.冲激响应不变法设计IIR数字滤波器的基本原理：,例4-37 MATLAB源程序如下：num=1;%模拟滤波器系统函数的分子den=1,sqrt(5),2,sqrt(2),1;%模拟滤波器系统函数的分母num1,den1=impinvar(num,den)%求数字低通滤波器的系统函数程序的执行结果如下：num1=-0.0000 0.0942 0.2158 0.0311den1=1.0000-2.0032 1.9982-0.7612 0.1069,MATLAB信号处理工具箱中的专用双线性变换函数bilinear()格式：numd,dendbilinear(num,den,Fs)功能：把模拟滤波器的传递函数模型转换成数字滤波器的传递函数模型。,双线性变换法,双线性变换利用频率变换关系：,例4-38 MATLAB源程序如下：num=1;%模拟滤波器系统函数的分子 den=1,sqrt(3),sqrt(2),1;%模拟滤波器系统函数的分母 num1,den1=bilinear(num,den,1)%求数字滤波器的传递函数运算的结果如下：num1=0.0533 0.1599 0.1599 0.0533den1=1.0000-1.3382 0.9193-0.1546,4.8.3 IIR数字滤波器的频率变换设计法,1.IIR数字滤波器的频率变换设计法的基本原理 根据滤波器设计要求，设计模拟原型低通滤波器，然后进行频率变换，将其转换为相应的模拟滤波器(高通、带通等)，最后利用冲激响应不变法或双线性变换法，将模拟滤波器数字化成相应的数字滤波器。,图4.39 IIR数字滤波器MATLAB设计步骤流程图,1MATLAB的典型设计,利用在MATLAB设计IIR数字滤波器可分以下几步来实现(1)按一定规则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标；(2)根据转换后的技术指标使用滤波器阶数函数，确定滤波器的最小阶数N和截止频率Wc；(3)利用最小阶数N产生模拟低通滤波原型；(4)利用截止频率Wc把模拟低通滤波器原型转换成模拟低通、高通、带通或带阻滤波器；(5)利用冲激响应不变法或双线性不变法把模拟滤波器转换成数字滤波器。,例4-39 设计一个数字信号处理系统，它的采样率为Fs100Hz，希望在该系统中设计一个Butterworth型高通数字滤波器，使其通带中允许的最小衰减为0.5dB，阻带内的最小衰减为40dB，通带上限临界频率为30Hz，阻带下限临界频率为40Hz。,MATLAB源程序设计如下：%把数字滤波器的频率特征转换成模拟滤波器的频率特征 wp=30*2*pi;ws=40*2*pi;rp=0.5;rs=40;Fs=100;N,Wc=buttord(wp,ws,rp,rs,s);%选择滤波器的最小阶数 Z,P,K=buttap(N);%创建Butterworth低通滤波器原型 A,B,C,D=zp2ss(Z,P,K);%零极点增益模型转换为状态空间模型 AT,BT,CT,DT=lp2hp(A,B,C,D,Wc);%实现低通向高通的转变 num1,den1=ss2tf(AT,BT,CT,DT);%状态空间模型转换为传递函数模型%运用双线性变换法把模拟滤波器转换成数字滤波器 num2,den2=bilinear(num1,den1,100);H,W=freqz(num2,den2);%求频率响应 plot(W*Fs/(2*pi),abs(H);grid;%绘出频率响应曲线 xlabel(频率/Hz);ylabel(幅值)程序运行结果如图4.40所示。,2MATLAB的直接设计,图4.39 IIR数字滤波器MATLAB设计步骤流程图,例4-41 试设计一个带阻IIR数字滤波器，其具体的要求是：通带的截止频率：wp1650Hz、wp2850Hz；阻带的截止频率：ws1700Hz、ws2800Hz；通带内的最大衰减为rp0.1dB；阻带内的最小衰减为rs50dB；采样频率为Fs2000Hz。MATLAB源程序设计如下：wp1=650;wp2=850;ws1=700;ws2=800;rp=0.1;rs=50;Fs=2000;wp=wp1,wp2/(Fs/2);ws=ws1,ws2/(Fs/2);%利用Nyquist频率频率归一化 N,wc=ellipord(wp,ws,rp,rs,z);%求滤波器阶数 num,den=ellip(N,rp,rs,wc,stop);%求滤波器传递函数 H,W=freqz(num,den);%绘出频率响应曲线 plot(W*Fs/(2*pi),abs(H);grid;xlabel(频率/Hz);ylabel(幅值)该程序运行后的幅频响应曲线如图4.42所示。,4.9 FIR数字滤波器设计,格式：w=boxcar(M)功能：返回M点矩形窗序列。MATLAB信号处理工具箱中的窗函数法设计FIR数字滤波器的专用命令fir1()。格式：Bfir1(N,wc)功能：设计一个具有线性相位的N阶(N点)的低通FIR数字滤波器，返回的向量B为滤波器的系数(单位冲激响应序列)，其长度为N+1。,窗函数设计法窗函数设计的基本原理：h(n)=w(n)hd(n)w(n)为窗函数，hd(n)理想数字滤波器的单位冲激响应。在MATLAB信号处理工具箱中为用户提供了Boxcar(矩形)、Bartlet(巴特利特)、Hanning(汉宁)等窗函数，这些窗函数的调用格式相同。,FIR数字滤波器的单位冲激响应h(n)满足偶(奇)对称 h(n)=h(N-n-1)或 h(n)=-h(N-n-1)FIR数字滤波器具有线性相位：,或,例4-43 用矩形窗设计线性相位FIR低通滤波器。该滤波器的通带截止频率wc=pi/4，单位脉冲响h(n)的长度M=21。并绘出h(n)及其幅度响应特性曲线。MATLAB源程序为：M=21;wc=pi/4;%理想低通滤波器参数n=0:M-1;r=(M-1)/2;nr=n-r+eps*(n-r)=0);hdn=sin(wc*nr)/pi./nr;%计算理想低通单位脉冲响应hd(n)if rem(M,2)=0,hdn(r+1)=wc/pi;end%M为奇数时，处理n=r点的0/0型wn1=boxcar(M);%矩形窗hn1=hdn.*wn1;%加窗subplot(2,1,1);stem(n,hn1,.);line(0,20,0,0);xlabel(n),ylabel(h(n),title(矩形窗设计的h(n);hw1=fft(hn1,512);w1=2*0:511/512;%求频谱subplot(2,1,2),plot(w1,20*log10(abs(hw1)xlabel(w/pi),ylabe