实验二用FFT作谱分析

1、实验二用FFT作谱分析时间：2021.03.09创作：欧阳法1. 实验目的(1)进一步加深DFT算法原理和基本性质的理解 (因为FFT只是DFT的一种快速算法，所以FFT的 运算结果必然满足DFT的基本性质)。熟悉FFT算法原理和FFT子程序的应用。(3)学习用FFT对连续信号和时域离散信号进行 谱分析的方法，了解可能出现的分析误差及其原因， 以便在实际中正确应用FFTo2. 实验步骤(1)复习DFT的定义、性质和用DFT作谱分析 的有关内容。复习FFT算法原理与编程思想，并对照DIT- FFT运算流图和程序框图，读懂本实验提供的FFT 子程序。(3)编制信号产生子程序，产生以下典型信号供 谱

2、分析用：编写主程序。/*DIT-FFT 函数(C 语言)*/fft基 2DITFFT 函数要求：指向复数数组指针X, FFT长度为2", m为正整数FFT输出结果放在输入复数数组中。/*计算N点FFT子程序*/* xr：=信号序列实部,xi：二信号序列虚部，N:=FFT变换区间长度N=2M*/*如果信号长度小于N,应该给xr, xi后面补()*/*计算如果X (K)的实部和虚部分别储存在数组xr和xi中*/Void Fft(doublc xr, double xi, int N, int M)int L, B,J,P, k, i;double rPartKB, iPartKB;dou

3、ble rCf128,iCf128/*计算旋转因子*/double PI2= 8.0 * atan(l.O);2021.03.09for(i=0; i<N; i+)rCfi=cos(i*PI2/N); iCfi=sin(i*PI2/N);Changc()rder( xr, xi, N );/*计算各级蝶形*/for(L=l;Lv=M;L+)B=(int)(pow(2, (L-l)+().5);for(J=();jv=B-l;J+)P=J*(int)(pow(2,(M-L)+0.5);for(k=J; k<=N-l; k+=(int)(pow(2,L)4-0.5)rPartKB =

4、xrk+B*rCfP-xik+B*iCfP;iPartKB = xik+B*rCfP+xrk+B*iCfPxrk+B = xrk- rPartKB;xik+B = xik- iPartKB;xrk = xrk + rPartKB;xik = xik + iPartKB;/*倒序子程序*/void Changc()rdor(doublc xrQ, double xi, int N )int LH,N1,I, J,K;double T;LH 二 N/2;J 二 LH; N1 二 N - 2;for(I=l;I<=Nl;I+)if(I<J)T = xrpj; xrl = xrJ; xrJ

5、 = T;K = LH;while(J>=I<)J 二 JKK = (int)(K/2+0.5);J=J+K;按实验内容要求，上机实验并写出实验报告。本实验采用的是MATLAB语言,因此FFT子程 序直接调用MATLAB语言中的FFT函数就可以实 现。下面给出完整的MATLAB程序%实验二，用FFT做谱分析b=mcnu(,请选择信号xl(n)x8(n)Txl(n)','x2(n)','x3(fi)','x4(n)','x5(ii)','x6(n)','x7=x4+x5：'x8=x

6、4+jx5：'Exif);if b=9b=();end=()；close all;while(b)if b=6tcmp=menu(,请选择FFT变换区间长度N7N=16,；N=32'；N=64,);if tcmp= = lN=16;elseif tcmp=2N=32;else N=64;endFs 二 64;n=():N-l;x=cos(8*pi*n/ fs)+cos(l 6*pi*n/ fs)+cos (20*pi*n / fs);elsetempmenuf请选择FFT变换区间长度N7N=8,；N=16,；N=32r);if tump二二 1N=8;elseif temp=2

7、N=16;else N 二32;endif b 二二 1x 二1 111() 0 0 ();else if b=2x 二1 2 3 4432 1;else iFb二二32021.03.09 x=4 3 2 1 1 2 34;else if b=4n=():N-l;x=cos(0.25*pi*n);else if b=5n=():N-l;x=sin(pi*n)/8);else if b=7n=():N-l;x=cos(n*pi/4)+sin(n*pi/8);else if b=8n=():N-l;x=cos(n*pi/4)+j*sin(n*pi/8);endendendend2021.03.09

8、endendendend%T() Calculate FFTf=fft(x,N);i=i+l;figuw(i);printf(x,abs (f) ,abs (N) ,abs (b);if N=16if b=7k=conj(f);X4= (f+k)/2;%ReX7(k) =x4(k)figurc(i+2);subplot(2,2,l);stcm(abs(x4)/.!);xlabelfk*);ylabtlC|X4(k)|*);tkl«恢复后的X4(k)*);x5=(f-k)/2;%jImX7(k)=X5(k) subplot(2,2,3);Stcm(abs(x5),.r);xlabelf

9、k*);ylabClC|X5(k)|*);tklW恢复后的X5(k)*);endif b=8k(l)=conj(f(l);for m=2:Nk(m)=conj(f(N-m+2);endfe=(x+k)/2;%求X8(k)的共辄对称分量fo=(x-k)/2;%求X8(k)的共毓反对称分量 xr=ifft(fe,N);%xr=x4(n)b 二 4;figure(i+l)printf(xr>abs(fc)>abs(N),abs(b);xi=ifft(fo,N)/j;%xi=x5(n)b 二 5;figure。+2)printf(xi,abs(Q,abs(N),abs(b);endendb

10、=mcnu(,请选择信号xl (n)x8(n)','xl(ii)','x2(n)Tx3(n)Tx4(n)','x5(n)','x6(n)Tx7=x4+x5：' x8=x4+jx5','Exif);if b=9b=();endclose all;end(5)按实验内容要求，上机实验，并写出实验报告。3. 上机实验内容(1)对2中所给出的信号逐个进行谱分析。下面 给出针对各信号的FFT变换区间N以及对连续信号 竝(z)的采样频率耳，供实验时参考。(n), x2(n), xii), *(坊：A-8, 16死(z)

11、:代二64Hz, A-16, 32, 64XiS)N=8N 二 16由以上两个图分析如下：离散傅里叶变换的N点变 换在频域范围内表现为对傅里叶变换即Z变换在单 位圆上的抽样。所以N取8点时,k=0,l,2,3,4,5,6,7与N取16点时,k=0,2,4,6,8,10,12,14的离散傅里叶变换 值对应相等，即它们都等于原信号在w=()、”/8、兀/4、3兀/8、4兀/8、5兀/8、6龙/8、6/8的傅里叶变 换，这在上面两图可以明显看出。所以，离散傅里叶 变换实际上是对该序列在频域范围内以2"/N的间隔 进行抽样。x2(n)N=8N=16xn) N=8N=16X4(n)N=8N=1

12、6从以上两组图可以看出，若按Xl(n)进行分析，则 明显不对。现分析如下：原信号周期为16,所以当 N=8时，未能取完一个周期的值，N=16则取完了一 个周期的值，所以这是两个不同的序列，所以按照 X1 (n)的分析方式是不对的，因为本身它们的傅里叶 变换就是不一样的。由于离散傅里叶变换是该序列周 期延拓后所对应的傅里叶级数变换的主值序列，所 以，当N=16时，所得的DFT值与X5(n)的傅里叶级 数变换的主值序列是一致的，而N二8时是X5(n)的部 分序列的周期延拓后的傅里叶级数变换的主值序列， 因此两者的值是不同的。X6(n)N=16N 二 32N=646*5 3 乂如 只曲也丁切 虫D致

13、 2 乃爲 0 口1/1 XIJ开苗I -U O O * I 制帀仗 | _3王1取町I 歩MMIQI3 t如C:Q?13原连续信号的周期为0.5,当采样频率二64Hz时, 所形成的序列周期为0.5*64=32。所以只有N>32,才 能取完一个周期的序列。这一点，从上面三个图可以 清晰看出。其中N二32和N二16的图形分析，可以参 考旳）的分析。N=32和N二64的图形分析，可以参 考花（劝的分析。 令双劝=无（劝+冯（力），用FFT计算8点和16点 离散傅里叶变换。二 DFT双功X7N二 8列的一部分，但是它们确们的傅里叶变换不同，即即DFT变换值不同。重复(2)。的频谱分析，从以上的图变换分为共辄对称部分合共个反对称部分，则可以得出原信号的实部对应离 散傅里叶变换的共馳对称部分，原信号的虚部对应离 散信号的共扼反对称部分。4. 思考题(1)在 g 时，花(坊和码(功的幅频特性会相同 吗？为什么？ A16呢？N= 8时一样，N= 1 6时不一样。因为DFT变换可以看成是将该序列进行周期延拓后的傅 里叶级数变换的主值序列。当N二8时，两序列进行 周期延拓后序列相同，所以其傅里叶级数变换的主值 序列也相同，进而DFT变换也相同。而当N二16 时，两序列进行周期延拓后序列不相同，所以其傅里 叶级数变换的主值