离散傅里叶变换和快速傅里叶变换

1、戶幵,戈丿、弟实验报告课程名称：彳指导老师 成绩：实验名称：离散傅里叶变换和快速傅里叶变换实验类型: 同组学生姓名：一、实验目的和要求(必填)二、实验内容和原理(必填)三、主要仪器设备(必填)四、操作方法和实验步骤五、实验数据记录和处理六、实验结果与分析(必填)七、讨论、心得一、实验目的和要求1. 掌握DFT的原理和实现2. 掌握FFT的原理和实现，掌握用FFT对连续信号和离散信号进行谱分析的方法。二、实验内容和原理2.1 DTFT 和 DFT序列x(n)的离散事件傅里叶变换(DTFT )表示为:X(ej )x( n)eN 1如果x(n)为因果有限长序列,n=0,1,.,N-1，则x(n)的D

2、TFT表示为：X(ej ) x(n)en 0x(n)的离散傅里叶变换(DFT)表达式为:X(k)x(n)en 0jnkN(k0,1,N 1)序列的N点DFT是DTFT在0,2 n上的N点等间隔采样，采样间隔为2 dN。通过DFT , 可以完成由一组有限个信号采样值x(n)直接计算得到一组有限个频谱采样值X(k)。X(k)的幅度谱为X(k) vxR(k) X|2(k) , XR(k)和Xi(k)分别为X(k)的实部和虚部。X(k)的相位谱为(k)列吩离散傅里叶反变换IDFT)定义为 x(n)丄 N (k)ej_Nnk (n 0,1，N 1)N n 02.2 FFT快速傅里叶变换(FFT)是DFT

3、的快速算法，它减少了DFT的运算量，使数字信号的处理速度大大提高。三、主要仪器设备PC 一台，matlab 软件四、实验内容4.1第一题求有限长离散时间信号x(n)的离散时间 傅里叶变换(DTFT ) X(ejQ)并绘图。(1)已知 x(n)0 其2他n 2; (2)已知 x(n) 2n 0 n 10。0其他4.1.1理论分析1) 由DTFT计算式,X( ) x(n)ejne2j15je1 e je2e20.5je0.5 jesin(2.5 )sin(0.5 )X (Q)是实数，可以直接作出它的图像。Figure 1 X (Q)曲线2) 由DTFT计算式:X(x(n)e jio2ne jn 0

4、1夕1111 2e jz11。111 2e jX(2111 2e j可以发现(Q)周期为2 n;而 X (Q)的相位在2n周期内有约十次振荡。4.1.2编程计算作图编写一个计算 DTFT的函数。fun cti on DTFT(x, n1,n2) w=-2*pi:2*pi/1000:2*pi; % X=zeros(size(w);for i=n 1: n2 %DTFT表示Q计算式X=X+x(i-n1+1)*exp(-1)*j*w*i);endan gle(X);subplot(2,1,1);plot(w, abs(X),r);xlabel( Omega );ylabel( subplot(2,1

5、,2);plot(w,a ngle(X),b);xlabel( Omega );ylabel( end|X(Omega)| );hold on; % 作幅频图an gle(Omega); % 作相频图输入序列x,和n的取值范围，即可计算其 DTFT。1) 输入：x=1 1 1 1 1;DTFT(x,-2,2);(因为X (Q)是实数，所以实际计算过程中对相频曲线取了绝对值) 结果：Figure 2 X (Q)的频谱可以看出，X (Q)的相位只有 0和n两种取值，X (Q)是实函数，而且其幅度频谱 与理论计算得到的相同。2) 输入：n=0:10; x=2.A n;DTFT(x,0,10)结果:2

6、500200015001000-10-5050Figure 3第1题(2)中X (Q)的频谱500-14.2第二题已知有限长序列x(n)=色7,9,5,1,7,9,5,试分别采用DFT和FFT求其离散傅里叶变换 X(k)的幅 度、相位图。4.2.1理论分析由 FFT 蝶形运算得到，X(k)= 51,7,-9-j4,7,3,7,-9+j4,74.2.2编程计算作图1. DFT编写一个计算 DTFT的函数。DFT (序列x，长度N )fun ctio nDFT(x,N)k=0:N-1;X=zeros(size(k);for n=0:N-1X=X+x( n+1)*exp(-1)*j*2*pi/N*

7、n*k);%DFT计算式endsubplot(2,1,1);stem(k,abs(X), .); xlabel( k );ylabel( |X(k)| );hold on; %幅频图subplot(2,1,2);stem(k,angle(X),* );xlabel( k );ylabel(Angle(k); %相频图end输入：x=8 7 9 5 1 7 9 5;DFT(x,8);结果：Figure 4 第2题DFT结果2. FFT编写一个利用matlab自带函数计算FFT并绘图的函数FFT1 (序列x，长度N)fun cti onFFT1(x,N)X=fft(x,8);%用自带的fft函数计

8、算k=0:N-1;subplot(2,1,1);stem(k,abs(X), .); xlabel( k );ylabel( |X(k)| );hold on; %幅频图subplot(2,1,2);stem(k,angle(X),* );xlabel( k );ylabel(Angle(k); %相目频图end输入：x=8 7 9 5 1 7 9 5;FFT1(x,8);结果：Figure 5 第2题FFT结果因为FFT只是DFT的一种快速算法，所以 FFT的结果与DFT结果相同。DFT和FFT的结果，符合理论计算得到的，X(k)= 51,7,-9-j4,7,3,7,-9+j4,74.3第三

9、题已知连续时间信号 x(t)=3cos8 n t, X( 3 )= 3 (8 )(8 )，该信号从t=0开始以采样周期Ts=0.1 s进行采样得到序列 x(n)，试选择合适的采样点数，分别采用DFT和FFT 求其离散傅里叶变换X(k)的幅度、相位图，并将结果与X(k)的幅度、相位图，并将结果与X( 3 )相比较。4.3.1理论分析1.原信号的频谱：X( 3 )=3 (8 )(8 ),只在土 8n不为0.且在土 8n处相2. 采样角频率Q S=20 n 2 X 8 n,满足采样定理。3. 采样后的信号，为 X( 3 )以20 n为周期的延拓。所以只在(土 8+20k )n( k为任意整数)处不为

10、0.如取区间0,20 n 内，只有8n和16 n处不为0。进行N点DFT后，将81220 n的区间映射为0,N区间。理想情况下仅在n=N和N两处不为0。20204. x(n) cos(8 nTs) cos(0.8 n)，周期为5,所以取采样点数为 5的倍数时，不会发生泄漏；而采样点数不是 5的倍数时，则会发生泄漏。101331111319-876)X 54310401Figure 6原始信号的频谱 X (Q)4.3.2编程计算作图编写一个获得信号的N点样本的函数 sample (点数N)fun cti on x=sample(N)t=0:0.1:(N-1)*0.1; %0.1s为间隔x=3*c

11、os(8*pi*t); %x即采样结果End输入：X=sample(N);FFT1(X,N);即可获得采样 N点的频谱图。因为 DFT结果与FFT是完全一样的，所以这里只使用FFT作图。取采样点数N=5 16 20 104获得以下频谱图： -.o51J2.5133.5Figure 7N=5Figure 8N=16Figure 9 N=2080Figure 10 N=104可以看出，N=5和20时，由于是周期的整数倍，频谱只有两条谱线，且满足前面812理论计算得出的公式n= N和宜 N，没有发生泄漏，且这两条谱线对应2020的相位是0.所以频谱与原信号频谱在形式上时相同的。而N=16和N=54时

12、，则都发生了频谱泄漏，频谱与原信号频谱就很不同了。但相比之下N=54时谱线更加接近原谱线。验证了 “为减小泄漏误差，如果待分析的信号实现不知道确切周期，则截取较长时间长度的样点进行分析”这个说法。同时也可以发现，虽然幅频图中显示幅值为0，但相频图中相应的位置仍有谱线。这可能是matlab浮点运算造成的误差，即本来为0处其实是一个非常小的复数， 所以仍有一定相位。4.4第四题4.4.1理论分析若噪声信号较小，则采样后的频谱仍能较准确地反映原信号的特征。4.4.2编程计算对原采样程序稍加改编，加入一个噪声信号p*randn(1,N)。p表示噪声信号的强度。fun cti onx=sample no

13、 ise(N,p)t=0:0.1:(N-1)*0.1;x=3*cos(8*pi*t)+p*ra ndn (1,N);end取采样点数N=20进行分析。输入：X=samplenoise(20,p); % 取P=1 和 10 两种情况。FFT1(X,20);30Figure 11N=20 噪声较小(p=1 )-02468101214161820k-2A4210-220e g AFigure 12N=20 噪声较大(p=10 )-4 IIIIIIIl_024681012141618k可见，较小的噪声对信号的频谱的影响不大,仍能较精确地获得频谱图。而噪声较大时则很难准确获得原信号的频谱。4.5第五题3

14、.5 已知序列 x(n) 4 (n)3 (n 1) 2 (n 2)(n 3)，X(k)是 x(n)的 6 点 DFT ,设(1)若有限长序列y(n)的6点DFT是Y(k) W64kX(k)，求y(n)。(2) 若有限长序列 w(n)的6点DFT W(k)是X (k)的实部，求w(n)。(3) 若有限长序列 q(n)的 3 点 DFT 是 Q(k) X(2k), k=0,1,2，求 q(n)。由题意得到：x( n)=4,3, 2, 1,0, 0 n=(0,1,2,3,4,5)4.5.1理论分析1)由DFT的性质可以得到，如果y (n)的DFT为Y(k) W64kX(k)，那么y(n)就是x(n)

15、圆周左移4位得到的。所以y (n) =0,0,4,3,2,12)3.5+j4.3301由题意，X(k)=103.5-j4.33012.5-8.66j22.5+j0.866取实部，则 W(k)=103.52.5 2.52 3.5按照IDFT计算式x(n)-N 1X(k)ejnk(n 0,1,., N 1)计算得到N n 0w(n )=4 1.5 1 1 1.5由DFT的性质，因为 W(k)是实数，所以对应的 w(n)也是实数。3) 由题意，Q(k)= 102.5-8.66j2.5+j0.866按照IDFT计算式x(n)1X(k)ePk0(nq(n )=5324.5.2编程计算1)x=4 3 2

16、1 0 0;X=fft(x,6);%求 DFTk=0:5;ifft(exp(j*4*pi/3*k).*X) %求 IFFT结果：arts =Calumns 1 Throueh 5D. QOOO 呻 O.ODOOi ChOODQ + OcQOOOi 4.C0Q0 - 0. OOOOi 3. (WOO - 0. OOQOi 2.QD0D - GQQMHColunn 0LOOOO + 0.00001y (n) =0,0,4,3,2,1 即x(n)圆周左移4位，与理论值相同。2)%接第(1)题的程序W=real(X);% 取实部ifft(W)结果：ans =4.0000 l-oOOO 1.0000L0

17、0001,00001-5000w( n)=41.5 1 1 1.5与理论计算值相同。3)%接第(1)题的程序Q=X(1) X(3) X(5); %Q 是 X( 2k),由于 matlab 的矩阵是从 X( 1)开始, ifft(Q)%所以对应的应该是第1、3、5个元素结果：305 =32q( n)=53 2 ，与理论值相同。已知信号x(t)4.6第六题sin(2f1t)sin(2f2t)sin(2f3t)，其中f1=4Hz、f2=4.02 Hz、f3=5 Hz，采用采样频率为20 Hz进行采样，求：(1) 当采样长度 N分别为512和2048情况下x(t)的幅度频谱；(2) 当采样长度 N为3

18、2,且增补N个零点、4N个零点、8N个零点、16N个零点情况下x(t)的幅度频谱。461理论分析1. 首先20Hz的采样频率是满足采样定理的。2. 频率分辨率是DFT中谱线间的最小间隔，单位是Hz。对于长度为N的序列，频率分辨率为fs/N,为采样频率。3. 因为采样点数N不是周期的整数倍，所以一定会存在频谱泄露情况。4. fs=20Hz时,N=512,则分辨率 20/512疋0.039Hz 0.02Hz所以不能区分开信号中频率为 4Hz禾口 4.02Hz的两个分量。5. N=2048,则分辨率 20/2048 0.01Hz0.02Hz 可以区分开 4Hz和4.02Hz的这两个分量。4.6.2编

19、程作图编写一个取 N个点并补充t*N个0的函数sample2(取样点数N,补零t)fun cti onx=sample2(N,t)n=0:0.02:0.02*(N-1); %N个点x=si n(2*pi*4* n)+si n(2*pi*4.02* n)+si n(2*pi*5* n); %取样过程x(N+1:N*(t+1)=0; %补零end1) 输入：x=sample2(N,0); %N=512,4096FFT1(x,N);结果：0100200300400500600k300250200150100500Figure 13N=512 幅频图N=512时，可以看到对 4Hz和5Hz的分量是明显区分开来的。对峰值附近放大来看:300Figure 14N=512幅频图放大看到k=103和104处有峰值。44 02理论计算中，512 102.4512 102.920 20由于频谱泄露，在 k=103和104处出现峰值，实际上并没有将4Hz和4.02Hz的两个分量区分开来。kFigure 15 N=2048 幅频图N=2048时，显然区分开了 4H