对fft的一些理解

例子1：对50Hz的正弦波做fft分析，默认加矩形窗

close all;

clc;

clear;

%%

Fs = 1e3; % Sampling frequency

T = 1/Fs; % Sample time

L = 1000; % Length of signal

t = (0:L-1)*T; % Time vector

x = 1*sin(2*pi*50*t)+0*sin(2*pi*50.5*t); % Signal

win = rectwin(L); % Windown

x_data = x .* win';

% plot Signal

figure('Position',[65 278 560 420]);

plot(t(1:L),x_data(1:L))

grid on;

title('Time Domain Signal')

xlabel('time (s)')

% fft

NFFT = 1*2^nextpow2(L); % Next power of 2 from length of x

X = fft(x_data,NFFT)/L;

f = (Fs/2)*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.

figure('Position',[0 278 560 420]);

% plot(f,2*abs(X(1:NFFT/2+1)))

plot(f,20*log10(2*abs(X(1:NFFT/2+1))))

grid on;

title('Single-Sided Amplitude Spectrum of Signal')

xlabel('Frequency (Hz)')

ylabel('|Y(f)| (dB)')

% plot TimeDomain Window

figure('Position',[65 78 560 420]);

plot(t(1:L),win(1:L))

grid on;

title('Time Domain Window')

xlabel('time (s)')

% plot FrequencyDomain Window

WIN = fft(win,NFFT)/L;

figure('Position',[0 78 560 420]);

plot(f,20*log10(1*abs(WIN(1:NFFT/2+1))))

grid on;

title('Single-Sided Amplitude Spectrum of Window')

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); % Next power of 2 from length of x

此时，频谱分辨率已经能分辨到1Hz

而且频率分辨率f0=(1000/1024/32)/1s= 0.9765625Hz/32，图中125Hz处的点其实是第128*32个频率点（0.9765625Hz/32*32*128=125Hz）

将fft的数据点补零扩大到*1024点

NFFT = *2^nextpow2(L); % Next power of 2 from length of x

可以看到，频谱分辨率依然为1Hz，并没有因为补零而增加频谱分辨率（可以理解为实际有用的信号点并没有增加，只是增加了冗余的补零点）

但是频率分辨率增加到了f0=(1000/1024/)/1s= 0.9765625Hz/，图中125Hz处的点其实是第128*个频率点（0.9765625Hz/**128=125Hz）

上述可见，补零可以无限增加fft频率分辨率（注意不是频谱分辨率），但是在频谱已经达到最高频谱分辨率后（理论上频谱最高分辨率为信号长度的倒数，上述例子为1Hz），补零并不能提高频谱分辨率，过高的频率分辨率会带来巨大的计算量，没有意义。

怎么能增加频谱分辨率，也即怎么增加我们肉眼能看到的一根一根的频谱呢？增加信号有效长度！

例子2：对50Hz的正弦波+50.2Hz正弦波做fft分析，默认加矩形窗

close all;

clc;

clear;

%%

Fs = 1e3; % Sampling frequency

T = 1/Fs; % Sample time

L = 1000; % Length of signal

t = (0:L-1)*T; % Time vector

x = 1*sin(2*pi*50*t)+1*sin(2*pi*50.2*t); % Signal

win = rectwin(L); % Windown

x_data = x .* win';

% plot Signal

figure('Position',[65 278 560 420]);

plot(t(1:L),x_data(1:L))

grid on;

title('Time Domain Signal')

xlabel('time (s)')

% fft

NFFT = 32*2^nextpow2(L); % Next power of 2 from length of x

X = fft(x_data,NFFT)/L;

f = (Fs/2)*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.

figure('Position',[0 278 560 420]);

% plot(f,2*abs(X(1:NFFT/2+1)))

plot(f,20*log10(2*abs(X(1:NFFT/2+1))))

grid on;

title('Single-Sided Amplitude Spectrum of Signal')

xlabel('Frequency (Hz)')

ylabel('|Y(f)| (dB)')

% plot TimeDomain Window

figure('Position',[65 78 560 420]);

plot(t(1:L),win(1:L))

grid on;

title('Time Domain Window')

xlabel('time (s)')

% plot FrequencyDomain Window

WIN = fft(win,NFFT)/L;

figure('Position',[0 78 560 420]);

plot(f,20*log10(1*abs(WIN(1:NFFT/2+1))))

grid on;

title('Single-Sided Amplitude Spectrum of Window')

xlabel('Frequency (Hz)')

可以看到，虽然频率分辨率f0=(1000/1024/32)/1s= 0.9765625Hz/32，但是由于信号有效长度为1s，频谱分辨率仍然为1Hz，无法分辨0.2Hz的信号，此时必须增加有效信号长度

L = 50000; % Length of signal

信号长度变成了50s，故频谱分辨率变成了0.02Hz，能够分辨出0.2Hz的信号了！

实际上从时域波形也能看出端倪，分析1s信号的时候，看起来似乎还是50Hz的信号，根本看不出有其他信号。

上述的另一个可见的问题是，频谱泄露！50Hz的原始信号，进行fft变换后，出现了其他频率的信号！

原因是采样时没有采用相关采样。

回顾例子1中的图，原始代码中，Fs=1000Hz，L=1000，但是后面做fft时，由于fft算法必须采用2^N个点，所以在后面做了补零操作，但这样并不满足信号的相关采样，导致频谱泄露！

所谓coherent sampling，就是相关采样，或者相干采样。进行FFT变换的时候，如果不采取相关采样(coherent sampling)的话，会产生一些问题，如频谱泄漏，为了避免频谱泄漏一般有两种方法，一种就是采取相关采样，另外一种则需要用到窗函数进行补偿，减少频谱泄漏，但是不能完全消除，一般窗函数有不同种类，可根据实际情况选择适合的窗函数。

关于频谱泄露和相关采样，可以参考：http://www.vfe.cc/NewsDetail-876.aspx

例子3：对50Hz正弦波做fft分析，采取相关采样，不加窗

close all;

clc;

clear;

%%

Fs = 1.024e3; % Sampling frequency

T = 1/Fs; % Sample time

L = 1024*32; % Length of signal

t = (0:L-1)*T; % Time vector

x = 1*sin(2*pi*50*t)+0*sin(2*pi*50.2*t); % Signal

win = rectwin(L); % Windown

x_data = x .* win';

% plot Signal

figure('Position',[65 278 560 420]);

plot(t(1:L),x_data(1:L))

grid on;

title('Time Domain Signal')

xlabel('time (s)')

% fft

NFFT = 2^nextpow2(L); % Next power of 2 from length of x

X = fft(x_data,NFFT)/L;

f = (Fs/2)*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.

figure('Position',[0 278 560 420]);

% plot(f,2*abs(X(1:NFFT/2+1)))

plot(f,20*log10(2*abs(X(1:NFFT/2+1))))

grid on;

title('Single-Sided Amplitude Spectrum of Signal')

xlabel('Frequency (Hz)')

ylabel('|Y(f)| (dB)')

% plot TimeDomain Window

figure('Position',[65 78 560 420]);

plot(t(1:L),win(1:L))

grid on;

title('Time Domain Window')

xlabel('time (s)')

% plot FrequencyDomain Window

WIN = fft(win,NFFT)/L;

figure('Position',[0 78 560 420]);

plot(f,20*log10(1*abs(WIN(1:NFFT/2+1))))

grid on;

title('Single-Sided Amplitude Spectrum of Window')

xlabel('Frequency (Hz)')

注意，上述代码中，为了做fft分析时，避免对有效信号再进行补零操作而违背相关采样，取Fs=1.024e3Hz，L = 1024的整数倍（这里取32）

可以清楚看到，频谱图中只出现了50Hz的频谱，没有其他无关频谱！相关采样成功！

例子4：对50Hz的正弦波做fft分析，不采取相关采样，加汉宁窗

close all;

clc;

clear;

%%

Fs = 1e3; % Sampling frequency

T = 1/Fs; % Sample time

L = 1000*32; % Length of signal

t = (0:L-1)*T; % Time vector

x = 1*sin(2*pi*50*t)+0*sin(2*pi*50.2*t); % Signal

win = hann(L); % Windown

x_data = x .* win';

% plot Signal

figure('Position',[65 278 560 420]);

plot(t(1:L),x_data(1:L))

grid on;

title('Time Domain Signal')

xlabel('time (s)')

% fft

NFFT = 2^nextpow2(L); % Next power of 2 from length of x

X = fft(x_data,NFFT)/L;

f = (Fs/2)*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.

figure('Position',[0 278 560 420]);

% plot(f,2*abs(X(1:NFFT/2+1)))

plot(f,20*log10(2*abs(X(1:NFFT/2+1))))

grid on;

title('Single-Sided Amplitude Spectrum of Signal')

xlabel('Frequency (Hz)')

ylabel('|Y(f)| (dB)')

% plot TimeDomain Window

figure('Position',[65 78 560 420]);

plot(t(1:L),win(1:L))

grid on;

title('Time Domain Window')

xlabel('time (s)')

% plot FrequencyDomain Window

WIN = fft(win,NFFT)/L;

figure('Position',[0 78 560 420]);

plot(f,20*log10(1*abs(WIN(1:NFFT/2+1))))

grid on;

title('Single-Sided Amplitude Spectrum of Window')

xlabel('Frequency (Hz)')

可见，加窗可以抑制频谱泄露，但是不能完全消除，完全消除频谱泄露的唯一办法，还是得采用相关采样！！

例子5：对50Hz的正弦波做fft分析，采取相关采样，并加汉宁窗

close all;

clc;

clear;

%%

Fs = 1.024e3; % Sampling frequency

T = 1/Fs; % Sample time

L = 1024*32; % Length of signal

t = (0:L-1)*T; % Time vector

x = 1*sin(2*pi*50*t)+0*sin(2*pi*50.2*t); % Signal

win = hann(L); % Windown

x_data = x .* win';

% plot Signal

figure('Position',[65 278 560 420]);

plot(t(1:L),x_data(1:L))

grid on;

title('Time Domain Signal')

xlabel('time (s)')

% fft

NFFT = 2^nextpow2(L); % Next power of 2 from length of x

X = fft(x_data,NFFT)/L;

f = (Fs/2)*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.

figure('Position',[0 278 560 420]);

% plot(f,2*abs(X(1:NFFT/2+1)))

plot(f,20*log10(2*abs(X(1:NFFT/2+1))))

grid on;

title('Single-Sided Amplitude Spectrum of Signal')

xlabel('Frequency (Hz)')

ylabel('|Y(f)| (dB)')

% plot TimeDomain Window

figure('Position',[65 78 560 420]);

plot(t(1:L),win(1:L))

grid on;

title('Time Domain Window')

xlabel('time (s)')

% plot FrequencyDomain Window

WIN = fft(win,NFFT)/L;

figure('Position',[0 78 560 420]);

plot(f,20*log10(1*abs(WIN(1:NFFT/2+1))))

grid on;

title('Single-Sided Amplitude Spectrum of Window')

xlabel('Frequency (Hz)')

由于窗函数本身就会带来频谱泄露，如果采取了相关采样，就不要再加窗了，用默认矩形窗即可。

Single-Sided Amplitude Spectrum of Signal0-50-100-150|Y(f)| (dB)-200-250-300-350-4000100200300Frequency (Hz)4005006个人理解，共同学习进步，如有谬误，望指出！

参考文献：

1. matlat help

2. 《数字信号处理》 （美）Richard G Lyons著；朱光明，程建远，刘保童等译

lifusu

2015.04.12