快速傅立叶变换FFT及其应用

一、实验目的

１. 利用MATLAB的快速傅立叶变换来计算信号的离散傅立叶变换。

２. 利用MATLAB程序，理解进一步离散傅立叶变换的物理意义。

3. 利用MATLAB程序，理解快速卷积算法。

二、实验原理

在MATLAB中，使用函数fft可以很容易地计算有限长序列x(n)的离散傅立叶变换X[k]。此函数有两种形式，fft(x)计算序列x(n) 的离散傅立叶变换X(k)，这里X(k)的长度与x(n)的长度相等。fft(x，L)计算序列x(n) 的L点离散傅立叶变换，其中L≥N。若L>N，在计算离散傅立叶变换之前，对x(n)尾部的L-N个值进行补零。同样，离散傅立叶变换序列X(k)的离散傅立叶逆变换x(n)用函数ifft计算，它也有两种形式。

（一）、基本序列的离散傅立叶变换计算

N点离散傅立叶变换的一种物理解释就是，X[k]是x(n)以N为周期的周期延拓序列的离散傅立叶级数系数的主值区间序列，即。例如序列,当N=16时，正好是的一个周期，所以的周期延拓序列就是这种单一频率的正弦序列。而当N=8时，正好是的半个周期，的周期延拓就不再是单一频率的正弦序列，而是含有丰富的谐波成分，其离散傅立叶级数的系数与N=16时的差别很大，因此对信号进行谱分析时，一定要截取整个周期，否则得到错误的频谱。

（二）、验证N点DFT的物理意义

假如x(n)非周期、有限长，则傅立叶变换存在，那么对在N个等间隔频率=2πk/N, k=0,1,…,,N-1取样,则可得X(k)

          序列x(n)的N点DFT的物理意义是对X(ω)在[0，2π]上进行N点的等间隔采样。

（三）、利用FFT计算序列的线性卷积

    直接计算线性卷积计算量大，并且计算机无法判断y(n)的长度，需要计算多少的y(n)值，若输入为无限长，就更无法计算，其运算量随长度成级数增长。由于可以利用FFT对DFT进行有效的计算，我们希望能够利用DFT来计算线性卷积。

设 x(n) 和 h(n) 是长度分别为M和N的有限长序列，

令 L=M+N-1，定义两个长度L的有限长序列：

 通过对x(n) 和 h(n)补充零样本值得到上面两个序列。那么：

                           (3.4.10)

上面的过程如下图所示：

计算线性卷积也可以直接调用函数con来计算，因为MATLAB中的计时比较粗糙，所以只有M和N较大的时候，才能比较两种方法的执行时间快慢。

三、实验内容与步骤

１. 对复正弦序列，利用MATLAB程序求当N=16和N=8时的离散傅立叶变换，并显示其图形。

程序：

       N=16；N1=8；

       n=0:N-1;k=0: N1-1;

       x=exp(j*pi*n/8);

       X1=fft(x,N);

       X2=fft(x,N1);

       subplot(2,1,1);

       stem(n,abs(X1));

       axis([0,20,0,20]);

       ylabel(‘∣X1(k)∣’);

       title(‘16点的DFT’)；

subplt(2,1, 2);

       stem(n, abs(X2));

       axis([0, 20, 0, 20]);

       ylabel(‘∣X2(k)∣’);

       title(‘8点的DFT’)；

２.已知,, 绘制相应的幅频和相频曲线，并计算N=8和N=16时的DFT。

  程序：

       N1=8;N2=16;

       n=0:N1-1;k1=0:N1-1;k2=0:N2-1;

       w=2*pi*(0:2047)/2048;

Xw=(1-exp(-j*4*w))./(1-exp(-j*));%对x(n)的频谱函数采样2048个点可以近似地看作是连续的频谱

xn=[(n>=0)&(n<4)];%产生x(n)

X1k=fft(xn,N1);

X2k=fft(xn,N2);

subplot(3,2,1);plot(w/pi,abs(Xw));xlabel(‘w/π’);

title(‘x(n)的幅频曲线’)；

subplot(3,2,2);plot(w/pi,angle(Xw));axis([0,2,-pi,pi]);

line([0,2],[0,0]);xlabel(‘w/π’);

title(‘x(n)的相频曲线’)；

subplot(3,2,3);stem(k1,abs(X1k),’.’);

xlabel(‘k(w=2πk/N1)’);ylabel(‘∣X1(k)∣’);hold on

plot(N1/2*w/pi,abs(Xw));%图形上叠加连续频谱的幅度曲线

subplot(3,2,4);stem(k1,angle(X1k)，’.’);

axis([0,N1,-pi,pi]);line([0,N1],[0,0]);

xlabel(‘k(w=2πk/N1)’);ylabel(‘Arg[X1(k)]’);hold on

plot(N1/2*w/pi,angle(Xw));%图形上叠加连续频谱的相位曲线

subplot(3,2,5);stem(k2,abs(X2k),’.’);

xlabel(‘k(w=2πk/N2)’);ylabel(‘∣X2(k)∣’);hold on

plot(N2/2*w/pi,abs(Xw));%图形上叠加连续频谱的幅度曲线

subplot(3,2,6);stem(k2,angle(X2k)，’.’);

axis([0,N2,-pi,pi]);line([0,N2],[0,0]);

xlabel(‘k(w=2πk/N2)’);ylabel(‘Ang[X2(k)]’);hold on

plot(N2/2*w/pi,angle(Xw));%图形上叠加连续频谱的相位曲线

3．分别利用快速卷积法以及conv函数计算下面两个序列的线性卷积。

h=[3  2  1  -2   1  0  -4  0  3];

              ↑

x=[1  -2  3  -4  3  2  1 ]

       ↑

程序1：快速卷积

           clf;

           h=[3  2  1  -2   1  0  -4  0  3];%冲激

           x=[1  -2  3  -4  3  2  1 ]; %输入序列

           L=pow2(nextpow2(length(x)+length(h)-1));

           Xk=fft(x,L);

           Hk=fft(h,L);

           Yk=Xk.*Hk;

           y=ifft(Yk,L);

           nh=0:8;nx=0:6;ny=0:L-1;

           subplot(3,1,1);stem(nx,x);title(‘x(n)’);

       subplot(3,1,2);stem(nh,h);title(‘h(n)’);

       subplot(3,1,3);stem(nx,x); xlabel(‘时间序号n’);ylabel(‘振幅’)；title(‘卷积y(n)’);

   程序2：

           clf;

           h=[3  2  1  -2   1  0  -4  0  3];%冲激

           x=[1  -2  3  -4  3  2  1 ]; %输入序列

           y=conv(h,x);

           n=0:14;

           stem(n,y);

           xlabel(‘时间序号n’);ylabel(‘振幅’)；

           title(‘用卷积函数conv得到的输出’)；grid;

四、实验仪器设备

计算机，MATLAB软件

五、实验结论和心得

计算该卷积的理论结果，与两个程序的输出一致，若改变L的长度（改为8时），两个程序的输出一致，这种在卷积的时候我们利用时域的卷积等于频域的相乘的结果去进行基于FFT的卷积操作,这是一种数学技巧，起利用了公式转换，使得乘法次数减少，从而加快了算法的运行速度,对于这种方法我们可以推广到二维卷积