傅里叶变换的卷积性质

1、PAGE PAGE 27信号与系统实验报告 题目：Matlab编程大作业 院（系） 电子与信息学院 专 业 信息工程 班 级 （一） 学生姓名 汪涛 座 号 42 提交日期 2007年 09月 03日一、实验环境 Matlab（Matrix Laboratory，又称“矩阵实验室”）是在电子与信息类教学和科研中最常用的一种信号处理和仿真软件。它是MathWorks公司于1984年正式推出的一种科学计算软件。（Matlab 4.2b for Win16） （Matlab 6.5 for Win32） 本次实验在Windows XP环境下进行，同时选用了Matlab 4.2b和6.5两个版本，以便

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

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

4、screte Fourier Transform）。DFT变换的定义为： k=1,2,N j=1,2,N其中可以看出，DFT是对有限长序列频谱的离散化。通过DFT使时域有限长序列与频域有限长序列相对应，从而可在频域用计算机进行信号处理。更重要的是DFT具有两种高效的快速算法FFT（Fast Fourier Transform）。当序列长度为2的整数次幂时，可以采用最快的基2算法，否则采用一般的混合算法。具体的算法可参阅Matlab 6.0或以上版本的Function Reference。一种更为通俗的解释是，DFT是对频谱采样得到的结果，即DFT所求到的X(k)，它只是一组等分频率点上的X(e

5、jw)，通过DFT并没有得到整个频谱。这也是它被称为DFT的原因（Discrete）。本实验要求计算02上X(ejw)和Y(ejw)在n=64的等分频率点上的值，这由DFT就可以得到。值得指出的是，当DFT得到的X(k)中k取遍1,2,N时，在频谱X(ejw)上正好对应了02。由于n=64并不小，所以还可以考虑由X(k)来近似拟合出X(ejw)，这就是能够近似画出在02上画出|X(ejw)|和|Y(ejw)|对的图的理论保证。 通过卷积法同样可以求得系统在某个输入下的输出，这一方法可以用来验证使用DFT算法求输出的正确性。在实验的最后还补充了一个关于两种方法在速度上的简单对比，以此说明DFT的

6、FFT算法的优越性（效率高）。三、实验步骤（包括分析、代码和波形） 首先来看看这个实验的要求。实验主要涉及LTI系统对信号的加工，需要我们通过卷积、FFT和IFFT来灵活计算系统的单位脉冲响应和输出：假设一个特定的LTI系统对输入xn=(-3/4)n的响应是: yn=2/5(1/2)nun+3/5(-3/4)nun要求用fft和ifft确定这个系统的单位脉冲响应.(1).将在区间0nN-1上xn和yn的值存入向量x和y中,这里n=64。画出这两个信号,并确信在这个区间以外,它们基本等于零,因此可以放心截断.(2).用fft计算在区间02内n=64的等分频率点上的X(ej)和Y(ej),并在02

7、上画出|X(ej)|和| Y(ej)|对的图;(3).在相同的64个频率样本点上计算出H(ej),用ifft计算在区间0nN-1内的hn;(4).用计算卷积y1=conv(h1,X)验证hn是正确的. (5).给该系统输入信号xn=n,-2=n=5,使用conv卷积函数，求系统输出。画出输出波形。 实验的思路是很明确的，结合原理中的讨论，我们只要按照题目的要求来编写代码、查看和记录波形并进行验证就可以了，不需要计算。 下面是第(1)小题的代码。这一部分需要注意的是stem函数，它可以用来画离散信号的波形：n=0:1:63; %生成序列的横坐标x=(-3/4).n; %生成序列xn，存放于一个下

8、标为1.64的一维向量中y=(2/5)*(1/2).n+(3/5)*(-3/4).n; %生成序列yn，存放于一个下标为1.64的一维向量中subplot(2,1,1);stem(n,x);xlabel(n);ylabel(xn); %画出xnsubplot(2,1,2);stem(n,y);xlabel(n);ylabel(yn); %画出yn 理解了Matlab中所有变量都是矩阵后，这段代码就更容易接受了。如果遇到不懂的函数，我们只要在Command Window中输入“help 函数名”就可以看到相关的帮助，十分方便。 下面的图给出了xn和yn的波形。 Matlab 4.2b提供了将wm

9、f格式矢量图复制到Windows剪切板的功能，可以将它直接粘贴到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); %求y

10、n的DFT的模subplot(2,1,1);plot(n,magdftx);%通过plot命令将离散序列近似拟合为连续频谱曲线xlabel(n);ylabel(xn Spectrum Magnitude);subplot(2,1,2);plot(n,magdfty);xlabel(n);ylabel(yn Spectrum Magnitude); 得到的波形如下图： 在实验原理部分已经提及，横坐标n从0到63的取值正好对应着频谱当中02上的频率。可以看到，这两个信号中模较大（也可以说是能量较大）的频率分量都集中在高频部分，这同第(1)小题绘制出的信号波形是吻合的。在时域x0和y0点附近，信号x

11、n和yn都呈现出高频的振荡。 第(3)小题的工作是紧接着第(2)小题展开的，H(ejw)正是Y(ejw)除以X(ejw)的结果。在这个DTFT对应的DFT中，dfth也就等于dfty除以dftx。在语句表达上，一定要将除法写成“./”，表示矩阵的对应元素相除。代码如下：dfth=dfty./dftx; %利用傅里叶变换的卷积性质求得系统频率响应的DFTh=ifft(dfth); %求hnmagdfth=abs(dfth); %求hn的DFT的模angdfth=angle(dfth); %求hn的DFT的相位%为了更加形象,把hn和它频谱的模和相位都画出来figure(3);subplot(3,

12、1,1);stem(n,h);xlabel(n);ylabel(hn);subplot(3,1,2);plot(n,magdfth);xlabel(n);ylabel(hn Spectrum Magnitude);subplot(3,1,3);plot(n,angdfth);xlabel(n);ylabel(hn Spectrum Angle); 得到的波形为： 至此我们就求得了系统的的单位脉冲响应hn、H(ejw)，如果此后要求系统在某一输入下的输出，我们就可以直接用hn去卷积输入。这个LTI系统的性质也通过频率响应很好地表现出来。第(4)小题的想法很好，要验证一下我们求到的hn是否正确？否

13、则再利用它来求系统输出可是不可靠的。方法已经说得很明白了，要求我们用hn*xn求一下y1n，看看是否等于yn。为了序列长度统一，我们对y1n进行截断得到y2n，不过被截断的部分是确实为零的，因此不用担心。得到y2n和yn之后，我们还把y2n-yn也画了出来。如果y2n-yn是一个恒为零的信号，就说明y2n和yn是完全相等的了。代码如下：y1=conv(h,x); %求频率响应为hn的系统在输入为xn下的输出y1n%如果y1n=yn,那么就可以验证hn就是正确的.所以把二者画在一起对比看看.%注意到y1共有127列，即y11.127，但是y165到y1127的值都是0.%故先将y1n其截断为和y

14、n等长的序列存于y2n中,以便比较.for i=1:64y2(i)=y1(i);end; %生成y1n的截断y2nfigure(4); %新开一个Figure窗口,将y2n和yn画在一起,再画出y2n-ynsubplot(3,1,1);stem(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-ynsubpl

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

16、括linspace函数的使用，以及如何确定y5n的第一个非零值并画出y5。代码如下：x5=linspace(-2,5,8); %产生信号x5n=n,-2=n=5y5=conv(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(x5n);subplot(3,1,2);stem(n,h),axis(-20,20,0,1);xlabel(n);ylabel(hn);fliph=fliplr(h); %求h-n=fliphn%y5n的起始值就是xn的起始值(-2)+h

17、n的起始值(0)=-2,长度即矩阵本身大小n5=-2:68;subplot(3,1,3);stem(n5,y5),axis(-20,20,-2,10);xlabel(n);ylabel(y5n); 这段代码所绘制出的体现整个卷积过程的波形绘制在下面的图中。在确定n5取值范围的时候，不妨参考一下卷积的这条性质：yn的起始/终止值就是xn和hn的起始/终止值之和。可以参考课本P58-P59。三、实验总结通过本次实验，验证了我们可以通过傅里叶变换的卷积性质将时域繁琐的卷积运算变换到频域后采用DTFT的相乘来简单地完成。在DSP中处理有限长信号（通常对应着FIR系统）时，DTFT又常常表现为它的一个离

18、散采样DFT。DFT使时域有限长序列与频域有限长序列相对应，而且拥有高效的FFT和IFFT算法。实验的基本原理可以用方框图表示为： XkFFTxnIFFT Yk ynFFT hn Hk最后比较一下使用传统的卷积算法和FFT算法在时间效率上的差别。就本实验而言，在Matlab 6.5环境下我们分别抽取了几条FFT/IFFT和卷积命令的代码并用tic/toc命令进行计时，时间如下：运算类型和命令时间(s)时间总和(s)dftx=fft(x); %DFTdfty=fft(y); %DFT0.07800.1090dfth=dfty./dftx;h=ifft(dfth); %IDFT0.0310y1=c

19、onv(h,x);0.15700.1570可以发现，一组FFT/IFFT命令的执行总时间比一条卷积命令的执行时间少了约30.6%。此外，FFT/IFFT命令可以在已知输入和输出的情况下求系统的频率响应，这是卷积无法做到的。要想实现“逆卷积”运算，往往需要考虑一个特殊的输入：sinc函数。具体的论述可在课本P390-P391找到。至此，DFT的原理、正确性和优势已经全部讨论完毕。当然，由DFT的定义可以看到，它的最大局限性在于只能处理时域离散的有限长序列。即使如此，在DSP的工程实践中所遇到的很多问题都是满足这个条件的，所以DFT的应用十分广泛。它的出现为时域数字信号的处理提供了一种新的思路：信

20、号可以转换到频域后进行各种复杂的处理，最后再还原回时域。实验参考书唐向宏，岳恒立，郑雪峰. MATLAB及在电子信息类课程中的应用. 电子工业出版社, 2006.Sanjit K. Miltra. Digital Signal Processing A Computer-based Approach (Third Edition). 清华大学出版社, 2006.A.V. Oppenheim, A.S.Willsky, S.Hamid Nawab著,刘树棠译. Singals and Systems (Second Edition). 西安交通大学出版社, 1998.附录资料：MATLAB的30

21、个方法1 内部常数pi 圆周率 exp(1)自然对数的底数ei 或j 虚数单位Inf或 inf 无穷大 2 数学运算符a+b 加法a-b减法a*b矩阵乘法a.*b数组乘法a/b矩阵右除ab矩阵左除a./b数组右除a.b数组左除ab 矩阵乘方a.b数组乘方-a负号 共轭转置.一般转置3 关系运算符=等于大于=大于或等于=不等于4 常用内部数学函数 指数函数exp(x)以e为底数对数函数log(x)自然对数，即以e为底数的对数log10(x)常用对数，即以10为底数的对数log2(x)以2为底数的x的对数开方函数sqrt(x)表示x的算术平方根绝对值函数abs(x)表示实数的绝对值以及复数的模三角

22、函数（自变量的单位为弧度）sin(x)正弦函数cos(x)余弦函数tan(x)正切函数cot(x)余切函数sec(x)正割函数csc(x)余割函数反三角函数 asin(x)反正弦函数acos(x)反余弦函数atan(x)反正切函数acot(x)反余切函数asec(x)反正割函数acsc(x)反余割函数双曲函数 sinh(x)双曲正弦函数cosh(x)双曲余弦函数tanh(x)双曲正切函数coth(x)双曲余切函数sech(x)双曲正割函数csch(x)双曲余割函数反双曲函数 asinh(x)反双曲正弦函数acosh(x)反双曲余弦函数atanh(x)反双曲正切函数acoth(x)反双曲余切函数

23、asech(x)反双曲正割函数acsch(x)反双曲余割函数求角度函数atan2(y,x)以坐标原点为顶点，x轴正半轴为始边，从原点到点（x，y）的射线为终边的角，其单位为弧度，范围为（ ， 数论函数gcd(a,b)两个整数的最大公约数lcm(a,b)两个整数的最小公倍数排列组合函数factorial(n)阶乘函数，表示n的阶乘 复数函数 real(z)实部函数imag(z)虚部函数abs(z)求复数z的模angle(z)求复数z的辐角，其范围是（ ， conj(z)求复数z的共轭复数求整函数与截尾函数ceil(x)表示大于或等于实数x的最小整数floor(x)表示小于或等于实数x的最大整数r

24、ound(x)最接近x的整数最大、最小函数max(a，b，c，)求最大数min(a，b，c，)求最小数符号函数 sign(x)5 自定义函数-调用时：“返回值列=M文件名（参数列）”function 返回变量=函数名（输入变量） 注释说明语句段（此部分可有可无）函数体语句 6进行函数的复合运算compose(f,g) 返回值为f(g(y)compose(f,g,z) 返回值为f(g(z)compose(f,g,x,.z) 返回值为f(g(z)compose(f,g,x,y,z) 返回值为f(g(z)7 因式分解syms 表达式中包含的变量 factor(表达式) 8 代数式展开syms 表达式

25、中包含的变量 expand(表达式)9 合并同类项syms 表达式中包含的变量 collect(表达式，指定的变量)10 进行数学式化简syms 表达式中包含的变量 simplify(表达式)11 进行变量替换syms 表达式和代换式中包含的所有变量 subs(表达式，要替换的变量或式子，代换式)12 进行数学式的转换调用Maple中数学式的转换命令，调用格式如下：maple(Maple的数学式转换命令) 即：maple(convert(表达式，form)将表达式转换成form的表示方式 maple(convert(表达式，form, x) 指定变量为x，将依赖于变量x的函数转换成form的表

26、示方式（此指令仅对form为exp与sincos的转换式有用） 13 解方程solve(方程，变元) 注：方程的等号用普通的等号： = 14 解不等式调用maple中解不等式的命令即可，调用形式如下： maple(maple中解不等式的命令)具体说，包括以下五种：maple( solve（不等式）) maple( solve（不等式，变元） ) maple( solve（不等式，变元） ) maple( solve（不等式，变元） ) maple( solve（不等式，变元） )15 解不等式组调用maple中解不等式组的命令即可，调用形式如下： maple(maple中解不等式组的命令) 即

27、：maple( solve（不等式组，变元组） )16 画图方法：先产生横坐标的取值和相应的纵坐标的取值，然后执行命令： plot(x,y) 方法2：fplot(f(x)，xmin,xmax) fplot(f(x)，xmin,xmax,ymin,ymax) 方法3：ezplot(f(x) ezplot(f(x) ,xmin,xmax) ezplot(f(x) ,xmin,xmax,ymin,ymax) 17 求极限（1）极限：syms x limit(f(x), x, a) （2）单侧极限：左极限：syms x limit(f(x), x, a,left)右极限：syms x limit(f(

28、x), x, a,right) 18 求导数diff(f(x) diff(f(x),x) 或者：syms x diff(f(x) syms x diff(f(x), x) 19 求高阶导数 diff(f(x)，n) diff(f(x),x，n)或者：syms x diff(f(x)，n)syms x diff(f(x), x，n) 20 在MATLAB中没有直接求隐函数导数的命令，但是我们可以根据数学中求隐函数导数的方法，在中一步一步地进行推导；也可以自己编一个求隐函数导数的小程序；不过，最简便的方法是调用Maple中求隐函数导数的命令，调用格式如下： maple(implicitdiff(f(x,y)=0,y,x) 在MATLAB中，没有直接求参数方程确定的函数的导数的命令，只能根据参数方程确定的函数的求导公式 一步一步地进行推导；或者，干脆自己编一个小程序，应用起来会更加方便。21 求不定积分 int(f(x) int (f(x),x)或者：syms x int(f(x) syms x int(f(x), x