matlab在信号处理中的应用课件

第4章MATLAB在信号处理中的应用4.1信号及其表示4.2信号的基本运算4.3信号的能量和功率4.4线性时不变系统4.5线性时不变系统的响应4.6线性时不变系统的频率响应

4.7傅里叶(Fourier)变换4.8IIR数字滤波器的设计方法4.9FIR数字滤波器设计1可编辑课件PPT4.1信号及其表示4.1.1连续时间信号的表示

在MATLAB中，对连续时间信号用采样点的数据来表示，当采样点很密时可看成连续信号。连续时间信号：时间变化连续。如y=x(t)离散时间信号(序列)：时间离散，如x(nT)=x(t)|t=nT.当T足够小时，可以认为是连续信号。例如：用MATLAB命令绘出关于t的曲线，t的范围为0~30s,并以0.1递增。2可编辑课件PPTt=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轴显示范围3可编辑课件PPT函数名格式功能sawtoothSawtooth（t,width）产生锯齿波或三角波信号（周期为2，幅值-1~1）squareSquare(t,duty)产生方波信号（周期为2，幅值-1~1）sincSinc(x)产生sinc函数波形chirpchirp(t,f0,t1,f1,’method’,phi)产生线性调频扫频信号gauspulsgauspuls（T,FC,BW,BWR）产生高斯正弦脉冲信号vcovoc(x,fc,fs)电压控制振荡器pulstranpulstran(t,d,’func’)产生冲激串rectpulerectpule(t,w)产生非周期的方波信号tripulstripuls(t,w)产生非周期的三角波信号diricdiric(x,n)产生Dirichlet或周期sinc函数gmonopulsgmonopuls(t,fc)产生高斯单脉冲信号4.1.2工具箱中的信号产生函数4可编辑课件PPTt=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%的方波5可编辑课件PPTt=0:0.001:2y=chirp(t,0,1,150)figure(1)plot(t,y)axis([0,0.5,0,1])产生线性调频信号：6可编辑课件PPTt=0:0.01:1;y1=tripuls(t);y2=tripuls(t,0.6);subplot(211);plot(t,y1);subplot(212);plot(t,y2);产生非周期三角波信号：7可编辑课件PPT4.1.3离散时间信号的表示在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('x[n]')运行结果如图4.11所示。图4.11离散时间信号图形8可编辑课件PPT4.1.4几种常用离散时间信号的表示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];9可编辑课件PPT3．实指数序列直接实现：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)复数运算的常用函数：10可编辑课件PPTfunctionfsxl(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)11可编辑课件PPT4.2信号的基本运算4.2.1信号的相加与相乘

y(n)=x1(n)+x2(n)y(n)=x1(n)×x2(n)MATLAB实现：y=x1+x2;y=x1.*x2两信号相加和相乘的MATLAB实现：首先将两序列时间变量延拓到相同长度，然后再逐点相加相乘4.2.2序列移位与周期延拓运算序列移位：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)12可编辑课件PPTn1=[-5:4];%序列x1(n)的起始及终止位置n1s=-5;n1f=4;x1=[231-13421-5-3];%序列x1(n)不同时间的幅度n2=[0:9];%序列x2(n)的起始及终止位置n2s=0;n2f=9;x2=[1111111111];%序列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<=n1f)==1))=x1;%给延拓序列y1赋值x1y2(find((n>=n2s)&(n<=n2f)==1))=x2;%给延拓序列y2赋值x2ya=y1+y2;%逐点相加yp=y1.*y2;%逐点相乘subplot(411);stem(n,y1,'.');subplot(412);stem(n,y2,'.');subplot(413);stem(n,ya,'.');subplot(414);stem(n,yp,'.');13可编辑课件PPT14可编辑课件PPTN=24;M=8;m=3;n=0:N-1;x1=(0.8).^n;%产生指数序列x2=[(n>=0)&(n<=M)];%形矩形序列RM(n)x=x1.*x2;%截取操作形成新序列x(n)xm=zeros(1,N);fork=m+1:m+Mxm(k)=x(k-m);%产生移位序列x(n-3)endxc=x(mod(n,M)+1);%产生x(n)的周期延拓x((n))8xcm=x(mod(n-m,M)+1);%产生x(n)的周期延拓x((n-3))8subplot(411);stem(n,x,'.');ylabel('x(n)');subplot(412);stem(n,xm,'.');ylabel('x(n-3)');subplot(413);stem(n,xc,'.');ylabel('x((n))_8');subplot(414);stem(n,xcm,'.');ylabel('x((n-3))_8');15可编辑课件PPT16可编辑课件PPT4.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)17可编辑课件PPT例如：求序列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)');18可编辑课件PPT4.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;19可编辑课件PPT非周期三角波信号能量的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.1667P=0.5001余弦信号的平均功率MATLAB计算：20可编辑课件PPT4.4线性时不变系统一个信号处理系统就是将输入信号变换成输出信号的运算过程，如图所示。在此过程中，输出的信号称为系统对输入信号作出的响应，输入信号称为系统的激励信号。当一个系统具有可加性和齐次性，则该系统称为线性系统。如果系统响应与激励加于系统的时刻无关，则称该系统为时不变系统。21可编辑课件PPT4.4线性时不变系统4.4.1系统的描述1．常系数线性微分/差分方程2．系统传递函数连续系统：

离散系统：

连续系统：

离散系统：

22可编辑课件PPT注意：在MATLAB中，传递函数由分子、分母两个多项式的系数来表示，系数为降幂排列。例如：系统传递函数的MATLAB表示：可以表示为：num=[1,0.2,1];den=[1,0.5,1]3．零－极点增益模型连续系统：

离散系统：

23可编辑课件PPT4．极点留数模型离散系统：

连续系统：

5．二次分式模型连续系统：

离散系统：

6．状态空间模型连续系统：

离散系统：

24可编辑课件PPT系统模型的MATLAB表示系统模型表示名称参数传递函数模型tfnum---传递函数的分子系数向量den---传递函数的分母系数向量零极点增益模型zpz---系统零点向量；p---系统极点向量k---增益系数二次分式模型sossos---二次分式的系数矩阵；g---比例系数状态空间模型ss[A,B,C,D]---状态空间模型参数部分分式模型residuezR---分子向量；P---分母系数向量；K---常数在MATLAB中，用sos、ss、tf、zp分别表示二次分式模型、状态空间模型、传递函数模型和零－极点增益模型。其中sos表示二次分式，g为比例系数，sos为L×6的矩阵，即

25可编辑课件PPT4.4.2系统模型的转换函数1．ss2tf函数

格式：[num,den]=ss2tf(A,B,C,D,iu)功能：将指定输入量iu的线性系统(A,B,C,D)状态空间模型转换为传递函数模型[num,den]。2．zp2tf函数格式：[num,den]=zp2tf(z,p,k)功能：将给定系统的零－极点增益模型转换为传递函数模型，z、p、k分别为零点列向量、极点列向量和增益系数。3．tf2sos函数格式：[sos,g]=tf2sos

(num,den)功能：将给定系统的传递函数模型[num,den]转换为二次分式模型sos,g为增益系数。26可编辑课件PPT线性系统模型的变换函数函数名功能说明函数名功能说明ss2tf状态空间模型转换为传递函数模型zp2tf零－极点增益模型转换为传递函数模型ss2zp状态空间模型转换为零－极点增益模型zp2ss零－极点增益模型转换为状态空间模型ss2sos状态空间模型转换为二次分式模型zp2sos零－极点增益模型转换为二次分式模型tf2ss传递函数模型转换为状态空间模型sos2tf二次分式模型转换为传递函数模型tf2zp传递函数模型转换为零－极点增益模型sos2zp二次分式模型转换为零－极点增益模型tf2sos传递函数模型转换为二次分式模型sos2ss二次分式模型转换为状态空间模型27可编辑课件PPT在命令窗口输入：>>z=[3];p=[1,2];k=2;>>[A,B,C,D]=zp2ss(z,p,k)屏幕显示结果：A=3.0000-1.41421.41420B=10C=2.0000-4.2426D=0[例4-18]将系统转换为状态空间模型[A,B,C,D]28可编辑课件PPT[例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=229可编辑课件PPT4.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.系统的级联30可编辑课件PPTMATLAB实现函数parallel()格式：[A,B,C,D]=parallel(A1,B1,C1,D1,A2,B2,C2,D2)或

[num,den]=parallel(num1,den1,num2,den2)2.系统的并联将系统1、系统2并联，可得到并联连接的传递函数形式为：31可编辑课件PPT3.两个系统的反馈连接函数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时表示负反馈。32可编辑课件PPT[例4-20]求两个单输入单输出子系统的级联、并联和反馈后系统的传递函数。MATLAB源程序为：

num1=1;den1=[1,1];%系统1

num2=2;den2=[1,2];%系统2

[nums,dens]=series(num1,den1,num2,den2)%实现两个系统级联

[nump,denp]=parallel(num1,den1,num2,den2)%实现两个系统并联

[numf,denf]=feedback(num1,den1,num2,den2)%实现两个系统反馈程序运行结果为：nums=002；dens=132nump=034；denp=132numf=012；denf=134因此，各系统的传递函数分别为：33可编辑课件PPT

在实际应用中，也可以把一个复杂的线性时不变（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源程序：34可编辑课件PPT级联型结构系数：sos=1.00000.536001.0000001.00000.00570.62191.000000g=2格型结构系数（反射系数）：K=0.25000.50000.333335可编辑课件PPTb=[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源程序：36可编辑课件PPT级联型结构系数：sos=1.0000-3.00002.00001.0000-0.2500-0.12501.00000.00009.00001.00001.00000.5000g=0.0625并联型结构系数：R=-5.0250-1.0750i-5.0250+1.0750i0.925027.1875P=-0.5000+0.5000i-0.5000-0.5000i0.5000-0.2500

K=-1837可编辑课件PPT38可编辑课件PPT4.5线性时不变系统的响应4.5.1线性时不变系统的时域响应1．连续LTI系统的响应2．离散LTI系统的响应用MATLAB中的卷积函数conv()来实现。用MATLAB中的卷积函数conv()来实现。39可编辑课件PPTdt=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;40可编辑课件PPTx=ones(1,10);n=[0:14];h=0.5.^n;y=conv(x,h);stem(y);xlabel('n');ylabel('y[n]');41可编辑课件PPTnx=-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('x[n]');xlabel('n');xlabel('x[n]');subplot(222);stem(nh,h);title('h[n]');xlabel('n');xlabel('h[n]');subplot(223);stem(ny,y);title('x[n]*h[n]');xlabel('n');xlabel('y[n]');h=get(gca,'position');h(3)=2.5*h(3)set(gca,'position',h)改写为：42可编辑课件PPT格式：[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()

43可编辑课件PPTnum=[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('方波响应');

44可编辑课件PPT45可编辑课件PPTnum=[2,-3.4,5.5];den=[1,-1.2,0.8]；u=randn(1,100);dlsim(num,den,u);title('随机噪声响应')

46可编辑课件PPT4.5.2LTI系统的单位冲激响应1.求连续LTI系统的单位冲激响应函数impulse()格式：[Y,T]=impulse(sys)或impulse(sys)功能：返回系统的响应Y和时间向量T，自动选择仿真的时间范围。其中sys可为系统传递函数、零极增益模型或状态空间模型。2.求离散系统的单位冲激响应函数dimpulse()格式：[y,x]=dimpulse(num,den)功能：返回多项式传递函数的单位冲激响应y向量和时间状态历史记录x向量。47可编辑课件PPTa=[-0.55,-0.78;0.78,0];b=[1;0];c=[1.96,6.45];d=[0];impulse(a,b,c,d);title('LTI系统的冲激响应')48可编辑课件PPTnum=[2,-3.5,1.5];den=[1,-1.7,0.3];dimpulse(num,den);title('离散LTI系统的冲激响应')49可编辑课件PPT4.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)表示的系统单位阶跃响应。50可编辑课件PPT4.6线性时不变系统的频率响应51可编辑课件PPT52可编辑课件PPT4.6线性时不变系统的频率响应1．求模拟滤波器Ha(s)的频率响应函数freqs()格式一：H＝freqs(B,A,W)

功能：计算由向量W(rad/s)指定的频率点上模拟滤器系统函数Ha(s)的频率响应Ha(jΩ)，结果存于H向量中。格式二：[H,w]＝freqs(B,A,M)功能：计算出M个频率点上的频率响应存于向量H中,M个频率存放在向量w中。freqs自动将这M个频率点设置在适当的频率范围。默认w和M时,freqs自动选取200个频率点计算。53可编辑课件PPT求该模拟滤波器的频率响应。MATLAB源程序如下：B=1;A=[12.61313.41422.61311];W=0:0.1:2*pi*5;freqs(B,A,W)[例4-31]已知某模拟滤波器的系统函数54可编辑课件PPT图4.30模拟滤波器的频率响应55可编辑课件PPT

[例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=[10000000–1];A=1;freqz(B,A)56可编辑课件PPT图4.31滤波器幅度和相位曲线

该程序运行所绘出的幅频与相频性曲线如图4.31所示。57可编辑课件PPT例：已知一系统传递函数为：求系统的单位冲激响应h(n),以及h(n)的幅频和相频响应。N=128;x=[1,zeros(1,N-1)];%产生单位冲激函数num=[0.01-0.030.06-0.030.01];den=[12.531.80.5];y=dlsim(num,den,x);%计算单位冲激响应figure(1);n=1:N;stem(n,y);gridon;title('单位冲激响应');figure(2)freqz(num,den,N);%绘制幅频和相频响应曲线gridon;dimpulse(num,den)单位冲激响应函数，不带输出变量时，直接绘出单位冲激响应58可编辑课件PPT单位冲激响应幅频和相频响应冲激响应59可编辑课件PPT3．滤波函数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)60可编辑课件PPT图4.32系统对信号的响应

该程序运行所得结果如图4.32所示61可编辑课件PPT时域信号特性频域特性变换名称非周期连续信号非周期连续频谱傅里叶变换周期性连续信号非周期离散频谱傅里叶级数非周期离散信号周期连续频谱序列傅里叶变换周期离散信号周期离散频谱离散傅里叶级数离散信号（有限样本点）周期离散频谱离散傅里叶变换傅里叶变换：是建立以时间为自变量的“信号”与以频率为自变量的“频谱函数”之间的某种对应关系。遵循以下原则：时域连续周期离散非周期非周期离散周期连续频域4.7傅里叶(Fourier)变换62可编辑课件PPT63可编辑课件PPT64可编辑课件PPT65可编辑课件PPT4.7傅里叶(Fourier)变换4.7.1连续时间、连续频率－傅里叶变换正变换：

逆变换：

这说明求和的问题可以用x(t)行向量乘以列向量来实现。式中是t的增量，在程序中用dt表示.例：求如图矩形脉冲信号f(t)在-40~40rad/s区间的频谱。66可编辑课件PPTclc;clearall;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)gridon;axis([0,10,0,1.1]);subplot(122)plot(w,abs(F),'linewidth',1.5);gridon;67可编辑课件PPT取时间分隔的点数N=64需求的频谱宽度wf=40需求的频谱点数Nf=256取时间分隔的点数N=256需求的频谱宽度wf=40需求的频谱点数Nf=64采样较密采样稀，有频率泄漏68可编辑课件PPT4.7.2连续时间、离散频率－傅里叶级数正变换：

逆变换：

式中为离散频率相邻两谱线间的角频率间隔，k为谐波序号。4.7.3离散时间、离散频率－离散傅里叶级数正变换：

逆变换：

69可编辑课件PPT4.7.5离散时间、离散频率－离散傅里叶变换(DFT)正变换：

逆变换：

4.7.4时间离散、连续频率－序列傅里叶变换正变换：

逆变换：

注意，对于序列傅里叶变换，如果x(n)为无限长，那么就不能用MATLAB直接利用上式计算，只可以用它对表达式在频率点上求值，再画出它的幅度和相位。如果x(n)为有限长，就可直接用MATLAB计算70可编辑课件PPT4.7.5离散时间、离散频率－离散傅里叶变换(DFT)在MATLAB中，旋转因子矩阵可以由dftmtx函数实现。71可编辑课件PPTdftmtx函数:调用格式：W=dftmtx(n)功能：函数调用后，返回一个n*n的离散傅里叶变换DFT的旋转因子矩阵W。这样对于给定长度为n的行向量信号x(n),利用Xk=x*W就可以获得x(n)的傅里叶变换X(k)。其实，旋转因子矩阵也可以自己通过表达式构造。例如：编写函数，首先利用MATLAB计算序列x(n)的DFT和IDFT函数.72可编辑课件PPTfunctionxn=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:functionxk=dft(xn)N=length(xn);WNnk=dftmtx(N);xk=xn*WNnk;DFT:WNnk=conj(dftmtx(N))73可编辑课件PPTDFT参数对频率分辨率的影响：频率分辨率就是，在频率轴上能分辨出的两个频率点的最小间隔要能区分两个频率点f1,f2,，有效数据长度N必须满足：可以对序列进行补零操作，以增加数据的长度，这样做虽然不能提高频率分辨率，但补零后对原来的X(k)起到了插值作用，可以平滑频谱的包络。但是，如果增加有效数据的长度，则可以提高频率分辨率。74可编辑课件PPT利用上面编写的DFT函数，求的频谱特性，其中取fs=20Hz75可编辑课件PPTN=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);gridon;title('x(n)(0=<n<128)')subplot(212);plot(k,abs(xk(1:N/2)));gridon;76可编辑课件PPTN=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);gridon;title('x(n)补零后');subplot(212);plot(m,abs(xk1(1:N1/2)));gridon;title(‘x(n)补零后’)77可编辑课件PPTN=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);gridon;title('x(n)(0=<n<401)')subplot(212);plot(k,abs(xk(1:N/2)));gridon;

78可编辑课件PPTN=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));gridon;subplot(212)k=(0:N-1)*fs/N;plot(k,abs(xk));gridon;79可编辑课件PPT1．一维快速正傅里叶变换函数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的左右两半进行交换，从而将零频分量移至频谱中心。,80可编辑课件PPTfs=1000;

%采样频率t=0:1/fs:0.5;%时间轴取值范围及间隔x=cos(160*pi*t);N=512;

%FFT点数subplot(311)；plot(t,x);gridon;y=fft(x,N);

%N点FFT[x(t)]xk=fftshift(abs(y));

%求幅度谱并将零频分量移至频谱中心位置f=(-N/2:N/2-1)*fs/512;

%频率范围subplot(312)plot(f,xk);%做频谱图gridon;%用IFFT恢复原始信号z=real(ifft(y,N));ti=[0:length(z)-1]/fs;subplot(313)plot(ti,z);gridon;例：求的FFT序列和IFFT序列81可编辑课件PPT82可编辑课件PPT[例4-36]用快速傅里叶变换FFT计算下面两个序列的卷积。图4.35快速卷积框图83可编辑课件PPTMATLAB程序(部分)：%线性卷积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));

ticXk=fft(xn,L);

Hk=fft(hn,L);

Yk=Xk.*Hk;

yn=ifft(Yk,L);

toc

图4.36x(n),h(n)及其线性卷积波形Elapsedtimeis0.005074seconds.Elapsedtimeis0.000072seconds.84可编辑课件PPT练习求x(t)在有噪声和无噪声条件下的FFT.(提示：噪声由randn函数生成)1.85可编辑课件PPT4.8IIR数字滤波器的设计方法1.数字滤波器的频率响应函数幅度响应：相位响应：图4.37理想低通、高通、带通、带阻数字滤波器幅度特性

86可编辑课件PPT2.滤波器的技术指标幅度响应指标、相位响应指标

图4.38数字低通滤波器的幅度特性

通带要求：

阻带要求：

通带最大衰减：阻带最小衰减：87可编辑课件PPT4.8.1冲激响应不变法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.00000.09420.21580.0311den1=1.0000-2.00321.9982-0.76120.106988可编辑课件PPTMATLAB信号处理工具箱中的专用双线性变换函数bilinear()格式：[numd,dend]＝bilinear(num,den,Fs)功能：把模拟滤波器的传递函数模型转换成数字滤波器的传递函数模型。4.8.2双线性变换法双线性变换利用频率变换关系：

[例4-38]MATLAB源程序如下：

num=1;%模拟滤波器系统函数的分子

den=[1,sqrt(3),sqrt(2),1];%模拟滤波器系统函数的分母

[num1,den1]=bilinear(num,den,1)%求数字滤波器的传递函数运算的结果如下：num1=0.05330.15990.15990.0533den1=1.0000-1.33820.9193-0.154689可编辑课件PPT4.8.3IIR数字滤波器的频率变换设计法1.IIR数字滤波器的频率变换设计法的基本原理根据滤波器设计要求，设计模拟原型低通滤波器，然后进行频率变换，将其转换为相应的模拟滤波器(高通、带通等)，最后利用冲激响应不变法或双线性变换法，将模拟滤波器数字化成相应的数字滤波器。

图4.39IIR数字滤波器MATLAB设计步骤流程图

90可编辑课件PPT1．MATLAB的典型设计利用在MATLAB设计IIR数字滤波器可分以下几步来实现(1)按一定规则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标；(2)根据转换后的技术指标使用滤波器阶数函数，确定滤波器的最小阶数N和截止频率Wc；(3)利用最小阶数N产生模拟低通滤波原型；(4)利用截止频率Wc把模拟低通滤波器原型转换成模拟低通、高通、带通或带阻滤波器；(5)利用冲激响应不变法或双线性不变法把模拟滤波器转换成数字滤波器。[例4-39]设计一个数字信号处理系统，它的采样率为Fs＝100Hz，希望在该系统中设计一个Butterworth型高通数字滤波器，使其通带中允许的最小衰减为0.5dB，阻带内的最小衰减为40dB，通带上限临界频率为30Hz，阻带下限临界频率为40Hz。

91可编辑课件PPTMATLAB源程序设计如下：%把数字滤波器的频率特征转换成模拟滤波器的频率特征

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所示。

92可编辑课件PPT2．MATLAB的直接设计图4.39IIR数字滤波器MATLAB设计步骤流程图

[例4-41]试设计一个带阻IIR数字滤波器，其具体的要求是：通带的截止频率：wp1＝650Hz、wp2＝850Hz；阻带的截止频率：ws1＝700Hz、ws2＝800Hz；通带内的最大衰减为rp＝0.1dB；阻带内的最小衰减为rs＝50dB；采样频率为Fs＝2000Hz。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');%求滤波器阶数