实验四离散傅里叶变换及其快速算法.doc

实验四 离散傅里叶变换及其快速算法一、 实验目的掌握快速傅立叶变换的应用方法； 二、 实验仪器：电脑一台，MATLAB6.5或更高级版本软件一套。三、 实验原理和实例分析 （一）离散傅里叶变换离散傅立叶级数变换是周期序列，仍不便于计算机计算。但离散傅立叶级数虽是周期序列，却只有N个独立的数值，所以它的许多特性可以通过有限长序列延拓来得到。对于一个长度为N的有限长序列，也即只在个点上有非零值，其余皆为零，即把序列以N为周期进行周期延拓得到周期序列，则有：所以，有限长序列的离散傅立叶变换(DFT)为：逆变换为：若将DFT变换的定义写成矩阵形式，则得到：X=Ax，其中DFT变换矩阵A为Dftmtx 函数：用来计算DFT变换矩阵A的函数调用方式Adftmta（n）：返回nn的DFT变换矩阵A。若x为给定长度的行向量，则yx*A，返回x的DFT变换y。Aiconj（dftmtx（n）/n；返回nn的IDFT变换矩阵Ai。【实例4-1】 A=dftmtx(4) Ai=conj(dftmtx(4)/4运行结果A = 1.0000 1.0000 1.0000 1.0000 1.0000 0 - 1.0000i -1.0000 0 + 1.0000i 1.0000 -1.0000 1.0000 -1.0000 1.0000 0 + 1.0000i -1.0000 0 - 1.0000iAi = 0.2500 0.2500 0.2500 0.2500 0.2500 0 + 0.2500i -0.2500 0 - 0.2500i 0.2500 -0.2500 0.2500 -0.2500 0.2500 0 - 0.2500i -0.2500 0 + 0.2500i【实例4-2】如果是一个N16的有限序列，用MATLAB求其DFT的结果，并画出其结果图，如图41所示。% 程序N=16;n=0:1:N-1; %时域采样xn=sin(n*pi/8)+sin(n*pi/4);k=0:1:N-1; %频域采样WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=xn*WNnk;subplot(2,1,1)stem(n,xn);subplot(2,1,2)stem(k,abs(Xk);运算结果Xk = Columns 1 through 5 0.0000 -0.0000 - 8.0000i -0.0000 - 8.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i Columns 6 through 10 -0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i Columns 11 through 15 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 8.0000i Column 16 0.0000 + 8.0000i图 4-1 有限长序列的DFT结果图思考：当N=8，16，24，48，分别执行程序，得出什么结论？频谱的误差主要是由于时域中对信号的非整周期截断产生的频谱泄漏。 （二） 快速傅立叶变换（FFT）快速离散傅里叶变换是计算离散傅里叶变换的一种快速算法，为了提高运算速度，FFT将DFT的计算逐次分解成较小点数的DFT。按时间抽取Decimation-In-Time(DIT)FFT算法把输入序列按其值为偶数或是奇数分解成越来越短的序列。按频域抽取（Decimation-In-Frequency(DIF)FFT算法是把输出序列按其值是偶数或是奇数来分解成越来越短的序列。本实验重点是针对的DIT FFT或DIF FFT。MATLAB为计算数据的离散快速傅立叶变换，提供了一系列丰富的数学函数，主要有fft、ifft、fft2 、ifft2, fftn、ifftn和fftshift、ifftshift等。当所处理的数据的长度为2的幂次时，采用基-2算法进行计算，计算速度会显著增加。所以，要尽可能使所要处理的数据长度为2的幂次或者用添零的方式来添补数据使之成为2的幂次。函数fft用来求序列的DFT,调用格式为：Xk=fft(xn,N)其中，xn为有限长序列，N为序列xn的长度，Xk为序列xn的DFT.函数ifft用来求IDFT，调用格式为：xn=ifft(Xk,N)其中，Xk为有限长序列，N为序列Xk的长度，xn为序列Xk的IDFT.1、快速傅里叶变换对信号进行频谱分析【实例4-3】已知序列，试绘制及它的离散傅里叶变换图。MATLAB实现程序：clear allN=100;n=0:N-1;xn=2*sin(0.48*pi*n)+cos(0.52*pi*n);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1)plot(n,xn)xlabel(n);ylabel(x(n);title(x(n) N=100);subplot(1,2,2)k=0:length(magXK)-1;stem(k,magXK,.);xlabel(k);ylabel(|X(k)|);title(X(k) N=100);运行结果：图4-2 序列x(n)及利用FFT求解的相应频谱思考：若N=80或N=120频谱有什么变换？2、用FFT实现两序列的线性卷积【实例4-5】有限长序列x(n)=2,1,0,1,3;h(n)=1,3,2,1，试利用FFT实现DFT计算线性卷积，并与线性卷积直接计算（conv）的结果进行比较。% 用FFT实现由DFT计算线性卷积x=1 2 0 1 3;h=1 3 2 1;L=length(x)+length(h)-1; XE=fft(x,L);HE=fft(h,L); y1=ifft(XE.*HE);%画出由圆周卷积计算线性卷积结果及误差k=0:L-1;subplot(211);stem(k,real(y1);axis(0 7 0 20);title(Result of Linear Convolution);xlabel(Time index k);ylabel(Amplitude);y2=conv(x,h);error=y1-y2;subplot(212);stem(k,abs(error);xlabel(Time index k);ylabel(Amplitude); title(Error Magnitude);若h=1 3 2 1 1; 如何修改程序？四、 实验内容1、对连续的单一频率周期信号, 按采样频率 采样，截取长度N分别选N =20和N =16，观察其DFT结果的幅度谱，并分析它们结果误差的原因。2、试用Mablab求其有限长序列与的圆周卷积，（N=20），并画出其结果图。五、 思考题1、利用DFT对连续信号进行傅里叶分析可能造成哪些误差？2、FFT的变换长度L必须满足什么条件，输出y(n)才等于x(n)和h(n)的线性卷积，这种快速卷积的方法是利用什么运算与线性卷积运算的关系得到的？x=1 2 0 1 3;h=1 3 2 1 1;%L=length(x)+length(h)-1; XE=fft(x,16);HE=fft(h,16); y1=ifft(XE.*HE);%k=length(y1)kn=0:k-1;subplot(211);stem(kn,real(y1);%axis(0 7 0 20);title(Result of Linear Convolution);xlabel(Time index k);yl