实验四 FFT的计算机实现及谱分析一、实验目的与要求(1) 了解DFT的性质及应用。(2) 熟悉MATLAB编程的特点。 (3)掌握有限离散傅里叶变换的计算机实现方法。()进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法，所以FPT的运算结果必然满足DFT的基本性质)。 ()熟悉FFT算法原理和FFT子程序的应用。 ()学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。二、实验原理与方法.的计算原理设，有限离散信号，则它的有限离散频谱为令则上式为把信号按下标偶数项和奇数项分成两部分，即令它们的有限离散频谱分别为经分析，我们可得如下的时域分解FFT算法若我们把按偶奇序号分成两部分，则可得如下的频域分解FFT算法、MATLAB提供的快速傅立叶变换函数（1）fftfft函数用于计算快速傅立叶变换，其语法格式为：B = fft(I)B = fft(I)返回信号I的fft变换矩阵，输入信号I和输出信号B大小相同。（2）fftshiftMATLAB提供的fftshift函数用于将变换后的信号频谱中心从矩阵的原点移到矩阵的中心，其语法格式为：B = fftshift(I)对于矩阵I，B = fftshift(I)将I的一、三象限和二、四象限进行互换。（3）ifftifft函数用于计算信号的傅立叶反变换，其语法格式为：B = ifft(I)B = ifft(A)返回信号I的傅立叶反变换矩阵，输入信号I和输出信号B大小相同。其语法格式含义与fft函数的语法格式相同，可以参考fft函数的说明。时域分解FFT算法的几点说明：1、对一般的，如何求信号的一项频谱排列？我们可通过二进制的逆序得到：假设的二进制为，则2、碟形运算两节点的距离：第一级（第一列）两节点间的距离为1；第二级（第二列）两节点间的距离为2；.第m级（第m列）两节点间的距离为；3、的确定因子最后一列有种，顺序为为；因子倒数第二列有种，为 ；因子倒数第三列有种，为；因子第一列有一种，为1频域分解FFT算法的几点说明：1、对一般的，最后得到的逆序频谱如何排列才能得到原信号的频谱？我们可通过二进制的逆序得到：假设的二进制为，则2、碟形运算两节点的距离：第一级（第一列，m1）两节点间的距离为；.第m级（第m列）两节点间的距离为；3、的确定因子第一列有种，顺序为为；因子第二列有种，为 ；因子第三列有种，为；因子最后一列有一种，为1三 实验内容及步骤 1、用三种不同的DFT程序计算的傅里叶变换,并比较三种程序计算机运行时间。（1）用for 语句的M函数文件dft1.m，用循环变量逐点计算X(k)；（2）编写时域分解FFT 算法或频域分解FFT 算法的M函数文件dft2.m；（3）调用FFT库函数，直接计算X(k)；（4）分别利用上述三种不同方式编写的DFT程序计算序列x(n)的傅立叶变换，并画出相应的幅频和相频特性，再比较各个程序的计算机运行时间。()复习DFT的定义、性质和用DFT作谱分析的有关内容。 ()复习FFT算法原理与编程思想，并对照DITFFT运算流图和程序框图，读懂本实验提供的FFT子程序。()编制信号产生子程序，产生以下典型信号供谱分析用： ()编写主程序。 ()按实验内容要求，上机实验，并写出实验报告。M函数文件如下： dft1.m:functionAm,pha=dft1(x)N=length(x);w=exp(-j*2*pi/N);for k=1:N sum=0; for n=1:N sum=sum+x(n)*w(k-1)*(n-1); end Am(k)=abs(sum); pha(k)=angle(sum);enddft3.m:functionAm,pha=dft3(x)Xk=fft(x);Am=abs(Xk);pha=angle(Xk); 源程序及运行结果：（1）N=256;x=ones(1,8),zeros(1,N-8);t=cputime;Am1,pha1=dft1(x);t1=cputime-tn=0:(length(x)-1);w=(2*pi/length(x)*n;figure(1)subplot(2,1,1), plot(w,Am1,b); grid;title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);subplot(2,1,2), plot(w,pha1,r); grid;title(Phase Part);xlabel(frequency in radians);ylabel(argXexp(jw)/radians); 运行时间：t1=0.288s（3）N=256;x=ones(1,8),zeros(1,N-8);t=cputime;Am3,pha3=dft3(x);t3=cputime-t;n=0:(length(x)-1);w=(2*pi/length(x)*n;figure(3)subplot(2,1,1), plot(w,Am3,b); grid;title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);subplot(2,1,2), plot(w,pha3,r); grid;title(Phase Part);xlabel(frequency in radians);ylabel(argXexp(jw)/radians);运行时间：t3=0从以上运行结果可以看出，调用FFT库函数直接计算X(k)速度最快，用循环变量逐点计算运行速度最慢。因此FFT算法大大提高了DFT的实用性。四、研究高密度频谱和高分辨率频谱设有连续信号Xa(t)=cos(26.5103t)+cos(27103t)+cos(29103t)以采样频率fs=32kHz对信号Xa(t)采样,分析下列三种情况的幅频特性。（1） 采集数据长度N=16点，做N=16点的DFT，并画出幅频特性。（2） 采集数据长度N=16点，补零到256点，做N=256点的DFT，并画出幅频特性。（3） 采集数据长度N=256点，做N=256点的DFT，并画出幅频特性。观察三幅不同频率特性图，分析和比较它们的特点以及形成的原因。源程序如下：fs=32000; %采样频率N=16; %采集16点n=0:N-1; %做16点DFTxa=cos(2*pi*6.5*103*n/fs)+cos(2*pi*7*103*n/fs)+cos(2*pi*9*103*n/fs);Am1,pha1=dft3(xa);n=0:(length(xa)-1);w=(2*pi/length(xa)*n;figure(1)plot(w,Am1,b); grid;title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);x=xa(1:16),zeros(1,240); %补零到256Am2,pha2=dft2(x); %做256点DFTn=0:(length(x)-1);w=(2*pi/length(x)*n;figure(2)plot(w,Am2,b); grid;title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);N0=256; %采集256点n0=0:N0-1; %做256点DFTxa0=cos(2*pi*6.5*103*n0/fs)+cos(2*pi*7*103*n0/fs)+cos(2*pi*9*103*n0/fs);Am3,pha3=dft3(xa0);n=0:(length(xa0)-1);w=(2*pi/length(xa0)*n;figure(3)plot(w,Am3,b); grid;title(Magnitude part);xlabel(frequency in radians); ylabel(|X(exp(jw)|); 运行结果：N=16点DFT幅频特性 高密度频谱： 高分辨率频谱：观察运行结果可知，N=16点时的DFT，由于点数太少，很难反映出频谱的细节特征，只能分辨出两个频率分量。将16点DFT补零到256点，做256点DFT,由结果可以看出，频率间隔缩小（由2p/16减小为2 p/256），得到了比较平滑的连续曲线，此为高密度频谱。它是由16点经过补零所得到的，并没有增加任何信息，频率分辨率并未提高，仍然只能看出两个频率分量。第三次做的是256点的DFT，采集点数大大增加，分辨率提高，可以清楚看到三根谱线的位置，是为高分辨率频谱。五、上机实验内容(1)对三中所给出的信号逐个进行谱分析。 (2)令，用FFT计算8点和16点离散傅里叶变换。 (3)令，重复（2）。(4)对第四部分给以实现并进行说明。六、思考题（1）在N=8时，和的幅频特性会相同吗？为什么？N=16呢？（2）如果周期信号的周期预先不知道，如何使用FFT进行谱分析？七、实验报告要求（1）简述实验原理及目的。（2）按实验步骤附上实验过程中的信号序列和幅频特性曲线，并对所得结果进行分析和解释。结合实验中所得给定典型序列幅频特性曲线，与理论结果比较，并分析说明误差产生的原因以及用FFT作谱分析时有关参数的选择方法。（3）总结实验所得主要结论。（4）简要回答思考题。