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

实验报告课程名称：信号分析与处理 指导老师 成绩： 实验名称：离散傅里叶变换和快速傅里叶变换 实验类型： 同组学生姓名： 一、实验目的和要求（必填）二、实验内容和原理（必填）三、主要仪器设备（必填）四、操作方法和实验步骤五、实验数据记录和处理六、实验结果与分析（必填）七、讨论、心得一、实验目的和要求1. 掌握DFT的原理和实现2. 掌握FFT的原理和实现，掌握用FFT对连续信号和离散信号进行谱分析的方法。二、实验内容和原理2.1 DTFT和DFT序列x(n)的离散事件傅里叶变换（DTFT）表示为：，如果x(n)为因果有限长序列，n=0,1,.,N-1，则x(n)的DTFT表示为：，x(n)的离散傅里叶变换（DFT）表达式为：，序列的N点DFT是DTFT在 0,2上的N点等间隔采样，采样间隔为2/N。通过DFT，可以完成由一组有限个信号采样值x(n)直接计算得到一组有限个频谱采样值X(k)。X(k)的幅度谱为，XR(k)和XI(k)分别为X(k)的实部和虚部。X(k)的相位谱为。离散傅里叶反变换（IDFT）定义为。2.2 FFT快速傅里叶变换（FFT）是DFT的快速算法，它减少了DFT的运算量，使数字信号的处理速度大大提高。三、主要仪器设备PC一台，matlab软件四、实验内容4.1第一题求有限长离散时间信号x(n)的离散时间傅里叶变换（DTFT）X(ej)并绘图。（1）已知；（2）已知。4.1.1理论分析1) 由DTFT计算式， X（）是实数，可以直接作出它的图像。Figure 1 X（）曲线2) 由DTFT计算式：可以发现X（）周期为2；而X（）的相位在2周期内有约十次振荡。4.1.2编程计算作图编写一个计算DTFT的函数。function DTFT(x,n1,n2)w=-2*pi:2*pi/1000:2*pi; %表示X=zeros(size(w); for i=n1:n2 %DTFT计算式 X=X+x(i-n1+1)*exp(-1)*j*w*i);endangle(X);subplot(2,1,1);plot(w, abs(X),r);xlabel(Omega);ylabel(|X(Omega)|);hold on; %作幅频图subplot(2,1,2);plot(w,angle(X),b);xlabel(Omega);ylabel(angle(Omega); %作相频图end输入序列x，和n的取值范围，即可计算其DTFT。1) 输入：x=1 1 1 1 1;DTFT(x,-2,2);（因为X（）是实数，所以实际计算过程中对相频曲线取了绝对值）结果：Figure 2 X（）的频谱可以看出，X（）的相位只有0和两种取值，X（）是实函数，而且其幅度频谱与理论计算得到的相同。2) 输入：n=0:10;x=2.n;DTFT(x,0,10)结果：Figure 3 第1题（2）中X（）的频谱4.2第二题已知有限长序列x(n)=8,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）function DFT(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*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）function FFT1(x,N)X=fft(x,8); %用自带的fft函数计算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,7 。4.3第三题已知连续时间信号x(t)=3cos8t, X()= ，该信号从t=0开始以采样周期Ts=0.1 s进行采样得到序列x(n)，试选择合适的采样点数，分别采用DFT和 FFT求其离散傅里叶变换X(k)的幅度、相位图，并将结果与X(k)的幅度、相位图，并将结果与X()相比较。4.3.1理论分析1. 原信号的频谱：X()=，只在8不为0. 且在8处相位为0 。2. 采样角频率s=2028，满足采样定理。3. 采样后的信号，为X()以20为周期的延拓。所以只在（8+20k） （k为任意整数）处不为0. 如取区间0,20内，只有8和16处不为0 。进行N点DFT后，将20的区间映射为0,N区间。理想情况下仅在 n= 和两处不为0 。4. ，周期为5，所以取采样点数为5的倍数时，不会发生泄漏；而采样点数不是5的倍数时，则会发生泄漏。Figure 6 原始信号的频谱X（）4.3.2编程计算作图编写一个获得信号的N点样本的函数 sample（点数N）function x=sample(N)t=0:0.1:(N-1)*0.1; %0.1s为间隔x=3*cos(8*pi*t); %x即采样结果。End输入:X=sample(N);FFT1(X,N);即可获得采样N点的频谱图。因为DFT结果与FFT是完全一样的，所以这里只使用FFT作图。取采样点数N=5 16 20 104 获得以下频谱图：Figure 7 N=5Figure 8 N=16Figure 9 N=20Figure 10 N=104可以看出，N=5和20时，由于是周期的整数倍，频谱只有两条谱线，且满足前面理论计算得出的公式n= 和，没有发生泄漏，且这两条谱线对应的相位是0.所以频谱与原信号频谱在形式上时相同的。而N=16和N=54时，则都发生了频谱泄漏，频谱与原信号频谱就很不同了。但相比之下N=54时谱线更加接近原谱线。验证了“为减小泄漏误差，如果待分析的信号实现不知道确切周期，则截取较长时间长度的样点进行分析”这个说法。同时也可以发现，虽然幅频图中显示幅值为0，但相频图中相应的位置仍有谱线。这可能是matlab浮点运算造成的误差，即本来为0处其实是一个非常小的复数，所以仍有一定相位。4.4第四题4.4.1理论分析若噪声信号较小，则采样后的频谱仍能较准确地反映原信号的特征。4.4.2编程计算对原采样程序稍加改编，加入一个噪声信号p*randn(1,N)。p表示噪声信号的强度。function x=samplenoise(N,p)t=0:0.1:(N-1)*0.1;x=3*cos(8*pi*t)+p*randn(1,N);end取采样点数N=20进行分析。输入：X=samplenoise(20,p); %取P=1和10两种情况。FFT1(X,20);Figure 11 N=20 噪声较小（p=1）Figure 12 N=20 噪声较大（p=10）可见，较小的噪声对信号的频谱的影响不大，仍能较精确地获得频谱图。而噪声较大时则很难准确获得原信号的频谱。4.5第五题3.5已知序列，X(k)是x(n)的6点DFT，设。（1） 若有限长序列y(n)的6点DFT是，求y(n)。（2） 若有限长序列w(n)的6点DFT W(k)是的实部，求w(n)。（3） 若有限长序列q(n)的3点DFT是，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(n)就是x(n)圆周左移4位得到的。所以y（n）=0,0,4,3,2,12) 由题意，X(k)=10 3.5-j4.3301 2.5-8.66j 2 2.5+j0.866 3.5+j4.3301取实部，则W(k)=10 3.5 2.5 2.5 2 3.5 按照IDFT计算式计算得到w(n)=4 1.5 1 1 1.5 由DFT的性质，因为W(k)是实数，所以对应的w(n)也是实数。3) 由题意，Q(k)= 10 2.5-8.66j 2.5+j0.866 按照IDFT计算式计算得到q(n)=5 3 24.5.2编程计算1)x=4 3 2 1 0 0;X=fft(x,6); %求DFTk=0:5;ifft(exp(j*4*pi/3*k).*X) %求IFFT 结果：y（n）=0,0,4,3,2,1 即x(n)圆周左移4位，与理论值相同。2)%接第（1）题的程序W=real(X); %取实部ifft(W)结果：w(n)=4 1.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个元素结果：q(n)=5 3 2 ，与理论值相同。4.6第六题已知信号，其中f1=4 Hz、f2=4.02 Hz、f3=5 Hz，采用采样频率为20 Hz进行采样，求:（1） 当采样长度N分别为512和2048情况下x(t)的幅度频谱；（2） 当采样长度N为32，且增补N个零点、4N个零点、8N个零点、16N个零点情况下x(t)的幅度频谱。4.6.1理论分析1. 首先20Hz的采样频率是满足采样定理的。2. 频率分辨率是DFT中谱线间的最小间隔，单位是Hz。对于长度为N的序列，频率分辨率为fs/N, 为采样频率。3. 因为采样点数N不是周期的整数倍，所以一定会存在频谱泄露情况。4. fs=20Hz时,N=512,则分辨率 20/5120.039Hz 0.02Hz 所以不能区分开信号中频率为4Hz和4.02Hz的两个分量。5. N=2048,则分辨率20/20480.01Hz0.02Hz 可以区分开4Hz和4.02Hz的这两个分量。4.6.2编程作图编写一个取N个点并补充t*N个0的函数sample2(取样点数N,补零t)function x=sample2(N,t)n=0:0.02:0.02*(N-1); %N个点x=sin(2*pi*4*n)+sin(2*pi*4.02*n)+sin(2*pi*5*n); %取样过程x(N+1:N*(t+1)=0; %补零end1) 输入：x=sample2(N,0); %N=512 , 4096FFT1(x,N);结果：Figure 13 N=512幅频图N=512时，可以看到对4Hz和5Hz的分量是明显区分开来的。对峰值附近放大来看：Figure 14 N=512 幅频图放大看到k=103和104处有峰值。理论计算中， 由于频谱泄露，在k=103和104处出现峰值，实际上并没有将4Hz和4.02Hz的两个分量区分开来。Figure 15 N=2048 幅频图N=2048时，显然区分开了4Hz和5Hz的分量。峰值附近放大看：在k=409、410和k=412有峰值。理论上 由于频谱泄露，在409.6附近的409和410处同时出现峰值，在411.6附近的412出现峰值。所以采样点数N=2048时一定程度上区分开了4Hz和4.02Hz的分量。2) 输入x=sample2(32,t); %t=取1,4,8,16FFT1(x,N); %N=32*(t+1)结果：Figure 16 补32个0 Figure 17 补128个0Figure 18 补256个0Figure 19 补512个0可以发现，采样点数相同，都是32。20/32=0.61Hz，频谱图可以区分4Hz和5Hz分量，但不可以区分4Hz和4.02Hz的分量。补不同个数的0，幅频图的总体形状基本上是不变的。