数字信号处理： MATLABdft对称性验证以及应用

数字信号处理： MATLABdft对称性验证以及应用 武汉理工大学《数字信号处理》课程说明书 目 录 1 MATLAB基本操作及常用命令介绍 .................................................. 1 1(1 MATLAB的启动 ..................................................................... 1 1(2桌面平台 .................................................................................. 1 1.3 基本平面图形绘制命令 plot ..................................................... 2 2 理论 ............................................................................................. 3 2.1实验内容 ..................................................................................... 3 2.2序列对称性的理论验证 .............................................................. 3 3 程序验证 ............................................................................................. 4 4 结果分析 ............................................................................................. 7 5 对称性的应用 ................................................................................... 10 5.1 FFT算法的基本思想 ................................................................ 10 5.2 对称性应用的程序实现 ........................................................... 11 6 心得体会 ........................................................................................... 15 参考文献 ............................................................................................... 16 武汉理工大学《数字信号处理》课程设计说明书 1 MATLAB基本操作及常用命令介绍 1(1 MATLAB的启动 启动MATLAB有多种方式，最常用的就是双击系统桌面的MATLAB图 标， 也可以在开始菜单的程序选项中选择MATLAB快捷方式。 初次启动MATLAB后，将进入MATLAB默认设置的桌面平台。 1(2桌面平台 默认设置情况下的桌面平台包括6个窗口，分别是MATLAB主窗口、命令窗口、历史窗口、当前目录窗口、发行说明书窗口和工作间管理窗口。下面分别对各窗口做简单介绍。 (1)MATLAB主窗口，与MATLAB的早期版本不同的是，MATLAB 6.1增加 了一个主窗口，其他的几个窗口都包含在这个大的主窗口中。 (2)命令窗口， 在主窗口的“View”菜单下选择“Command Window”命令，可以 打开或关闭MATLAB的命令窗口。其中“>>”为运算提示符，示MATLAB 正处在准备状态。当在提示符后输入一段运算式并按Enter键后，MATLAB 将给出计算结果，然后再次进入准备状态。 (3)工作间窗口 在主窗口的“View”菜单下选择“Workspace”命令，可以打开或关闭MATLAB 的工作间窗口，工作间窗口是MATLAB的重要组成部分。 MATLAB的程序类型有三种 1)脚本M文件: 在命令窗口中输入并执行，它所用的变量都要在工作空 间中获取，不需要输入输出参数的调用，退出MATLAB后就释放了。 2)程序M文件: 以.m格式进行存取，包含一连串的MATLAB指令和必 要的注解。需要在工作空间中创建并获取变量，也就是说处理的数据为 命令窗口中的数据，没有输入参数，也不会返回参数。程序运行时只需 1 武汉理工大学《数字信号处理》课程设计说明书 在命令窗口中键入其名称即可(不需要后缀名)。 3)函数M文件(常称为M函数):与在命令窗口中输入命令一样，函数接 受输入参数，然后执行并输出结果。用help命令可以显示它的注释说明。 4)M函数具有的基本结构: 1)函数定义行(关键字function) function[out1,out2,..]=filename(in1,in2,..) 输入和输出(返回)的参数个数分别由nargin和nargout两个MATLAB 保的变量来给出。 2)第一行帮助行以(%)开头，作为lookfor指令搜索的行 3)函数体说明及有关注解以(%)开头，用以说明函数的作用及有关内容 4)函数体语句:函数体内使用的除返回和输入变量这些在function语句中 直接引用的变量以外的所有变量都是局部变量，即在该函数返回之后， 这些变量会自动在MATLAB的工作空间中清除掉。如果希望这些中间 变量成为在整个程序中都起作用的变量，则可以将它们设置为全局变 量。 1.3 基本平面图形绘制命令 plot 1)plot(X,Y) %当X,Y均为实数向量，且为同维向量，则以x的元素为横 坐标，y的元素为纵坐标，绘制曲线。若X，Y均为同维同型实数矩阵，X = [X(i)]，Y = [Y(i)]，其中X(i),Y(i)为列向量，则plot(X,Y)依次画出plot(X(i),Y(i))， 矩阵有几列就有几条线; 2)plot(X1,Y1,X2,Y2,…)，其中Xi与Yi成对出现，plot(X1,Y1,X2,Y2,…)将 分别按顺序取两数据Xi与Yi进行画图。 3)plot(X1,Y1,LineSpec1,X2,Y2,LineSpec2…) 将按顺序分别画出由三参数定 义Xi,Yi,LineSpeci的线条。其中参数LineSpeci指明了线条的类型，标记符 2 武汉理工大学《数字信号处理》课程设计说明书 2 理论分析 2.1实验内容 1,编写程序验证DFT的正确性 2,利用DFT的对称性用一次FFT实现两个序列的FFT变换 2.2序列对称性的理论验证 (1) 有限长共轭对称序列和共轭反对称序列 定义: 共轭对称序列:x ( n ) = x *( N,n ) , 0 ?n? N ,1 e p e p 共轭反对称序列:x (n) =,x*( N,n ) ,0?n?N,1 o p o p 如下图所示: 图2-1 共轭对称与共轭反对称序列示意图 3 武汉理工大学《数字信号处理》课程设计说明书 (2) DFT的共轭对称性 1)设x(n)=x(n)+jx(n) ri 其中x(n)=Re[x(n)]=[x(n)+x*(n)]/2 , r x(n)=Im[x(n)]=[x(n),x*(n)]/2 i 则 DFT[x(n)]=Xep(k) r DFT[jx(n)]=Xop(n) i 2) 设x(n)=x(n)+jx(n) epop 其中x(n)=[x(n)+x*(N,n)]/2 , ep x(n)=[x(n),x*(N,n)]/2 i 则 DFT[x(n)] =Re[X(k)] ep DFT[x(n)]=Im[X(k)] op 结论:如果序列x(n)的DFT为X(k),则x(n)的实部和虚部(包括j)的DFT分 别为X(k)的共轭反对称分量和共轭对称分量;而x(n)的共轭反对称分量和共 轭对称分量的DFT分别为X(k)的实部和虚部乘以j 。 3 程序验证 %定义分解函数 function[xec,xoc]=xhfenjie(x)% % %该函数实现将信号分解为奇分量和偶分量的功能 %x为输入信号序列 %xec为奇分量序列 %xoc为偶分量序列 % 4 武汉理工大学《数字信号处理》课程设计说明书 %首先进行判别 if any(imag(x)~=0) error('x不是实序列'); end %进行信号分解 N=length(x); n=0:N-1; xec=0.5*(x+x(mod(-n,N)+1)); xoc=0.5*(x-x(mod(-n,N)+1)); %定义傅里叶变换DFT函数 function[Xk]=dft(xn,N) % %计算离散傅里叶变换 %Xk为在0<=n<=N-1间的DFT系数数组 %xn为N点有限持续时间序列 %N为DFT的长度 % n=0:N-1; k=0:N-1; WN=exp(-j*2*pi/N); nk=n'*k; WNnk=WN.^nk; Xk=xn*WNnk; %初始化序列 n=0:10; x=5*(0.4).^n; [xec,xoc]=xhfenjie(x); %绘制计算结果 figure(1) stem(n,xec); title('循环偶分量'); xlabel('n'); ylabel('xec'); axis([-0.5 10.5 -0.5 5.5]); hold; plot([0 10],[0 0]); hold off; figure(2) 5 武汉理工大学《数字信号处理》课程设计说明书 stem(n,xoc); title('循环奇分量'); xlabel('n'); ylabel('xoc'); axis([-0.5 10.5 -2.3 2.3]); hold; plot([0 10],[0 0]); hold off; %对各个序列进行DFT变换 X=dft(x,11); Xec=dft(xec,11); Xoc=dft(xoc,11); %绘制计算结果 figure(1) subplot(2,1,1) stem(n,real(X)); title('输入序列DFT的实部'); xlabel('k'); axis([-0.5 10.5 -0.5 10.5]); hold; plot([0 10],[0 0]); hold off; subplot(2,1,2) stem(n,imag(X)); title('输入序列DFT的虚部'); xlabel('k'); axis([-0.5 10.5 -4 4]); hold; plot([0 10],[0 0]); hold off; figure(2) subplot(2,1,1) stem(n,real(Xec)); title('偶分量的实部'); xlabel('k'); axis([-0.5 10.5 -0.5 10.5]); hold; plot([0 10],[0 0]); hold off; subplot(2,1,2); stem(n,real(Xec)); 6 武汉理工大学《数字信号处理》课程设计说明书 title('偶分量的虚部'); xlabel('k'); axis([-0.5 10.5 -0.5 10.5]); hold; plot([0 10],[0 0]); hold off; figure(3) subplot(2,1,1) stem(n,real(Xoc)); title('奇分量的实部'); xlabel('k'); axis([-0.5 10.5 -15.5e-15 5.5e-15]); hold; plot([0 10],[0 0]); hold off; subplot(2,1,2) stem(n,imag(Xoc)); title('奇分量的虚部'); xlabel('k'); axis([-0.5 10.5 -4 4]); hold; plot([0 10],[0 0]); hold off; 4 结果分析 如下图所示，任一有限长序列，都可以表示成一个共轭对称序列与共轭反对 称序列的和，也称循环偶分量和循环奇分量。其中xec表示x(n)的共轭对称分 量，xoc表示x(n)的共轭反对称分量。 图4-1 循环偶分量示意图 7 武汉理工大学《数字信号处理》课程设计说明书 图4-2 循环奇分量示意图 共轭对称序列满足:x ( n ) = x *( N,n ) , 0 ?n? N ,1 e p e p 共轭反对称序列满足:x (n) =,x*( N,n ) ,0?n?N,1 o p o p 图4-3 输入序列的实部与虚部 上图说明:共轭对称序列与其共轭序列以n=N/2成偶对称，共轭反对称序列与其 共轭序列n=N/2成奇对称。 8 武汉理工大学《数字信号处理》课程设计说明书 图4-4 偶分量的实部与虚部 比较图a和上图中的复序列实数部分的离散傅里叶变换是原来序列离散傅里叶变换的共轭对称分量;复序列虚数部分的离散傅里叶变换是原来序列离散傅里叶变换的共轭反对称分量。 图4-5 奇分量的实部与虚部 复序列共轭对称分量序列的离散傅里叶变换是原来序列离散傅里叶变换的实数部分;复序列共轭对称分量的离散傅里叶变换是原来序列离散傅里叶变换的虚数部分。 9 武汉理工大学《数字信号处理》课程设计说明书 5 对称性的应用 5.1 FFT算法的基本思想 仔细考察DFT与IDFT的运算发现，利用以下两个特性可减少运算量: 1)系数是一个周期函数，它的周期性和对称性可利用来改 进运算，提高计算效率。 例 ，又如 因此 利用这些周期性和对称性，DFT运算中有些项可合并; nk 2)利用W的周期性和对称性，把长度为N点的大点数的DFT运算分解为N 2若干个小点数的DFT。因为DFT的计算量正比于N，N小，计算量也就小。 FFT算法正是基于这样的基本思想发展起来的。它有多种形式，但基本上可 分为两类，时间抽取法和频率抽取法。 按时间抽取的FFT(N点DFT运算的分解): M 先从一个特殊情况开始，假定N是2的整数次方，N=2，M:正整数 将N点的DFT分解为两个N/2点的DFT: 首先将序列x(n)分解为两组，一组为偶数项，一组为奇数项 r=0,1，…，N/2-1 将DFT运算也相应分为两组: (1) 其中X(k)和X(k)分别是x(r)和x(r)的N/2点DFT。 可见，一个N1212 点的DFT可以分解为两个N/2点的DFT，这两个N/2点的DFT再按照上面(1) 10 武汉理工大学《数字信号处理》课程设计说明书 式合成为一个N点DFT，注意到，X(k)，X(k)有N/2个点，即k=0，1，…，12 N/2-1，由(1)式得到X(k)只有N/2点，而实际上X(k)有N个点，即k=0，1，…，N-1，要用X(k)，X(k)表示全部X(K)值，还必须应用系数w的周12 期性和对称性。 快速傅里叶变换FFT是实现DFT的一种快速算法，利用FFT可以大大减少运算量。然而一般讨论的FFT都是基于复数序列的，在实际应用中，经常遇到的是求实数序列的DFT.利用DFT的共轭对称性，用一个N点FFT的程序来计算两个N点实数序列的DFT.这样不仅是一种高效的FFT算法，而且还减少了FFT运算时所必须的存储量和运算时间。具体做法如下: 设x1(n)与x2(n)都为N点实数序列，它们的离散傅里叶变换分别为X1(k)与X2(k).构造一个新的复数序列为 x(n)=x1(n)+jx2(n) 然后利用复数序列FFT程序计算 X(k)=DFT[x(n)]=DFT[x1(n)]+DFT[jx2(n)] =X1(k)+jX2(k) 利用DFT的共轭对称性得 X1(k)=Xep(k)=1/2[X(k)+X*(N-k)] X2(k)=-jXep(k)-1/2[X(k)-X*(N-k)] 5.2 对称性应用的程序实现 程序实现如下: N=8;n=0:N-1;k=0:N-1; x=[0,3,6,5,4,3,2,1]; n1=mod(N-n,N); xep=(x+x(n1+1))/2; xop=(x-x(n1+1))/2; figure(1); subplot(311);stem(n,x);xlabel('n');ylabel('x(n)'); subplot(312);stem(n,xep);xlabel('n');ylabel('xep(n)'); subplot(313);stem(n,xop);xlabel('n');ylabel('xop(n)'); xk=fft(x,N); rexk=real(xk); 11 武汉理工大学《数字信号处理》课程设计说明书 imxk=imag(xk); xn1=ifft(rexk); xn2=ifft(imxk); figure(2); subplot(211);stem(k,abs(rexk));xlabel('k');ylabel('|reX(k)|'); subplot(212);stem(k,abs(imxk));xlabel('k');ylabel('|imX(k)|'); figure(3); subplot(211);stem(n,xn1);xlabel('n');ylabel('xn1'); subplot(212);stem(n,xn2);xlabel('n');ylabel('xn2'); 在MATLAB中仿真，其仿真结果如下图: 图5-1 实数序列的奇偶分量表示 12 武汉理工大学《数字信号处理》课程设计说明书 图5-2 两个实数序列 图5-3 两个序列的DFT变换 13 武汉理工大学《数字信号处理》课程设计说明书 上述程序的思想是:将一个实数序列分解成共轭对称分量和共轭反对称分量，分别表示原序列x(n)的实部和虚部，用快速傅里叶变换FFT求出原序列的傅里叶变换，则傅里叶变换的实部表示的是共轭对称分量的傅里叶变换，虚部表示共轭反对称分量的傅里叶变换，再用逆傅里叶变换求出其时域表达式x1(n)和 x2(n)，即用一次FFT变换实现了两个序列的DFT变换。 14 武汉理工大学《数字信号处理》课程设计说明书 6 心得体会 通过这两个星期数字信号处理的课程设计，我对这门课程有了更进一步的了解。这次课程设计利用了MATLAB编程对傅里叶变换DFT的对称性进行了验证，并利用对称性这一性质说明DFT在计算FFT时可大大减少计算量，用一次FFT即可实现两个序列的FFT运算，从而实现了整体设计。 由于对MATLAB编程还不是很熟悉，所以在编写程序中遇到了很多问题。刚开始在做这个设计时，我先将课本上DFT对称性的相关章节进行了复习，也再次熟悉了MATLAB编程的相关语法。然后通过查阅了网上很多的资料和老师的指导，了解了对称性的原理及其应用，对整个课程设计也有了总体的轮廓，通过学过的知识以及做过的实验，最后有了整体的思路。 通过本次课程设计我学到了很多新的东西，极大地拓宽了我的知识面,锻炼了能力，综合素质也得到较大提高，我感到收获不小。但在设计中也发现了大量问题，有些在设计过程中已经解决，有些还需要在今后慢慢学习，只要学习就会有更多的问题，有更多的难点，但也会有更多的收获。 本次数字信号处理的课程设计，是培养和锻炼在校大学生综合应用所学理论知识解决实际问题能力、进行实际锻炼能力的重要环节，它具有动手、动脑、理论联系实际的特点，是培养学生理论联系实际、敢于动手、善于动手和独立自主解决设计实践中遇到的各种问题能力的一种较好方法。 从这次的课程设计中，我真真正正的意识到，在以后的学习中，要理论联系实际，把我们所学的理论知识用到实际当中，由于此次课设我们用到了单片机的相关知识，虽然我们还未学过相关课程。但我们通过此次课设积极主动去查阅相关资料并向其他同学和学长请教，学到了不少东西，在这个过程中与其他同学积极探讨，培养了我的团队协作精神。这就是我在这次课程设计中的最大收获。 15 武汉理工大学《数字信号处理》课程设计说明书 参考文献 1 陈怀琛 吴大正高西全等编著 MATLAB及在电子信息课程中的应用.第二版. 2003年7月 2 (美)Vinay K.Ingle John G.Proakis编著 陈怀琛 王朝英等译 数字信 号处理及其MATLAB实现 电子工业出版社 2003年5月 3 刘泉 阙大顺编著 数字信号处理原理与实现 电子工业出版社 2005年8 月 4 姚天任 孙洪编著 数字信号处理理论及其应用 华中理工大学出版社 1999年11月 5 程佩青编著 数字信号处理教程 第二版 清华大学出版社 2001年8月 16