傅里叶变换的卷积性质.docx

傅里叶变换的卷积性质.docx

《傅里叶变换的卷积性质.docx》由会员分享，可在线阅读，更多相关《傅里叶变换的卷积性质.docx（15页珍藏版）》请在冰点文库上搜索。

傅里叶变换的卷积性质

信号与系统实验报告

题目：

Matlab编程大作业

院（系）电子与信息学院

专业信息工程

班级

（一）

学生姓名汪涛

座号42

提交日期2007年09月03日

一、实验环境

Matlab（MatrixLaboratory，又称“矩阵实验室”）是在电子与信息类教学和科研中最常用的一种信号处理和仿真软件。

它是MathWorks公司于1984年正式推出的一种科学计算软件。

（Matlab4.2bforWin16）

（Matlab6.5forWin32）

本次实验在WindowsXP环境下进行，同时选用了Matlab4.2b和6.5两个版本，以便对这款软件十多年的开发、完善过程有一个直观的认识。

两个版本发布于1994年和2002年，它们分别是16位和32位Windows系统上最为风靡的Matlab版本。

Matlab6.x的一个重要新特性是LaunchPad，它为用户打开工具箱提供了方便，其中包括了帮助文件、演示范例、实用工具和Web文档等等。

经过简单的使用后，我发现6.5相对于4.2b在易用性方面的优势包括：

（1）Workspace可以让用户方便地看到当前的工作空间，而4.2b只能通过whos等命令查看，占用了CommandWindow，不方便；

（2）6.5内置更多常用的函数，比如y=sinc（x），在4.2b中只能通过手工定义产生；

（3）6.5的帮助文档比4.2b更详尽，而且图文并茂、易于阅读。

而4.2b的优势也是十分明显的，4.2b版本软件体积小而基本功能齐备、占用内存小、运行速度快、界面简洁明了、易于上手。

本次实验主要在4.2b环境中完成，同时也参考了6.5的帮助文档。

二、实验原理

在离散时间、连续频率的傅里叶变换中，由于卷积性质知道，对系统输出的计算可以通过求x[n]和h[n]的DTFT，将得到的X（ejw）和H（ejw）相乘就可以得到Y（ejw），进而再通过反变换得到y[n]。

这就避免了在时域进行繁琐的卷积求解。

在数字信号处理（DSP）中，对于有限长序列存在一种离散时间、离散频率的傅里叶变换，称为DFT（DiscreteFourierTransform）。

DFT变换的定义为：

k=1,2,…,N

j=1,2,…,N

其中

可以看出，DFT是对有限长序列频谱的离散化。

通过DFT使时域有限长序列与频域有限长序列相对应，从而可在频域用计算机进行信号处理。

更重要的是DFT具有两种高效的快速算法FFT（FastFourierTransform）。

当序列长度为2的整数次幂时，可以采用最快的基2算法，否则采用一般的混合算法。

具体的算法可参阅Matlab6.0或以上版本的FunctionReference。

一种更为通俗的解释是，DFT是对频谱采样得到的结果，即DFT所求到的X（k），它只是一组等分频率点上的X（ejw），通过DFT并没有得到整个频谱。

这也是它被称为DFT的原因（Discrete）。

本实验要求计算0≤ω<2π上X（ejw）和Y（ejw）在n=64的等分频率点上的值，这由DFT就可以得到。

值得指出的是，当DFT得到的X（k）中k取遍1,2,…,N时，在频谱X（ejw）上正好对应了0≤ω<2π。

由于n=64并不小，所以还可以考虑由X（k）来近似拟合出X（ejw），这就是能够近似画出在0≤ω≤2π上画出|X（ejw）|和

|Y（ejw）|对ω的图的理论保证。

通过卷积法同样可以求得系统在某个输入下的输出，这一方法可以用来验证使用DFT算法求输出的正确性。

在实验的最后还补充了一个关于两种方法在速度上的简单对比，以此说明DFT的FFT算法的优越性（效率高）。

三、实验步骤（包括分析、代码和波形）

首先来看看这个实验的要求。

实验主要涉及LTI系统对信号的加工，需要我们通过卷积、FFT和IFFT来灵活计算系统的单位脉冲响应和输出：

假设一个特定的LTI系统对输入x[n]=（-3/4）n的响应是:

y[n]=2/5（1/2）nu[n]+3/5（-3/4）nu[n]

要求用fft和ifft确定这个系统的单位脉冲响应.

（1）.将在区间0≤n≤N-1上x[n]和y[n]的值存入向量x和y中,这里n=64。

画出这两个信号,并确信在这个区间以外,它们基本等于零,因此可以放心截断.

（2）.用fft计算在区间0≤ω<2π内n=64的等分频率点上的X（ejω）和Y（ejω）,并在0≤ω≤2π上画出|X（ejω）|和|Y（ejω）|对ω的图;

（3）.在相同的64个频率样本点上计算出H（ejω）,用ifft计算在区间0≤n≤N-1内的h[n];

（4）.用计算卷积y1=conv（h1,X）验证h[n]是正确的.

（5）.给该系统输入信号x[n]=n,-2<=n<=5,使用conv卷积函数，求系统输出。

画出输出波形。

实验的思路是很明确的，结合原理中的讨论，我们只要按照题目的要求来编写代码、查看和记录波形并进行验证就可以了，不需要计算。

下面是第

（1）小题的代码。

这一部分需要注意的是stem函数，它可以用来画离散信号的波形：

n=0:

1:

63;%生成序列的横坐标

x=（-3/4）.^n;%生成序列x[n]，存放于一个下标为[1..64]的一维向量中

y=（2/5）*（1/2）.^n+（3/5）*（-3/4）.^n;%生成序列y[n]，存放于一个下标为[1..64]的一维向量中

subplot（2,1,1）;stem（n,x）;

xlabel（'n'）;ylabel（'x[n]'）;%画出x[n]

subplot（2,1,2）;stem（n,y）;

xlabel（'n'）;ylabel（'y[n]'）;%画出y[n]

理解了Matlab中所有变量都是矩阵后，这段代码就更容易接受了。

如果遇到不懂的函数，我们只要在CommandWindow中输入“help函数名”就可以看到相关的帮助，十分方便。

下面的图给出了x[n]和y[n]的波形。

Matlab4.2b提供了将wmf格式矢量图复制到Windows剪切板的功能，可以将它直接粘贴到Word文档中而不发生任何失真，十分好用。

下面是第

（2）小题的代码。

应该注意到abs和angle函数可以用来分解出复数的模和幅角。

另外，plot函数可以将离散序列画成连续曲线：

dftx=fft（x）;%求x[n]的DFT

dfty=fft（y）;%求y[n]的DFT

%特别注意:

DFT算出的只是对连续频谱的64个采样值,而非整个频谱

%但是可以利用这些采样值来近似拟合出频谱,这就是下面的工作了

figure

（2）;%新开一个Figure窗口

magdftx=abs（dftx）;%求x[n]的DFT的模

magdfty=abs（dfty）;%求y[n]的DFT的模

subplot（2,1,1）;plot（n,magdftx）;

%通过plot命令将离散序列近似拟合为连续频谱曲线

xlabel（'n'）;ylabel（'x[n]SpectrumMagnitude'）;

subplot（2,1,2）;plot（n,magdfty）;

xlabel（'n'）;ylabel（'y[n]SpectrumMagnitude'）;

得到的波形如下图：

在实验原理部分已经提及，横坐标n从0到63的取值正好对应着频谱当中0≤ω<2π上的频率。

可以看到，这两个信号中模较大（也可以说是能量较大）的频率分量都集中在高频部分，这同第

（1）小题绘制出的信号波形是吻合的。

在时域x[0]和y[0]点附近，信号x[n]和y[n]都呈现出高频的振荡。

第（3）小题的工作是紧接着第

（2）小题展开的，H（ejw）正是Y（ejw）除以X（ejw）的结果。

在这个DTFT对应的DFT中，dfth也就等于dfty除以dftx。

在语句表达上，一定要将除法写成“./”，表示矩阵的对应元素相除。

代码如下：

dfth=dfty./dftx;%利用傅里叶变换的卷积性质求得系统频率响应的DFT

h=ifft（dfth）;%求h[n]

magdfth=abs（dfth）;%求h[n]的DFT的模

angdfth=angle（dfth）;%求h[n]的DFT的相位

%为了更加形象,把h[n]和它频谱的模和相位都画出来

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]SpectrumMagnitude'）;

subplot（3,1,3）;plot（n,angdfth）;

xlabel（'n'）;ylabel（'h[n]SpectrumAngle'）;

得到的波形为：

至此我们就求得了系统的的单位脉冲响应h[n]、H（ejw），如果此后要求系统在某一输入下的输出，我们就可以直接用h[n]去卷积输入。

这个LTI系统的性质也通过频率响应很好地表现出来。

第（4）小题的想法很好，要验证一下我们求到的h[n]是否正确？

否则再利用它来求系统输出可是不可靠的。

方法已经说得很明白了，要求我们用h[n]*x[n]求一下y1[n]，看看是否等于y[n]。

为了序列长度统一，我们对y1[n]进行截断得到y2[n]，不过被截断的部分是确实为零的，因此不用担心。

得到y2[n]和y[n]之后，我们还把y2[n]-y[n]也画了出来。

如果y2[n]-y[n]是一个恒为零的信号，就说明y2[n]和y[n]是完全相等的了。

代码如下：

y1=conv（h,x）;%求频率响应为h[n]的系统在输入为x[n]下的输出y1[n]

%如果y1[n]=y[n],那么就可以验证h[n]就是正确的.所以把二者画在一起对比看看.

%注意到y1共有127列，即y1[1..127]，但是y1[65]到y1[127]的值都是0.

%故先将y1[n]其截断为和y[n]等长的序列存于y2[n]中,以便比较.

fori=1:

64

y2（i）=y1（i）;

end;%生成y1[n]的截断y2[n]

figure（4）;%新开一个Figure窗口,将y2[n]和y[n]画在一起,再画出y2[n]-y[n]

subplot（3,1,1）;stem（n,y2）,axis（[0,70,0,1]）;%指定坐标的最小、最大值,方便统一观察

xlabel（'n'）;ylabel（'y2[n]'）;

subplot（3,1,2）;stem（n,y）,axis（[0,70,0,1]）;%指定坐标的最小、最大值,方便统一观察

xlabel（'n'）;ylabel（'y[n]'）;

chk=y2-y;%生成检查信号chk[n]=y2[n]-y[n]

subplot（3,1,3）;stem（n,chk）,axis（[0,70,0,1]）;

xlabel（'n'）;ylabel（'y2[n]-y[n]'）;

%如果y2[n]-y[n]是一个恒等于0的信号,说明y2[n]和y[n]相等,事实上也确实如此.

%这就验证了使用fft和ifft函数来计算频率响应是正确的.

这段代码所绘制出的波形在下一页的图中。

由于y2[n]-y[n]画出来确实是一个恒等于0的信号，我们也可以通过在CommandWindow中输入“chk”来查看这个矩阵的元素确实都是0，所以现在我们可以放心地使用h[n]来求输出了。

第（5）小题就提出了用h[n]来求输出的问题，这一次不要求我们用DFT，而用最传统的卷积法。

代码中值得注意的技巧包括linspace函数的使用，以及如何确定y5[n]的第一个非零值并画出y[5]。

代码如下：

x5=linspace（-2,5,8）;%产生信号x5[n]=n,-2<=n<=5

y5=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（'x5[n]'）;

subplot（3,1,2）;stem（n,h）,axis（[-20,20,0,1]）;

xlabel（'n'）;ylabel（'h[n]'）;

fliph=fliplr（h）;%求h[-n]=fliph[n]

%y5[n]的起始值就是x[n]的起始值（-2）+h[n]的起始值（0）=-2,长度即矩阵本身大小

n5=-2:

68;

subplot（3,1,3）;stem（n5,y5）,axis（[-20,20,-2,10]）;

xlabel（'n'）;ylabel（'y5[n]'）;

这段代码所绘制出的体现整个卷积过程的波形绘制在下面的图中。

在确定n5取值范围的时候，不妨参考一下卷积的这条性质：

y[n]的起始/终止值就是x[n]和h[n]的起始/终止值之和。

可以参考课本P58-P59。

三、实验总结

通过本次实验，验证了我们可以通过傅里叶变换的卷积性质将时域繁琐的卷积运算变换到频域后采用DTFT的相乘来简单地完成。

在DSP中处理有限长信号（通常对应着FIR系统）时，DTFT又常常表现为它的一个离散采样——DFT。

DFT使时域有限长序列与频域有限长序列相对应，而且拥有高效的FFT和IFFT算法。

实验的基本原理可以用方框图表示为：

X[k]

x[n]

Y[k]y[n]

h[n]H[k]

最后比较一下使用传统的卷积算法和FFT算法在时间效率上的差别。

就本实验而言，在Matlab6.5环境下我们分别抽取了几条FFT/IFFT和卷积命令的代码并用tic/toc命令进行计时，时间如下：

运算类型和命令

时间（s）

时间总和（s）

dftx=fft（x）;%DFT

dfty=fft（y）;%DFT

0.0780

0.1090

dfth=dfty./dftx;

h=ifft（dfth）;%IDFT

0.0310

y1=conv（h,x）;

0.1570

0.1570

可以发现，一组FFT/IFFT命令的执行总时间比一条卷积命令的执行时间少了约30.6%。

此外，FFT/IFFT命令可以在已知输入和输出的情况下求系统的频率响应，这是卷积无法做到的。

要想实现“逆卷积”运算，往往需要考虑一个特殊的输入：

sinc函数。

具体的论述可在课本P390-P391找到。

至此，DFT的原理、正确性和优势已经全部讨论完毕。

当然，由DFT的定义可以看到，它的最大局限性在于只能处理时域离散的有限长序列。

即使如此，在DSP的工程实践中所遇到的很多问题都是满足这个条件的，所以DFT的应用十分广泛。

它的出现为时域数字信号的处理提供了一种新的思路：

信号可以转换到频域后进行各种复杂的处理，最后再还原回时域。

[实验参考书]

唐向宏，岳恒立，郑雪峰.MATLAB及在电子信息类课程中的应用.电子工业出版社,2006.

SanjitK.Miltra.DigitalSignalProcessing–AComputer-basedApproach（ThirdEdition）.清华大学出版社,2006.

A.V.Oppenheim,A.S.Willsky,S.HamidNawab著,刘树棠译.SingalsandSystems（SecondEdition）.西安交通大学出版社,1998.