傅里叶变换的卷积性质

1、#詢理N久穿South China University of Technology信号与系统实验报告题目：Matlab编程大作业院（系）电子与信息学院专 业信息工程班 级（一）学生姓名汪涛座 号 42提交日期 2007年09月03日ii3、实验环境Matlab （Matrix Laboratory，又称“矩阵实验室”）是在电子与信息类教学和 科研中最常用的一种信号处理和仿真软件。它是MathWorks公司于1984年正式推出的一种科学计算软件。About MATLABMATLABJ hf LatignuRC n/ lechnicitl (. orup ttfimCopjnlght 71984

2、-1994 The MatihWorks, Inc.FI3TR0SRESERWD0S4 2bMay 26,1994Cracked by BurJassailINGENIE RIA INFOR MATICASerial Number: 666(Matlab 4.2b for Win 16)Version 6.5.01 8091 3a Release 13June fB, 2002License Number Q Wee anvvecanl InternationalShow Licenseli 二巫二 X(Matlab 6.5 for Win 32)本次实验在 Windows XP环境下进行，同

3、时选用了 Matlab 4.2b和6.5两个版 本，以便对这款软件十多年的开发、完善过程有一个直观的认识。两个版本发布 于1994年和2002年，它们分别是16位和32位Windows系统上最为风靡的Matlab 版本。Matlab 6.x的一个重要新特性是Launch Pad,它为用户打开工具箱提供了方 便，其中包括了帮助文件、演示范例、实用工具和Web文档等等。经过简单的使用后，我发现6.5相对于4.2b在易用性方面的优势包括：（1） Workspace可以让用户方便地看到当前的工作空间，而4.2b只能通过 whos等命令查看，占用了 Comma nd Win dow，不方便；（2）6.5

4、内置更多常用的函数，比如 y=sinc（x），在4.2b中只能通过手工定 义产生；（3）6.5的帮助文档比4.2b更详尽，而且图文并茂、易于阅读。而4.2b的优势也是十分明显的，4.2b版本软件体积小而基本功能齐备、占 用内存小、运行速度快、界面简洁明了、易于上手。本次实验主要在4.2b环境中完成，同时也参考了 6.5的帮助文档。二、实验原理在离散时间、连续频率的傅里叶变换中，由于卷积性质知道，对系统输出的 计算可以通过求xn和hn的DTFT，将得到的X(ejw)和H(ejw)相乘就可以得到 丫(ejw)，进而再通过反变换得到yn。这就避免了在时域进行繁琐的卷积求解。在数字信号处理(DSP)中

5、，对于有限长序列存在一种离散时间、离散频率)。DFT变换的定义为:的傅里叶变换，称为 DFT( Discrete Fourier Transform口 r 4 小匚'k=1,2,NJ = 1亠班力=(1/N) Y三 入j=1,2,N其中e )/N可以看出，DFT是对有限长序列频谱的离散化。 通过DFT使时域有限长序 列与频域有限长序列相对应，从而可在频域用计算机进行信号处理。更重要的是 DFT具有两种高效的快速算法 FFT (Fast Fourier Transform)。当序列长度为2的 整数次幕时，可以采用最快的基 2算法，否则采用一般的混合算法。具体的算 法可参阅Matlab 6

6、.0或以上版本的Function Referenc一种更为通俗的解释是，DFT是对频谱采样得到的结果，即 DFT所求到的 X(k)，它只是一组等分频率点上的 X(ejw)，通过DFT并没有得到整个频谱。这 也是它被称为DFT的原因(Discrete) 0本实验要求计算0W 3 <2 n上X(ejw)和Y(ejw)在n=64的等分频率点上的值， 这由DFT就可以得到。值得指出的是，当DFT得到的X(k)中k取遍1,2,N时, 在频谱X(ejw)上正好对应了 0W® <2n。由于n=64并不小，所以还可以考虑 由 X(k)来近似拟合出X(ejw)，这就是能够近似画出在 0W

7、3 W 2 n上画出|X(ejw)|和 |Y(ejw)|对3的图的理论保证。通过卷积法同样可以求得系统在某个输入下的输出，这一方法可以用来验证使用DFT算法求输出的正确性。在实验的最后还补充了一个关于两种方法在速 度上的简单对比，以此说明 DFT的FFT算法的优越性(效率高)。三、实验步骤(包括分析、代码和波形)首先来看看这个实验的要求。 实验主要涉及 LTI 系统对信号的加工， 需要我 们通过卷积、 FFT 和 IFFT 来灵活计算系统的单位脉冲响应和输出：假设一个特定的LTI系统对输入xn=(-3/4)n的响应是：yn=2/5(1/2)nun+3/5(-3/4)nun要求用 fft 和 i

8、fft 确定这个系统的单位脉冲响应 .(1) .将在区间Ow nWN1上xn和yn的值存入向量x和y中,这里n=64。 画出这两个信号 , 并确信在这个区间以外 , 它们基本等于零 , 因此可以放心截断 .(2) .用fft计算在区间0=3 <2n内 n=64的等分频率点上的X(e j。和Y(e j 并在0W3W 2上画出|X(e jw)|和| Y(e j 3)|对3的图；(3) . 在相同的 64 个频率样本点上计算出 H(e j 3), 用 ifft 计算在区间0w nwN-1 内的 hn；(4) . 用计算卷积 y1=conv(h1,X) 验证 hn 是正确的 .(5) . 给该系

9、统输入信号 xn=n,-2<=n<=5, 使用 conv 卷积函数，求系统输 出。画出输出波形。实验的思路是很明确的， 结合原理中的讨论， 我们只要按照题目的要求来编 写代码、查看和记录波形并进行验证就可以了，不需要计算。下面是第(1)小题的代码。这一部分需要注意的是stem函数，它可以用来画离散信号的波形：，存放于一个下标为 1.64 的一维向量中%生成序列 yn ，存放于一个下标为 1.64画出 xn画出 ynn=0：1：63； % 生成序列的横坐标 x=(-3/4).An; % 生成序列 xn y=(2/5)*(1/2).A n+(3/5)*(-3/4).A n;的一维向量中

10、subplot(2,1,1);stem(n,x); xlabel('n');ylabel('xn'); % subplot(2,1,2);stem(n,y); xlabel('n');ylabel('yn'); %理解了 Matlab 中所有变量都是矩阵后，这段代码就更容易接受了。如果遇 到不懂的函数，我们只要在 Comma nd Win dow中输入“ help函数名”就可以看 到相关的帮助，十分方便。下面的图给出了 xn和yn的波形。Matlab 4.2b提供了将wmf格式矢量图复制到 Windows剪切板的功能，可以 将它直

11、接粘贴到 Word文档中而不发生任何失真，十分好用。下面是第(2)小题的代码。应该注意到abs和angle函数可以用来分解出复数 的模和幅角。另外，plot函数可以将离散序列画成连续曲线：dftx=fft(x); %求 xn 的 DFTdfty=fft(y); %求 yn的 DFT%持别注意:DFT算出的只是对连续频谱的64个采样值，而非整个频谱%但是可以利用这些采样值来近似拟合出频谱，这就是下面的工作了figure(2); %新开一个 Figure 窗口magdftx=abs(dftx); %求 xn的 DFT 的模magdfty=abs(dfty); %求 yn的 DFT 的模subplo

12、t(2,1,1);plot( n,magdftx);%通过plot命令将离散序列近似拟合为连续频谱曲线xlabel(' n');ylabel('x n Spectrum Mag ni tude');subplot(2,1,2);plot( n,magdfty);xlabel(' n');ylabel('y n Spectrum Mag ni tude');得到的波形如下图:在实验原理部分已经提及，横坐标n从0到63的取值正好对应着频谱当中 OW® <2n上的频率。可以看到，这两个信号中模较大(也可以说是能量较大)

13、的频率分量都集中在高频部分，这同第 (1)小题绘制出的信号波形是吻合的。在 时域x0和y0点附近，信号xn和yn都呈现出高频的振荡。第(3)小题的工作是紧接着第(2)小题展开的，H(ejw)正是丫(ejw)除以X(ejw)的 结果。在这个DTFT对应的DFT中，dfth也就等于dfty除以dftx。在语句表达 上，一定要将除法写成“ /”，表示矩阵的对应元素相除。代码如下：dfth=dfty./dftx; %利用傅里叶变换的卷积性质求得系统频率响应的DFTh=ifft(dfth); %求 hnmagdfth=abs(dfth); %求 hn的 DFT 的模angdfth=angle(dfth)

14、; %求 hn 的 DFT 的相位%为了更加形象，把hn和它频谱的模和相位都画出来figure(3);subplot(3,1,1);stem( n,h);xlabel(' n');ylabel('h n');subplot(3,1,2);plot( n,magdfth);xlabel(' n');ylabel('h n Spectrum Mag ni tude');subplot(3,1,3);plot( n,an gdfth);xlabel(' n');ylabel('h n Spectrum An gl

15、e');得到的波形为:至此我们就求得了系统的的单位脉冲响应 hn、H(ejw)，如果此后要求系统 在某一输入下的输出，我们就可以直接用hn去卷积输入。这个LTI系统的性质 也通过频率响应很好地表现出来。第(4)小题的想法很好，要验证一下我们求到的hn是否正确？否则再利用它 来求系统输出可是不可靠的。方法已经说得很明白了，要求我们用hn*xn求一 下y1n，看看是否等于yn。为了序列长度统一，我们对y1n进行截断得到y2n，不过被截断的部分是确实为零的，因此不用担心。得到y2n和yn之后，我们还把y2n-yn也画了出来。如果y2n-yn是一个恒为零的信号，就说明y2n和 yn 是完全相等

16、的了。代码如下：y1=conv(h,x); % 求频率响应为 hn 的系统在输入为 xn 下的输出 y1n%如果 y1n=yn, 那么就可以验证 hn 就是正确的 . 所以把二者画在一起对比看看 . %注意到 y1 共有 127 列，即 y11.127，但是 y165 到 y1127 的值都是 0.%故先将 y1n 其截断为和 yn 等长的序列存于 y2n中 , 以便比较 .for i=1:64y2(i)=y1(i);end; % 生成 y1n 的截断 y2nfigure(4);嘶开一个 Figure 窗口，将 y2n和 yn画在一起，再画出 y2n-ynsubplot(3,1,1);stem

17、(n,y2),axis(0,70,0,1); %指定坐标的最小、最大值 , 方便统一观察xlabel('n');ylabel('y2n');subplot(3，1，2);stem(n，y)，axis(0，70，0，1); %指定坐标的最小、最大值 ， 方便统一观察xlabel('n');ylabel('yn');chk=y2-y; % 生成检查信号 chkn=y2n-ynsubplot(3，1，3);stem(n，chk)，axis(0，70，0，1);xlabel('n');ylabel('y2n-yn&

18、#39;);%如果 y2n-yn 是一个恒等于 0 的信号 ， 说明 y2n 和 yn 相等 ， 事实上也确实如此 . %这就验证了使用 fft 和 ifft 函数来计算频率响应是正确的 .这段代码所绘制出的波形在下一页的图中。由于y2n-yn画出来确实是一个恒等于0的信号，我们也可以通过在Comma nd Window中输入“ chk”来查看这个矩阵的元素确实都是0，所以现在我们可以放心地使用h n来求输出了。第(5)小题就提出了用hn来求输出的问题，这一次不要求我们用DFT，而用最传统的卷积法。代码中值得注意的技巧包括lin space函数的使用，以及如何确定y5n的第一个非零值并画出y5

19、。代码如下：x5=linspace(-2,5,8); %产生信号 x5n=n,-2<=*=5y5=c on v(h,x5);figure(5);subplot(3,1,1); stem(-2,-1,0,1,2,3,4,5,x5),axis(-20,20,-2,5);xlabel(' n');ylabel('x5 n');subplot(3,1,2);stem( n,h),axis(-20,20,0,1);xlabel(' n');ylabel('h n');fliph=fliplr(h); %求 h-n=fliphn%y5n

20、的起始值就是xn的起始值(-2)+hn的起始值(0)=-2, 长度即矩阵本身大小n5=-2:68;subplot(3,1,3);stem( n5,y5),axis(-20,20,-2,10);xlabel(' n');ylabel('y5 n');这段代码所绘制出的体现整个卷积过程的波形绘制在下面的图中。在确定 n5取值范围的时候，不妨参考一下卷积的这条性质：yn的起始/终止值就是xn 和hn的起始/终止值之和。可以参考课本P58-P59。三、实验总结通过本次实验，验证了我们可以通过傅里叶变换的卷积性质将时域繁琐的卷 积运算变换到频域后采用DTFT的相乘来简单地完成。在DSP中处理有限长信 号(通常对应着FIR系统)时，DTFT又常常表现为它的一个离散采样一一 DFT。 DFT使时域有限长序列与频域有限长序列相对应，而且拥有高效的FFT和IFFT算法。实验的基本原理可以用方框图表示为：Xk丫kynIFFT最后比较一下使用传统的卷积算法和 FFT算法在时间效率上的差别。就本 实验而言，在Matlab 6.5环境下我们分别抽取了几条FFT/IFFT和卷积命令的代 码并用tic/toc命令进行计时，时间如下：运算类型和命令时间(s)时间总和(s)dftx=fft