基于matlab的FFT分析

离散傅里叶变换DFT的计算公式如下，FFT为DFT的一种快速算法。

matlab代码

%===============
%   N = 64时
%===============
fs = 100;   %采样频率
N = 64;    %数据点数
n = 0 : N-1;
% 抽样间隔Ts=1/fs,所以t=n*Ts=n/fs为时间序列
t = n / fs; %时间序列
x = 0.5*sin(2*pi*15*t) + 2*sin(2*pi*40*t);
y = fft(x, N);  %对信号进行64点快速傅里叶变换
mag1 = abs(y);   %求傅里叶变换后的振幅
% 在fs(100Hz)内做N(64)点FFT,频率分辨率为fs/N=1.5625
f = n * fs / N; %频率序列
subplot(2,2,1);
plot(f, mag1,'LineWidth',1);   %绘制随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');
title('N=64');
subplot(2,2,2);
plot(f(1:N/2), mag1(1:N/2)); %绘制奈奎斯特频率之间随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');
title('N=64');

%===============
%   N = 1024时
%===============
fs = 100;   %采样频率
N = 1024;    %数据点数
n = 0 : N-1;
% 抽样间隔Ts=1/fs,所以t=n*Ts=n/fs为时间序列
t = n / fs; %时间序列
x = 0.5*sin(2*pi*15*t) + 2*sin(2*pi*40*t);
y = fft(x, N);  %对信号进行1024点快速傅里叶变换
mag2 = abs(y);   %求傅里叶变换后的振幅
% 在fs(100Hz)内做N(1024)点FFT,频率分辨率为fs/N=0.097656
f = n * fs / N; %频率序列
subplot(2,2,3);
plot(f, mag2,'LineWidth',1);   %绘制随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');
title('N=1024');
subplot(2,2,4);
plot(f(1:N/2), mag2(1:N/2)); %绘制奈奎斯特频率之前随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');
title('N=1024');

FFT结果具有对称性，通常我们只用前半部分的结果，也就是小于采样频率一半的结果。同时，只有在采样频率一半以内、具有一定幅值的信号频率才是真正的信号频率。

原信号为
$$0.5sin(2\ast pi\ast 15\ast t)+2sin(2\ast pi\ast 40\ast t)$$
f1=15Hz，f2=40Hz，最高频率fm为40Hz，根据奈奎斯特抽样定理，抽样频率fs必须大于等于原信号最高频率的2倍，即 fs>=2fm。本例中，fs取100Hz。

说明（一）

由上图可以看出，N越大，频谱就越逼真。且只有在f1=15Hz，f2=40Hz时，其模值才会出现峰值，而其他频率点上模值接近于0。
我们分别来查看N=64和N=128时15Hz和40Hz处的模值。频率f = n * fs / N，所以n = f * N / fs。
N = 64，f1 = 15Hz时n = 15 * 64 /100 =9.6，向上取整加1为11

N = 64，f2 = 40Hz时n = 40 * 64 /100 =25.6，向上取整加1为27

N = 1024，f1 = 15Hz时n = 15 * 1024 /100 =153.6，向上取整加1为155

N = 1024，f2 = 40Hz时n = 40 * 1024 /100 =409.6，向上取整加1为411

说明（二）

另外，细心的读者会发现（图中有标注）N=1024时频率比N=64时更接近50Hz。这是由于N=64时的频率分辨率为fs / N = 100 / 64 = 1.5625，而N = 1024时的频率分辨率为fs / N = 100 / 1024 = 0.097656，N = 1024时的频率分辨率明显高于N = 64时的频率分辨率。