2022年数字信号处理基本知识点Matlab实现.pdf

《2022年数字信号处理基本知识点Matlab实现.pdf》由会员分享，可在线阅读，更多相关《2022年数字信号处理基本知识点Matlab实现.pdf（29页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、数字信号处理（第二版）绪论1.4 MATLAB 在信号处理中的应用简介MATLAB是美国 Mathworks 公司于 1984 年推出的一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算、信号处理、系统仿真和图形显示于一体，从而被广泛地应用于科学计算、控制系统、信息处理等领域的分析、仿真和设计工作。MATLAB软件包括五大通用功能：数值计算功能（Numeric ），符号运算功能（Symbolic ）；数据可视化功能（Graphic），数据图形文字统一处理功能（Notebook）和建模仿真可视化功能（ Simulink ）。该软件有三大特点：一是功能强大；二是界面友善、语言自然；三是开放性

2、强。目前， Mathworks 公司已推出30 多个应用工具箱。MATLAB 在线性代数、矩阵分析、数值及优化、数理统计和随机信号分析、电路与系统、系统动力学、信号和图像处理、控制理论分析和系统设计、过程控制、建模和仿真、通信系统、 以及财政金融等众多领域的理论研究和工程设计中得到了广泛应用。2.10 离散时间信号与系统的Matlab 表示2.10.1 离散时间信号的表示和运算1、基本序列的Matlab 表示单位采样序列在 MATLAB 中，单位采样序列可以通过编写以下的DTimpulse.m 文件来实现，即function y=DTimpulse (n) y=(n=0); %当参数为 0 时

3、冲激为1,否则为 0 调用该函数时n 必须为整数或整数向量。单位阶跃序列在 MATLAB 中，单位阶跃序列可以通过编写DTu.m 文件来实现，即function y=DTu (n) y=n=0; %当参数为非负时输出1调用该函数时n 必须为整数或整数向量。矩形序列用 MATLAB 表示矩形序列可根据公式( )( )()NRnu nu nN并利用 DTu 函数生成，即function y=DTR(n,N) y=DTu(n)-DTu(n-N); 调用该函数时n 必须为整数或整数向量，N 必须为整数。实指数序列用 MATLAB 表示实指数序列( )( ),nx na u nnN aR，即精品资料 -

4、 - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 1 页，共 29 页 - - - - - - - - - - n=0:N-1; x=a.n; 复指数序列用 MATLAB 表示实指数序列0()0( ),jnx nenNR，即n=0:N-1; x=exp(sigma+j*w0)*n); 正弦序列用 MATLAB 表示实指数序列00( )sin(), ,x nAnnN AR，即n=0:N-1; x=A*sin(w0*n+phi); 2、基本序列运算的Matlab 表示序列相加用 MATLAB 实现序列12( )( )( )

5、y nx nxn，即y=x1+x2; 当 x1 和 x2 序列长度不等或其起始位置不对应时，序列相加前需要调整x1 和 x2，使之具有相同长度并具有相同起始位置。序列相乘用 MATLAB 实现序列12( )( )*( )y nx nx n，即y=x1.*x2; 当 x1 和 x2 序列长度不等或其起始位置不对应时，序列相乘前需要调整x1 和 x2，使之具有相同长度并具有相同起始位置。此处需要利用Matlab 中的点乘（ .*）实现。序列反转用 MATLAB 实现序列( )()y nxn，即y=flipr(x); 序列移位用 MATLAB 实现序列0( )()y nx nn，即n=-N:N; n

6、1=n-n0; subplot(211);stem(n,x, fill ); subplot(211);stem(n1,x, fill ); 2.10.2 离散时间信号的频域变换1、Matlab 实现离散时间傅里叶变换（DTFT ）10()( )NjjnnX ex n e序列的离散时间傅里叶变换（DTFT ）由于频率为连续实数变量，所以严格的讲其在精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 2 页，共 29 页 - - - - - - - - - - Matlab 下是无法实现的，但是为了说明DTFT

7、 在 Matlab 下的运算过程，这里将其频率变量设置为频率向量0121M，只计算频率向量指定的有限个频率点上的频谱结果。设序列(0)(1)(1)xxxxN，0121nN，0121M，则 DTFT 的矩阵表达形式为：011011011011000111(1)(1)(1)()()()(0)(1)(1)MMMMjjjjjjjjjjNjNjNX eX eX eeeeeeexxx Neee其中：011011012101100001111()2(1)(1)(1)1MMMMjjjjjjjjNjNjNNDTFT 在 Matlab 中可以利用矩阵乘法的形式表示如下：exp ()TXxjn【例 2-1】 用 M

8、atlab 实现5( )( )x nR n的离散时间傅里叶变换。解： DTFT 的 MA TLAB 源程序为：N=5; n=0:N-1; k=-200:200; w=(4*2*pi/400)*k; x= DTR(n,N) DTFT_X=x*exp(-j*(n*w); magX=abs(DTFT_X); plot(w/pi,magX); xlabel(w/pi); ylabel(|X(jw)|); title(x(n) 的离散时间傅里叶变换的幅频响应曲线); 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第

9、3 页，共 29 页 - - - - - - - - - - 图 2-1 离散时间傅里叶变换幅度响应曲线2、Matlab 实现 Z 变换Z 正反变换的符号表示MATLAB符号数学工具箱提供了计算离散时间信号单边Z 变换的函数ztrans 和 Z 反变换函数 iztrans，其语句格式分别为z=ztrans(x) x=iztrans(z) 上式中的 x 和 z 分别为时域表达式和z 域表达式的符号表示，可通过sym 函数来定义。【例 2-2】 试用 ztrans 函数求下列函数的z 变换。（1）)()cos()(nunanxn；（2）)()2(2)(11nunxnn。解： （1）Z 变换 MAT

10、LAB源程序为x=sym(an*cos(pi*n); Z=ztrans(x); simplify(Z) ans= z/(z+a) （2）Z 变换 MATLAB源程序为x=sym(2(n-1)-(-2)(n-1); Z=ztrans(x); simplify(Z) ans= z2/(z-2)/(z+2) 【例 2-3】 试用 iztrans 函数求下列函数的Z 反变换。（1）65198)(2zzzzX（2）32)2)(1()12112()(zzzzzzX解： （1）Z 反变换 MA TLAB 源程序为Z=sym(8*z-19)/(z2-5*z+6); x=iztrans(Z); 精品资料 - -

11、 - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 4 页，共 29 页 - - - - - - - - - - simplify(x) ans= -19/6*charfcn0(n)+5*3(n-1)+3*2(n-1)其中， charfcn0(n) 是)(n函数在 MA TLAB 符号工具箱中的表示，反变换后的函数形式为)()2335()(619)(11nunnxnn。（2）Z 反变换 MATLAB 源程序为Z=sym(z*(2*z2-11*z+12)/(z-1)/(z-2)3); x=iztrans(Z); simplif

12、y(x) ans= -3+3*2n-1/4*2n*n-1/4*2n*n2 其函数形式为)()241241233()(2nunnnxnnn。Z 反变换的留数表示如果信号的z 域表示式)(zX是有理函数，进行z 反变换的另一个方法是对)(zX进行部分分式展开，然后求各简单分式的z 反变换。设)(zX的有理分式表示为)()(1)(221122110zAzBzazazazbzbzbbzXnnmm（4-3）MATLAB信号处理工具箱提供了一个对)(zX进行部分分式展开的函数residuez，其语句格式为R,P,K=residuez(B,A) 其中， B，A 分别表示)(zX的分子与分母多项式的系数向量；

13、R 为部分分式的系数向量；P为极点向量； K 为多项式的系数。若)(zX为有理真分式，则K 为零。【例 2-4】 试用 MATLAB 命令对函数321431818)(zzzzX进行部分分式展开，并求出其Z 反变换。解： MATLAB源程序为B=18; A=18,3,-4,-1; R,P,K=residuez(B,A) R= 0.3600 0.2400 0.4000 P= 0.5000 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 5 页，共 29 页 - - - - - - - - - - -0.333

14、3 -0.3333 K= 从运行结果可知，32pp，表示系统有一个二重极点。所以，)(zX的部分分式展开为2111)3330.314.03333.0124.05 .0136.0)(zzzzX（因此，其 Z 反变换为)()3333.0)(1(4. 0)3333.0(24.0)5.0(36.0)(nunnxnnnZ 变换的零极点分析离散时间系统的系统函数定义为系统零状态响应的Z 变换与激励的Z 变换之比，即)()()(zXzYzH（4-4）如果系统函数)(zH的有理函数表示式为11211121)(nnnnmmmmazazazabzbzbzbzH（4-5）那么，在 MATLAB中系统函数的零极点就可

15、通过函数roots 得到，也可借助函数tf2zp 得到，tf2zp 的语句格式为Z,P,K=tf2zp(B,A) 其中， B 与 A 分别表示)(zH的分子与分母多项式的系数向量。它的作用是将)(zH的有理分式表示式转换为零极点增益形式，即)()()()()(2121nmpzpzpzzzzzzzkzH（4-6）若要获得系统函数)(zH的零极点分布图，可直接应用zplane 函数，其语句格式为zplane(B,A) 其中，B 与 A 分别表示)(zH的分子和分母多项式的系数向量。它的作用是在Z 平面上画出单位圆、零点与极点。【例 2-5】 已知某离散因果LTI 系统的系统函数为16.032.0)

16、(2zzzzH试用 MATLAB命令求该系统的零极点并绘制其零极点分布图。解：用 tf2zp 函数求系统的零极点，MATLAB源程序为精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 6 页，共 29 页 - - - - - - - - - - B=1,0.32; A=1,1,0.16; R,P,K=tf2zp(B,A) zplane(B,A),grid on; legend(零点 ,极点 ); title( 零极点分布图 ); 零极点计算结果如下所示。R= -0.3200 P= -0.8000 -0.20

17、00 K= 1 零极点分布如图2-2 所示。可见，该因果系统的极点全部在单位圆内，故系统是稳定的。图 2-2 系统零极点分布图2.10.3 离散时间系统的表示1、离散时间系统的时域表示方法离散时间LTI 系统在时域中常常可通过线性常系数差分方程和单位脉冲响应两种方式描述。线性常系数差分方程MjjNiijnxbinya00)()(其中，ia（0i，1，, ， N）和jb（0j，1， , ， M）为实常数。MATLAB中函数 filter 可对差分方程在指定时间范围内的输入序列所产生的响应进行求解。函数filter 的语句格式为y=filter(b,a,x) 精品资料 - - - 欢迎下载 - -

18、 - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 7 页，共 29 页 - - - - - - - - - - 其中， x 为输入的离散序列；y 为输出的离散序列；y 的长度与 x 的长度一样； b 与 a 分别为差分方程右端与左端的系数向量。【例 2-6】 已知某 LTI 系统的差分方程为)1(2)()2(2)1(4)(3nxnxnynyny用 MATLAB 命令绘出当输入激励信号为)()2/1()(nunxn时，该系统的零状态响应。解： MATLAB源程序为a=3 -4 2; b=1 2; n=0:30; x=(1/2).n; y=filt

19、er(b,a,x); stem(n,y,fill),grid on xlabel(n),title( 系统响应 y(n) 程序运行结果如图2-3 所示。图 2-3 系统零状态响应单位脉冲响应系统的单位脉冲响应定义为系统在)(n激励下系统的零状态响应，用)(nh表示。MATLAB求解单位脉冲响应可利用函数filter ， 并将激励设为前面所定义的DTimpulse 函数。例如用 filter 函数求解例2-6 中系统的单位脉冲响应的MA TLAB 程序为a=3 -4 2; b=1 2; n=0:30; x=impDT(n); h=filter(b,a,x); stem(n,h,fill),gri

20、d on xlabel(n),title( 利用 filter 函数求解的系统单位脉冲响应h(n) 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 8 页，共 29 页 - - - - - - - - - - MATLAB另一种求单位脉冲响应的方法是利用控制系统工具箱提供的函数impz 来实现。 impz 函数的常用语句格式为impz(b,a,N) 其中，参数 N 通常为正整数，代表计算单位脉冲响应的样值个数。例如用impz 函数求解例2-6 中系统的单位脉冲响应的MATLAB 程序为a=3 -4 2;

21、b=1 2; n=0:30; impz(b,a,30),grid on xlabel(n),title( 利用 impz 函数求解的系统单位脉冲响应h(n) 单位脉冲响应的因果稳定性判断单位脉冲响应)(nh作为时域中描述系统的一种方法，可以用它来分析系统的因果性和稳定性。系统的因果性可以直接根据)(nh满足因果关系的成立条件直接进行判定；系统的稳定性一般不太容易直接检查系统的单位脉冲响应是否满足绝对可和，可行的方法是在系统的输入端加入单位阶跃序列，如果系统的输出趋近一个常数（包括零），就可以判定系统的稳定性。这种判定方式的思想是，如果系统稳定，信号加入系统后，系统输出的开始一段成为暂态响应，随

22、着n 的加大，幅度趋于稳定，达到稳态输出。【例 2-7】设一谐振器的差分方程为00( )1.8237 (1)0.9801 (2)( )(2)y ny ny nb x nb x n其中01100.49b，谐振器的谐振频率为0.4rad。试用 Matlab 实验方法检查系统是否稳定？解： Matlab 源程序为N=400; b0=1/100.49; n=0:N-1; B=b0,0,-b0; A=1,-1.8237,0.9801; x=DTu(n); y=filter(B,A,x); subplot(211); stem(n,x,fill); title( 输入信号 -单位阶跃序列u(n); axi

23、s(0,N,0,1.2); subplot(212); stem(n,y,fill); title( 输出信号 -u(n) 系统响应 ); axis(0,N,-0.05,0.05); 程序执行结果如图2-4 所示，可见，当此系统输入单位阶跃序列后，系统的响应呈现出幅值递减状态，可以预见当n时，幅值趋近于一个常数（零常数），故而可以判定此系统稳定。精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 9 页，共 29 页 - - - - - - - - - - 图 2-4 输入单位阶跃序列的系统响应2、离散时间系

24、统的频域表示方法离散时间系统的频率响应定义为)(|)(|)(jjjeeHeH（ -7）其中，|)(|jeH称为离散时间系统的幅频响应，在Matlab 中利用函数abs对频率响应进行求解；)(称为离散时间系统的相频响应，在Matlab 中利用函数angle 对频率响应进行求解。由于)(jeH是以2为周期的周期函数，故只需分析)(jeH在一个周期范围内的特征便可分析出系统的整个频率特性。MATLAB为求解离散时间系统频响特性的提供了函数freqz，调用freqz 的格式主要有两种。第一种方式为：H,w=freqz(B,A,N) 其中， B 与 A 分别表示)(zH的分子和分母多项式的系数向量；N

25、为正整数， 默认值为512；返回值w 包含,0范围内的N 个频率等分点；返回值H 则是离散时间系统频率响应)(jeH在0范围内 N 个频率处的值。第二种方式为：H,w=freqz(B,A,N, whole ) 与第一种方式不同之处在于角频率的范围由,0扩展到2,0。【例 2-7】 用 MATLAB命令绘制系统8109.056.19028.096.0)(22zzzzzH的频率响应曲线。精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 10 页，共 29 页 - - - - - - - - - - 解： MAT

26、LAB源程序为b=1 -0.96 0.9028; a=1 -1.56 0.8109; H,w=freqz(b,a,400,whole); Hm=abs(H); Hp=angle(H); subplot(211) plot(w,Hm),grid on xlabel(omega(rad/s),ylabel(Magnitude) title( 幅频响应曲线 ) subplot(212) plot(w,Hp),grid on xlabel(omega(rad/s),ylabel(Phase) title( 相频响应曲线 ) 程序运行结果如图2-5 所示。图 2-5 系统的频率响应2.10.4 离散时间

27、系统的输入输出关系1、 用 filter 函数实现系统输出利用 2.10.3 节中介绍的filter 函数可以对输入数据进行滤波处理，此函数适用于IIR 和 FIR滤波器。滤波器的z 域表示为11(1)(2)(1)( )( )(1)(2)(1)NMbbzb NzY zX zaaza Mz其中(1)1a，如果(1)1a，则 Matlab 将自动进行归一化处理；如果(1)0a，则 Matlab会给出出错提示。filter 函数的调用格式有：y=filter(b,a,x) ，此方式利用给定的描述系统的系数向量a 和 b，对输入信号x 中的数据精品资料 - - - 欢迎下载 - - - - - - -

28、 - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 11 页，共 29 页 - - - - - - - - - - 进行滤波，将结果存入y 输出向量中。y,zf=filter(b,a,x) ，此方式调用方式与前一个方式基本一致，不同之处在于除了返回滤波处理结果y 向量外，还得到输入信号x 的最终状态向量zf。y=filter(b,a,x,zi) ，此方式调用方式还可以增加对输入信号x 的初始状态的表示，用 zi 描述。2、 用 conv 函数实现系统输出离散时间系统的输出可以利用输入信号和描述系统的单位脉冲响应的线性卷积表示，即mmnhmxnhnxny)()()(*

29、)()(（3-2）MATLAB求离散时间信号线性卷积的命令为conv，其语句格式为y=conv(x,h) 其中，x 与 h表示离散时间信号的输入向量和单位脉冲响应向量；y 为输出向量。 用 MATLAB进行线性卷积运算时，只能计算信号有限项的卷积，无法实现无限项的卷积。【例 2-8】 已知 LTI 系统的单位脉冲响应为88 .0nununhn，用 MATLAB求当输入信号为)4()()(nununx时，系统的零状态响应。解：依题可知，)(nh向量的长度至少为8，)(nx向量的长度至少为4，为了使图形显示美观完整，将)(nh向量和)(nx向量加上一些附加的零值。MATLAB源程序为nx=-1:5

30、; %x(n) 向量显示范围 (添加了附加的零值) nh=-2:10; %h(n)向量显示范围 (添加了附加的零值) x=uDT(nx)-uDT(nx-4); h=0.8.nh.*(uDT(nh)-uDT(nh-8); y=conv(x,h); ny1=nx(1)+nh(1); %卷积结果起始点%卷积结果长度为两序列长度之和减1,即 0 到(length(nx)+length(nh)-2) %因此卷积结果的时间范围是将上述长度加上起始点的偏移值ny=ny1+(0:(length(nx)+length(nh)-2); subplot(311) stem(nx,x,fill),grid on xl

31、abel(n),title(x(n) axis(-4 16 0 3) subplot(312) stem(nh,h,fill),grid on xlabel(n),title(h(n) axis(-4 16 0 3) subplot(313) stem(ny,y,fill),grid on xlabel(n),title(y(n)=x(n)*h(n) 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 12 页，共 29 页 - - - - - - - - - - axis(-4 16 0 3)2.10.5

32、离散时间系统的网络结构数字信号处理中一般涉及到的网络结构可分为无限长单位脉冲响应（IIR ）网络结构和有限长单位脉冲响应（FIR）网络结构。 FIR 网络中一般不存在输出对输入的反馈支路，即其n时刻的输出只与n时刻及其附近的输入有关，而与n时刻附近的输入无关，而IIR 网络与 FIR 正好相对， 其网络中一般存在反馈支路，并且n时刻的输出不仅与n时刻及其附近的输入有关， 还与n时刻附近的输出也相关。不同的系统网络结构实际上代表不同的算法，而算法直接影响了系统的运算误差、运算速度和运算复杂度，因此系统实现过程需要根据具体的应用和技术指标确定选择哪种具体的网络结构。1、IIR 网络结构IIR 系统

33、的基本网络结构通常有：直接型、级联型和并联型。直接型网络在 Matlab 中，直接型结构有分别描述系统函数的分子、分母系数的两个向量B 和 A 表示。如果数字滤波器的直接型系统函数关系为：图 2-6 利用 conv 函数求解系统的零状态响应精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 13 页，共 29 页 - - - - - - - - - - 01( )MiiiNiiib zH za z则向量012,NAa a aa，向量012,MBb b ab。一般默认01a，如果01a，需要将0a转换为 1。M

34、atlab 中，实现IIR 滤波器的直接型可通过调用filter 函数来实现。如果滤波器输入信号为向量IX，输出信号向量为OY，则( ,)OIYfilter B A X按照直接型结构实现对IX的滤波，计算输入信号IX的零状态响应输出信号OY。级联型网络结构在 Matlab 中，级联型的网络结构可以通过调用tf2sos 函数实现直接型到级联型的转换，同样也可以调用sos2tf 函数实现级联型到直接型的转换。除了tf2sos/sos2tf 函数外，还可以利用 tf2zp/zp2tf 函数，两者的区别主要是生成的级联型结构，tf2sos 函数转换实现了二阶级联型结构，而tf2zp 函数转换实现了零极

35、点级联型结构。下面以tf2sos 函数为例说明。tf2sos 函数调用方式为,2(,)S Gtfsos B A，其中向量B，A 分别表示直接型系统函数的分子、分母多项式系数。当1A时，表示 FIR 滤波器系统函数，否则表示IIR 滤波器系统函数。 tf2sos 函数调用返回L 级二阶级联型结构的系数矩阵S和增益常数G。0111211121010101122201212111LLLLLbbbaabbbaaSbbbaaS 为 L6 的矩阵，其中每一行表示一个二阶级联子系统函数的分子、分母系数向量，例如第 k 行所对应的二阶子系统函数为：120121212( )1, 1kkkkkkbb zb zHz

36、kLa za ztf2sos 函数最终输出的级联型结构系统函数为：12( )( )( )( )LH zGHz HzHz其中 G 表示不为零的增益常数，也可理解为一个系统函数为常数的子系统。并联型网络结构在Matlab中，直接型到并联型的转换目前官方还没有提供转换函数，而数字信号处理教程 MATLAB 释义与实现一书中提供的tf2par函数实现了直接型到并联型转换，并且还提供了 par2tf完成并联型到直接型的转换。2、FIR 网络结构FIR 系统的基本网络结构通常有：直接型、级联型、频率采样型和线性相位型。其中直精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下

37、载 名师归纳 - - - - - - - - - -第 14 页，共 29 页 - - - - - - - - - - 接型、级联型与IIR 系统的直接型、级联型的实现函数相同，分别是filter 函数和 tf2sos 函数，它们的区别主要是系统函数分母多项式系数向量A 是否等于1。 如果1A， 则表示 FIR系统，否则表示IIR 系统。FIR 系统中还有一种易于模块化设计实现的网络结构格型网络结构，此结构适用于对有限字长效应不敏感的时域离散系统，并且适合递推运算，它在一般的数字滤波器、自适应滤波器、线性预测系统及谱估计中都有广泛的应用。在Matlab 中，直接型到格型网络结构的转换可以通过调

38、用tf2latc 函数实现。由于本书中未涉及格型网络结构，在此仅给出调用的函数名称，具体实现可参考其它资料。3.5 离散傅里叶变换的Matlab 实现1、离散傅里叶变换离散傅立叶变换不仅在理论上有着重要的意义，而且因为其具有快速计算的方法，所以在各种数字信号处理的实际应用中也具有重要的实用意义。有限长序列)(nx的离散傅立叶变换（DFT ）为10,)()()(10NnWnxnxDFTkXNnknN逆变换为10,)(1)()(10NnWkXNkXIDFTnxNnknN若将 DFT变换的定义写成矩阵形式，则得到x，其中 DFT变换矩阵为2)1(111.1.11.11NNNNNNNWWWWAMatl

39、ab 中可借助 Dftmtx 函数计算DFT变换矩阵 A的函数，调用方式如下（1）Adftmtx （n） ：返回 nn 的 DFT变换矩阵 A。若 x 为给定长度的行向量，则yx*A，返回 x 的 DFT变换 y。（2）Ai conj （dftmtx （n） ）/n ：返回 nn 的 IDFT 变换矩阵Ai。【例 3-1】如果)4/sin()8/sin()(nnnx是一个长度N16 的序列，用MATLAB 求其DFT的结果，并画出其结果图。解： Matlab程序如下N=16; n=0:1:N-1; %时域采样xn=sin(n*pi/8)+sin(n*pi/4); 精品资料 - - - 欢迎下载

40、 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 15 页，共 29 页 - - - - - - - - - - Yk=dftmtx(N)*xn; subplot(2,1,1) stem(n,xn,fill);title(序列 x(n); subplot(2,1,2) stem(k,abs(Yk),fill);title(序列 x(n) 的 DFT结果 ); 图 3-1 调用 Dftmtx 函数计算DFT 2、离散傅里叶变换的性质离散傅里叶变换的性质有很多，这里主要针对离散傅里叶变换与其它变换最主要的两个不同特性进行说明，它们分别是圆周

41、移位和圆周卷积。圆周移位一个有限长序列)(nx的圆周移位定义( )mNNxxnmRn式中，Nxnm表示将)(nx的周期延拓序列( )x n的移位 N个单位。【例 3-2】求有限长序列,)4 .0(8)(nnx200n的圆周移位)()10()(2020nRnxnxm，并画出其结果图。解： Matlab 程序如下N=20; m=10; n=0:1:N-1; x=8*(0.4).n; n1=mod(n+m),N); xm=x(n1+1); 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 16 页，共 29 页

42、- - - - - - - - - - subplot(2,1,1) stem(n,x); title(原始序列 ); xlabel(n); ylabel(x(n); subplot(2,1,2) stem(n,xm); title(圆周移位后序列); xlabel(n); ylabel(x(n+10)mod20); 图 3-2 圆周移位前后比较圆周卷积MATLAB 中构造用于计算圆周卷积的函数circonvt程序如下function y=circonvt(x1,x2,N) if(length(x1)N) error(N should bigger than or equal to the l

43、ength of x1!) end if(length(x2)N) error(N should bigger than or equal to the length of x2!) x1=x1,zeros(1,N-length(x1) x2=x2,zeros(1,N-length(x2) m=0:1:N-1 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 17 页，共 29 页 - - - - - - - - - - for n=1:1:N end y=x1*H4.7 FFT 的 Matlab 实现快速

44、傅立叶变换（FFT）并不是与DFT 不同的另外一种变换，而是为减少DFT 计算次数的一种快速有效的算法。快速傅立叶变换算法形式很多，但是基本上可以分为两大类，即按时间抽取（ Decimation-In-Time ，简称 DIT ）法和按频率抽取（Decimation-In-Frequency ，简称 DIF）法。MATLAB为计算 FFT 提供的函数主要有fft 、 ifft 、fft2 、ifft2, fftn 、ifftn 和 fftshift 、ifftshift等，其中 fft、ifft 主要针对1 维数据进行处理，而fft2 、ifft2, fftn 、ifftn 主要针对2 维或

45、n 维数据进行处理。当所处理的数据的长度为2 的指数幂时，采用基-2 算法进行计算，计算速度会显著提高。 所以，要尽可能使所要处理的数据长度为2 的指数幂或者用添零的方式来添补数据使之成为2 的指数幂。1、 FFT 变换和逆变换FFT 变换和逆变换在Matlab 中通过调用fft 和 ifft 函数实现，其中FFT 函数调用过程有三种情况：（1）Y=fft(X) 如果 X 是向量，则采用傅立叶变换来求解X 的离散傅立叶变换；如果 X 是矩阵，则计算该矩阵每一列的离散傅立叶变换；如果 X 是（ N)D维数组，则是对第一个非单元素的维进行离散傅立叶变换；（2）Yfft(X,N) N 是进行离散傅立

46、叶变换的X 的数据长度，可以通过对X 进行补零或截取来实现。（3）Yfft(X,dim) 或 Yfft(X,N,dim) 在参数 dim 指定的维上进行离散傅立叶变换；当 X 为矩阵时， dim 用来指定变换的实施方向：dim=1,表明变换按列进行；dim=2 表明变换按行进行。函数ifft 的参数应用与函数fft 完全相同。【例 4-1】fft 在信号分析中的应用，使用频率分析方法从受噪声污染的信号x(t)中鉴别出有用的信号。解： Matlab 程序如下：t=0:0.001:1; %采样周期为0.001s,即采样频率为1000Hz; 精品资料 - - - 欢迎下载 - - - - - - -

47、 - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 18 页，共 29 页 - - - - - - - - - - %产生受噪声污染的正县正弦波信号；x=sin(2*pi*100*t)+sin(2*pi*200*t)+rand(size(t); subplot(2,1,1) plot(x(1:50); %画出时域内的信号；Y=fft(x,512); %对 X 进行 512 点的傅立叶变换；f=1000*(0:256)/512; %设置频率轴（横轴）坐标，1000 为采样频率；subplot(2,1,2) plot(f,Y(1:257); %画出频域内的信号运行结果图

48、 4-1 时域信号和频域信号的比较由上面的图4-1 可以看出，从受噪声污染信号的时域形式中，很难看出正弦波的成分。但是通过对x(t) 作傅立叶变换， 把时域信号变换到频域进行分析，可以明显看出信号中100Hz和 200Hz 的两个主要频率分量。2、 FFT 频率搬移在 Matlab 中， FFT 变换后默认的频谱图频率成分并不直观，可以通过fftshift 和 ifftshift函数进行频率搬移，其中fftshift 函数的调用方式为：Z=fftshift(Y) 此函数可用于将傅立叶变换结果Y（频域数据）中的直流成分（即频率为0 处得值）移到频谱的中间位置。当fftshift 函数调用参数不同

49、时会有以下三种情况：如果 X 是向量，则变换Y 的左右两边；如果 X 是矩阵，则交换Y 的一、三象限和二、四象限；精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 19 页，共 29 页 - - - - - - - - - - 如果 Y 是多维数组，则在数组的每一维交换其“半空间”。 函数 ifftshift 的参数应用与函数 fftshift 完全相同。【例 4-2】设序列( )2sin(4)5cos(8)x nnn，求64N点的 FFT 的幅度谱，并画出其图。解： Matlab 程序如下：N=64; n

50、=0:N-1; w=2*pi*n/N; x=2*sin(4*pi*n)+5*cos(8*pi*n); y=fft(x,N); z=fftshift(y); subplot(211); plot(w/pi,abs(y) title(fft 频谱计算结果 ) axis(0,2,-10,350) subplot(212); plot(w/pi,abs(z) title(fftshift频谱搬移后结果) axis(0,2,-10,350) 图 4-2 FFT 频谱搬移前后比较精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - -