对fft的一些理解

例子1：对50Hz的正弦波做fft分析，默认加矩形窗closeall;clc;clear;%%Fs=1e3;%SamplingfrequencyT=1/Fs;%SampletimeL=1000;%Lengthofsignalt=(0:L-1)*T;%Timevectorx=1*sin(2*pi*50*t)+0*sin(2*pi*50.5*t);%Signalwin=rectwin(L);%Windownx_data=x.*win';%plotSignalfigure('Position',[65278560420]);plot(t(1:L),x_data(1:L))gridon;title('TimeDomainSignal')xlabel('time(s)')%fftNFFT=1*2^nextpow2(L);%Nextpowerof2fromlengthofxX=fft(x_data,NFFT)/L;f=(Fs/2)*linspace(0,1,NFFT/2+1);%Plotsingle-sidedamplitudespectrum.figure('Position',[640278560420]);%plot(f,2*abs(X(1:NFFT/2+1)))plot(f,20*log10(2*abs(X(1:NFFT/2+1))))gridon;title('Single-SidedAmplitudeSpectrumofSignal')xlabel('Frequency(Hz)')ylabel('|Y(f)|(dB)')%plotTimeDomainWindowfigure('Position',[6578560420]);plot(t(1:L),win(1:L))gridon;title('TimeDomainWindow')xlabel('time(s)')%plotFrequencyDomainWindowWIN=fft(win,NFFT)/L;figure('Position',[64078560420]);plot(f,20*log10(1*abs(WIN(1:NFFT/2+1))))gridon;title('Single-SidedAmplitudeSpectrumofWindow')xlabel('Frequency(Hz)')由于采用fft算法，对原始L=1000点数据进行补零处理，得到NFFT=1024点数据，故频率分辨率f0=(1000/1024)/1s=0.9765625Hz，图中125Hz处的点其实是第128个频率点（0.9765625Hz*128=125Hz）由于栅栏效应，本来应该体现出1Hz的频谱分辨率（注意区别于频率分辨率），上图并未体现出1Hz的频谱分辨率将fft的数据点补零扩大到32*1024点NFFT=32*2^nextpow2(L);%Nextpowerof2fromlengthofx此时，频谱分辨率已经能分辨到1Hz而且频率分辨率f0=(1000/1024/32)/1s=0.9765625Hz/32，图中125Hz处的点其实是第128*32个频率点（0.9765625Hz/32*32*128=125Hz）将fft的数据点补零扩大到64*1024点NFFT=64*2^nextpow2(L);%Nextpowerof2fromlengthofx可以看到，频谱分辨率依然为1Hz，并没有因为补零而增加频谱分辨率（可以理解为实际有用的信号点并没有增加，只是增加了冗余的补零点）但是频率分辨率增加到了f0=(1000/1024/64)/1s=0.9765625Hz/64，图中125Hz处的点其实是第128*64个频率点（0.9765625Hz/64*64*128=125Hz）上述可见，补零可以无限增加fft频率分辨率（注意不是频谱分辨率），但是在频谱已经达到最高频谱分辨率后（理论上频谱最高分辨率为信号长度的倒数，上述例子为1Hz），补零并不能提高频谱分辨率，过高的频率分辨率会带来巨大的计算量，没有意义。怎么能增加频谱分辨率，也即怎么增加我们肉眼能看到的一根一根的频谱呢？增加信号有效长度！例子2：对50Hz的正弦波+50.2Hz正弦波做fft分析，默认加矩形窗closeall;clc;clear;%%Fs=1e3;%SamplingfrequencyT=1/Fs;%SampletimeL=1000;%Lengthofsignalt=(0:L-1)*T;%Timevectorx=1*sin(2*pi*50*t)+1*sin(2*pi*50.2*t);%Signalwin=rectwin(L);%Windownx_data=x.*win';%plotSignalfigure('Position',[65278560420]);plot(t(1:L),x_data(1:L))gridon;title('TimeDomainSignal')xlabel('time(s)')%fftNFFT=32*2^nextpow2(L);%Nextpowerof2fromlengthofxX=fft(x_data,NFFT)/L;f=(Fs/2)*linspace(0,1,NFFT/2+1);%Plotsingle-sidedamplitudespectrum.figure('Position',[640278560420]);%plot(f,2*abs(X(1:NFFT/2+1)))plot(f,20*log10(2*abs(X(1:NFFT/2+1))))gridon;title('Single-SidedAmplitudeSpectrumofSignal')xlabel('Frequency(Hz)')ylabel('|Y(f)|(dB)')%plotTimeDomainWindowfigure('Position',[6578560420]);plot(t(1:L),win(1:L))gridon;title('TimeDomainWindow')xlabel('time(s)')%plotFrequencyDomainWindowWIN=fft(win,NFFT)/L;figure('Position',[64078560420]);plot(f,20*log10(1*abs(WIN(1:NFFT/2+1))))gridon;title('Single-SidedAmplitudeSpectrumofWindow')xlabel('Frequency(Hz)')可以看到，虽然频率分辨率f0=(1000/1024/32)/1s=0.9765625Hz/32，但是由于信号有效长度为1s，频谱分辨率仍然为1Hz，无法分辨0.2Hz的信号，此时必须增加有效信号长度L=50000;%Lengthofsignal信号长度变成了50s，故频谱分辨率变成了0.02Hz，能够分辨出0.2Hz的信号了！实际上从时域波形也能看出端倪，分析1s信号的时候，看起来似乎还是50Hz的信号，根本看不出有其他信号。上述的另一个可见的问题是，频谱泄露！50Hz的原始信号，进行fft变换后，出现了其他频率的信号！原因是采样时没有采用相关采样。回顾例子1中的图，原始代码中，Fs=1000Hz，L=1000，但是后面做fft时，由于fft算法必须采用2^N个点，所以在后面做了补零操作，但这样并不满足信号的相关采样，导致频谱泄露！所谓coherentsampling，就是相关采样，或者相干采样。进行FFT变换的时候，如果不采取相关采样(coherentsampling)的话，会产生一些问题，如频谱泄漏，为了避免频谱泄漏一般有两种方法，一种就是采取相关采样，另外一种则需要用到窗函数进行补偿，减少频谱泄漏，但是不能完全消除，一般窗函数有不同种类，可根据实际情况选择适合的窗函数。关于频谱泄露和相关采样，可以参考：例子3：对50Hz正弦波做fft分析，采取相关采样，不加窗closeall;clc;clear;%%Fs=1.024e3;%SamplingfrequencyT=1/Fs;%SampletimeL=1024*32;%Lengthofsignalt=(0:L-1)*T;%Timevectorx=1*sin(2*pi*50*t)+0*sin(2*pi*50.2*t);%Signalwin=rectwin(L);%Windownx_data=x.*win';%plotSignalfigure('Position',[65278560420]);plot(t(1:L),x_data(1:L))gridon;title('TimeDomainSignal')xlabel('time(s)')%fftNFFT=2^nextpow2(L);%Nextpowerof2fromlengthofxX=fft(x_data,NFFT)/L;f=(Fs/2)*linspace(0,1,NFFT/2+1);%Plotsingle-sidedamplitudespectrum.figure('Position',[640278560420]);%plot(f,2*abs(X(1:NFFT/2+1)))plot(f,20*log10(2*abs(X(1:NFFT/2+1))))gridon;title('Single-SidedAmplitudeSpectrumofSignal')xlabel('Frequency(Hz)')ylabel('|Y(f)|(dB)')%plotTimeDomainWindowfigure('Position',[6578560420]);plot(t(1:L),win(1:L))gridon;title('TimeDomainWindow')xlabel('time(s)')%plotFrequencyDomainWindowWIN=fft(win,NFFT)/L;figure('Position',[64078560420]);plot(f,20*log10(1*abs(WIN(1:NFFT/2+1))))gridon;title('Single-SidedAmplitudeSpectrumofWindow')xlabel('Frequency(Hz)')注意，上述代码中，为了做fft分析时，避免对有效信号再进行补零操作而违背相关采样，取Fs=1.024e3Hz，L=1024的整数倍（这里取32）可以清楚看到，频谱图中只出现了50Hz的频谱，没有其他无关频谱！相关采样成功！例子4：对50Hz的正弦波做fft分析，不采取相关采样，加汉宁窗closeall;clc;clear;%%Fs=1e3;%SamplingfrequencyT=1/Fs;%SampletimeL=1000*32;%Lengthofsignalt=(0:L-1)*T;%Timevectorx=1*sin(2*pi*50*t)+0*sin(2*pi*50.2*t);%Signalwin=hann(L);%Windownx_data=x.*win';%plotSignalfigure('Position',[65278560420]);plot(t(1:L),x_data(1:L))gridon;title('TimeDomainSignal')xlabel('time(s)')%fftNFFT=2^nextpow2(L);%Nextpowerof2fromlengthofxX=fft(x_data,NFFT)/L;f=(Fs/2)*linspace(0,1,NFFT/2+1);%Plotsingle-sidedamplitudespectrum.figure('Position',[640278560420]);%plot(f,2*abs(X(1:NFFT/2+1)))plot(f,20*lo