深入理解fft变换-SciDog

Fourier 变换是牛人傅里叶提出的一种数学变换，可以用来做频率分析。现在利用matlab已经可以简单快速的实现一个fft变换了，但是要正确的理解结果，还是有几个点需要特别注意：看下面这个比较傻的matlab fft程序：

clc; clear; close all;
fs = 1000;
t = 0: 1/fs: 1-1/fs;
f1=2;
x1 = 0.5*cos(2*pi*f1*t +0.2);
subplot(2,1,1)
plot(x1)
xlabel('Sample');
ylabel('Amplitude');
X1 = fft(x1);
subplot(2,1,2)
plot([0:length(X1)-1], abs(X1));
ylabel('magnitude'); xlabel('Bins');
angle(X1(3))

可以发现几点：

一个单频信号，提取出了2各频点；
fft变换前后时间域和频率域的点数是相同的，这里都是1000点；

但是这个结果，我们并不满意：我们希望得到的频率分析结果是在频率f=2处，有一个幅度为1的频谱，显然现在的结果并不对。我们将程序调整以下，fft后的幅度除以N，然后在将零频率之外的幅度谱乘以2，得到正确的单边频谱如下：

clc; clear; close all;
fs = 1000;
t = 0: 1/fs: 1-1/fs;
f1=2;
x1 = 0.5*cos(2*pi*f1*t +0.2);
subplot(2,1,1)
plot(x1)
xlabel('time');
ylabel('Amplitude');

X1 = fft(x1);
L = length(X1);
P2 = abs(X1/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(L/2))/L;

subplot(2,1,2)
plot(f, P1);
ylabel('magnitude'); xlabel('Frequency');

结果似乎是对的，还比较满意，但是，我们把频率改成f=2.5Hz，结果如下：

由于频率分辨率是1Hz，对于2.5Hz的信号无法分辨，出现了频谱泄露，原来2.5处的单频泄露到周围的频率范围内。

怎么解决这个问题？

通过zero padding 补零的方法提高频率分辨率（频率分辨率=fs/N=1/tlen），采样频率不变的情况下，增加N，可以在最高频率不变的情况下提高频率分辨率；

x1 = 0.5*cos(2*pi*f1*t +0.2);
x1 = [x1, zeros(1,11000)];

主频附近的主瓣能量变强，但是引入和很多旁瓣的干扰。通过对信号进行加窗的方法，我们可以去除旁瓣的干扰。

clc; clear; close all;
fs = 1000;
t = 0: 1/fs: 1-1/fs;
f1=2.5;
x1_org = 0.5*cos(2*pi*f1*t +0.2);
x1 = x1_org.*hanning(length(x1_org))';

hanwin = hanning(1000);
plot(hanwin);
hold on 
plot(x1,'r');
hold on
plot(x1_org,'g')
x1 = [x1, zeros(1,11000)];

subplot(2,1,1)
plot(x1)
xlabel('time');
ylabel('Amplitude'); title('f=2.5Hz')

X1 = fft(x1);
L = length(X1);
P2 = abs(X1/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(L/2))/L;

subplot(2,1,2)
plot(f, P1);
ylabel('magnitude'); xlabel('Frequency');

这时我们通过 zero-padding + hanning window的方法提高了对2.5Hz频率的分辨率，但是，我们注意到幅度比0.5 小。下面通过考察2个正弦信号叠加的信号，继续考虑频率分辨率的问题：

clc; clear; close all;
fs = 1000;
t = 0: 1/fs: 1-1/fs;
f1=2.5;
f2=3;
x1_org = 0.5*cos(2*pi*f1*t +0.2)+1*cos(2*pi*f2*t+0.8);
x1=x1_org;
% x1 = [x1, zeros(1,11000)];

subplot(2,1,1)
plot(x1)
xlabel('time');
ylabel('Amplitude'); title('f=2.5Hz')

X1 = fft(x1);
L = length(X1);
P2 = abs(X1/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(L/2))/L;

subplot(2,1,2)
plot(f, P1);
ylabel('magnitude'); xlabel('Frequency');

可以看到，频域分辨率是1Hz，不足以分辨，2.5Hz 和 3Hz 的两个正弦信号：

在补零以后，得到下面的曲线：仍然不能去区分2.5Hz 和 3Hz，为什么？补零后的频率分辨率达到1000/12000=0.0833Hz，但是我们注意到，补零只能增加频域的视觉分辨率，时域补零相当于频域插值，这也就是“栅栏效应“。而波形的频域分辨率由原始数据的长度决定，并不随补零发生改变，仍然是1Hz。

要想区分这两个信号，需要增加原始数据的时长；