用FFT作谱分析(附matlab程序)

一、实验目的学习用FFT对连续信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便正确应用FFT。二、实验原理1.对连续时间信号进行谱分析，首先应该进行时域采样，采样频率Fs要满足乃奎斯特采样定理，及Fs大于等于信号最高频率的两倍。2.对信号进行谱分析的重要问题是频谱分辨率D和分析误差，频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现的频谱分辨率是2π/N，因此要求2π/N≤D，可以据此是选择FFT的变换区间N。3.用FFT对连续信号进行谱分析的流程框图如下：三、实验内容时域采样连续时间信号X(t)时域离散信号X(n)截断后的信号xN(n)分析不同长度矩形窗截断对信号的频谱的影响分析不同采样频率对信号频谱的影响时域补零后的信号分析时域补零对信号频谱有何影响矩形窗时域截断时域补零归纳总结开始1.对信号x=cos(8πt)+cos(16πt)+cos(20πt)进行谱分析，将cos(8πt)、cos(16πt)、cos(20πt)分别记作x1、x2、x3，那么x=x1+x2+x3。2.观察信号在采样点数一定，分析不同采样频率对信号频谱的影响。3.观察信号在矩形窗下，不同的截取长度对信号频谱的影响。4.观察信号补零前后频谱的变化，以及补零后信号与截取相同长度数据下信号频谱有何不同。四、实验步骤1.选取采样频率Fs=15、20、32、64Hz，采样点数N一定时，对信号x进行采样，观察采样信号频谱有何不同，并分析原因。2.分别取矩形窗函数的长度N=16、32、64，截取信号后观察信号的频谱有何不同，并分析原因。3.分别取窗函数长度N=15、96截取信号。将N=15点的信号补零至96点，观察其与N=96点的截取信号频谱有何不同，并分析原因。五、实验结果及分析1.不同采样频率下信号的频谱比较⑴取采样频率Fs=15、20、32、64Hz，采样点数N=64固定，信号频谱如下图5.1所示（实现程序见附录1）。-8-6-4-20246802040幅频特性Fs=15频率f/Hz幅度-10-8-6-4-20246810050100幅频特性Fs=20频率f/Hz幅度-20-15-10-50510152002040幅频特性Fs=32频率f/Hz幅度-40-30-20-1001020304002040幅频特性Fs=64频率f/Hz幅度图5.1⑵结果分析从图5.1可知，不同采样频率下信号的频谱是不同的。①当采样频率Fs=15Hz时，经过显然信号频谱出现混叠失真现象，很难观察出信号x频谱所含的频率成分。究其原因，这主要是采样频率Fs=15Hz小于信号x所含频率成分最高频率f3=10Hz的两倍（即20Hz），导致信号频谱产生混叠现象，不能准确地反映出信号x的频谱信息。②当采样频率Fs=20Hz时，较采样频率Fs=15Hz时的频谱来说，此时频谱混叠现象减轻了很多，但还是不能清晰地反映出信号的频谱成分。可见采样频率很是不够高，为了更好地观察信号x的频谱成分应加大采样频率。③当采样频率Fs=32Hz时，信号频谱没有混叠失真的产生，可以很清楚地观察到信号x所包含的三种频率成分。④当Fs=64Hz时，此时采样频率已经足够高了，频谱混叠失真虽然没有了，但信号x中所包含的f2=8Hz与f3=10Hz的频谱成分没有分开，可见此时信号x的频率分辨率降低了。究其原因，虽然采样频率足够高了，但采样点数N=64一定，频谱的固有分辨率F=Fs/N降低，导致这种现象的出现。2.不同长度矩形窗截取信号的频谱比较⑴矩形窗函数长度N=16、32、64下，截取后信号的频谱如下图5.2（实现程序见附录2）。-40-30-20-10010203001020幅频特性N=16频率f/Hz幅度-40-30-20-10010203001020幅频特性N=32频率f/Hz幅度-40-30-20-1001020304002040幅频特性N=64频率f/Hz幅度图5.2⑵结果分析根据图5.2的结果分析如下。①当N=16时，信号x中包含的信号x1、x2、x3的频率成分无法分辨出来，即频率f1,f2,f3无法分辨出来；当N=32时，仅能分辨出信号x中包含的信号x1的频率成分，其他的两根谱线无法分辨出来；当N=64时，信号x中包含的信号x1、x2、x3的频率成分都能够较好地分辨出来。究其原因，信号x中含有x1、x2、x3三种信号的频率成分，其中f1=4Hz,f2=8Hz,f3=10Hz,频率之间的最小间隔minf=f3-f2=2Hz。N=16时，频谱辨率F=Fs/N=64/16=4Hz,可见Fminf。故此时，信号x中的x2与x3的频率成分可以分辨出来，此时F=f2-f1=4Hz,x1与x2的频率差恰好处在临界点上，反映在幅频特性上正如图5.2最上方的图所示，频谱出现两个拐点。②当N=32时，两个频率之间的最小间隔minf=f3-f2=2Hz。频谱分辨率F=Fs/N=64/32=2Hz,可见F=minf。故此时，信号x的频谱可以将f1的频率成分显示出来，而f2与f3的频率间隔恰好处在临界处，故信号x的频谱表现出图5.2中间图的情况。③当N=64时，两个频率之间的最小间隔minf=f3-f2=2Hz。而FFT之后的信号频谱分辨率F=Fs/N=64/64=1Hz,可见Fminf。故此时，信号x的频谱可以将f1、f2与f3的频率成分全都显示出来，x的频谱表现出图5.2下方图的结果。另外，我们也应注意到，N=64时在x的频谱中，虽然x中含有的三种频率成分全都显示了出来，但各频率成分的频谱宽度很大，没有达到理论上的尖脉冲，效果不佳。这主要是因为所取数据长度有限，不能很好的反映出原始信号的特性。也可以理解为，时域所取数据有限相当于时域截断，反映在频域上就会产生频谱扩展、拖尾现象，难以出现理想的尖脉冲效果。综上所述，时域截断使原信号x频谱中的单频谱线展宽，即出现频谱泄露。同时也可以观察到，随着截取长度的增加，频谱泄露越来越小，频率分辨率越高。3.时域补零后信号频谱的比较⑴补零前后信号频谱如下图5.3（实现程序见附录3）。-30-20-10010203001020幅频特性N=15，无补零频率f/Hz幅度-40-30-20-1001020304001020幅频特性N=15,补零81点频率f/Hz幅度-40-30-20-10010203040050幅频特性N=96，无不零频率f/Hz幅度图5.3⑵结果分析从图5.3可知，①当截取信号长度N=15，无不零时，其信号频谱几乎无法显示出相关频谱的信息，即有栅栏效应的出现，如图5.3上方的图所示。②将N=15时截取信号补81个零，可以观察到此时的频谱谱线很密集，但很难看出信号的频谱成分，如图5.3中间图所示。③对原信号加长取样数据，得到N=100的截取信号，可以观察到，此时可以清晰地获取信号的频谱成分，足见此时信号具有较高的频率分辨力，如图5.3下方图所示。究其原因，时域补零之后增加了频域的抽样点数，谱线变密，谱线变密后原来看不到的频谱分量也就有可能看到了。但时域数据的实际长度并没没有改变，也就不能提高频率分辨力，那只有增加数据的有效长度才能提高频率分辨力。六、归纳总结综上所述，用FFT对连续信号进行谱分析时应该注意以下几点要求：1.要对连续信号进行采样，根据奎斯特采样定理确定的采样频率Fs要大于信号最高频率的两倍，一般应取两倍以上，以避免频谱的混叠失真。同时在采样点数一定的情况下，并不算采样频率越高得到的信号频谱越好，还要考虑到频率分辨率的要求。2.要考虑截取时域数据的长度，若信号为周期信号，则所取时域长度应为周期的整数倍；若信号为非周期信号，那么所取信号时域长度应尽量长，减小频谱泄露，以便频谱较好地反映出信号的频谱特性。3.还要考虑到频率分辨力F的大小，F一定要小于信号中不同频率成分最小差值minf。采样频率一定时，F与所取信号的时域长度有直接关系，因此截取的数据长度N是一定要达到F所满足的条件，即F=1/NT≤minf。4.在截取信号长度一定的情况下，为了更好地观察信号的频谱分量(减小栅栏效应的影响)，应该采取时域补零的办法，但时域补零不能提高频率分辨力。时域补零之后增加了频域的抽样点数，谱线变密，谱线变密后原来看不到的频谱分量也就有可能看到了。但时域数据的实际长度并没没有改变，也就不能提高频率分辨力，那只有增加截取数据的有效长度才能提高频率分辨力。附录1%%%--------不同采样频率下信号频谱分析----------%%%clc;clearall;closeall;%%%采样频率Fs=15%%%%%%%%%%%%%%Fs=15;Ts=1/Fs;N1=64;t1=0:Ts:(N1-1)*Ts;x1=cos(8*pi*t1)+cos(16*pi*t1)+cos(20*pi*t1);x_fft1=fftshift(fft(x1,N1));%将零频移至中间f1=-Fs*round((N1-1)/2)/N1:Fs/N1:Fs*fix((N1-1)/2)/N1;%%%%%%采样频率Fs=20%%%%%%%%%%%%%%Fs=20;Ts=1/Fs;N2=64;t2=0:Ts:(N2-1)*Ts;x2=cos(8*pi*t2)+cos(16*pi*t2)+cos(20*pi*t2);x_fft2=fftshift(fft(x2,N2));%将零频移至中间f2=-Fs*round((N2-1)/2)/N2:Fs/N2:Fs*fix((N2-1)/2)/N2;%%%%%%采样频率Fs=32%%%%%%%%%%%%%%Fs=32;Ts=1/Fs;N3=64;t3=0:Ts:(N3-1)*Ts;x3=cos(8*pi*t3)+cos(16*pi*t3)+cos(20*pi*t3);x_fft3=fftshift(fft(x3,N3));%将零频移至中间f3=-Fs*round((N3-1)/2)/N3:Fs/N3:Fs*fix((N3-1)/2)/N3;%%%%%%采样频率Fs=64%%%%%%%%%%%%%%Fs=64;Ts=1/Fs;N4=64;t4=0:Ts:(N4-1)*Ts;x4=cos(8*pi*t4)+cos(16*pi*t4)+cos(20*pi*t4);x_fft4=fftshift(fft(x4,N4));%将零频移至中间f4=-Fs*round((N4-1)/2)/N4:Fs/N4:Fs*fix((N4-1)/2)/N4;%%%%频谱图对比%%%%%%%%%%%%figure;subplot(411);plot(f1,abs(x_fft1));title('幅频特性Fs=15');xlabel('频率f/Hz');ylabel('幅度');gridon;subplot(412);plot(f2,abs(x_fft2));title('幅频特性Fs=20');xlabel('频率f/Hz');ylabel('幅度');gridon;subplot(413);plot(f3,abs(x_fft3));title('幅频特性Fs=32');xlabel('频率f/Hz');ylabel('幅度');gridon;subplot(414);plot(f4,abs(x_fft4));title('幅频特性Fs=64');xlabel('频率f/Hz');ylabel('幅度');gridon;附录2%%%--------矩形窗下不同截取长度信号频谱分析——————%%%clc;clearall;closeall;Fs=64;Ts=1/Fs;%%%截取长度N=16%%%%%%%%%%%%%%N1=16;t1=0:Ts:(N1-1)*Ts;x1=cos(8*pi*t1)+cos(16*pi*t1)+cos(20*pi*t1);x_fft1=fftshift(fft(x1,N1));%将零频移至中间f1=-Fs*round((N1-1)/2)/N1:Fs/N1:Fs*fix((N1-1)/2)/N1;%%%%%%截取长度N=32%%%%%%%%%%%%%