Matlab FFT计算方法及幅值问题详解
Matlab FFT计算方法及幅值问题详解
快速傅里叶变换（FFT）是将时域信号转换到频域的重要工具，在信号处理、音频分析等领域有着广泛的应用。本文将详细介绍如何在Matlab中计算FFT，并探讨采样频率、采样点数等参数对FFT结果的影响。
1. Matlab FFT函数介绍
FFT（Fast Fourier Transform）中文为快速傅里叶变换，其作用是将离散的时域信号变换到频域。在一般情况下，时域信号都比较难看出特征，但是转换成频域后，就比较容易看出来，因此，很多信号分析都会采用FFT变换，然后对信号进行频谱分析。
Matlab中，FFT函数的语法如下：
- Y = fft(X)
- Y = fft(X,n)
- Y = fft(X,n,dim)
推荐直接在命令行使用help fft看到更加详细的描述信息与例程，这里主要记录一下用法以及注意事项。
N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方，但是这里直接取读取出来的样本数作为N。
2. Matlab FFT程序示例
下面是一个使用Matlab进行FFT计算的示例程序：
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 功能：FFT 运用方法（例程） %
% 时间：2020.04.12 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 参数定义
A1 = 2; % 信号1 幅值（A）
f1 = 11; % 信号1 频率（f）
A2 = 0.5; % 信号2 幅值（A)
f2 = 29; % 信号2 频率（f)
A3 = 1; % 信号3 幅值（A)
f3 = 51; % 信号3 频率（f)
Dc = 0; % 直流分量
Noise = 0; % 噪声大小
fs = 128; % 采样频率 （fs）
N = 1024; % 采样点数（样本数量）
%% 运算（生成时域曲线、FFT计算）
n = 0:N-1; % 等差生成序列，
t = n / fs; % 时间序列，用于下式
RA = rand(1,N); % 随机噪声（0~1）
RA = RA - mean(RA); % 减去噪声均值（-0.5~0.5）,去除噪声中的直流分量。
x = Dc + A1 * sin(2*pi*f1*t) + A2 * sin(2*pi*f2*t) + A3 * sin(2*pi*f3*t)+ Noise * RA ; % 两个正弦信号、直流分量、噪声 相叠加
y = abs(fft(x,N)); % 对上式进行 N 点 FFT 计算 ，并取模值
A = y * 2 / N; % 模值转换为幅值
A(1) = A(1) / 2; % 根据公式 ，X[0] 不用乘以 2
f = n * fs / N; % 转换为频率区间
%% 绘图
subplot(2,1,1);
plot(t,x);
xlabel('时间/s');ylabel('振幅');
title('时域曲线')
subplot(2,1,2);
plot(f(1:N/2),A(1:N/2)); % 描绘图像
xlabel('频率/Hz');ylabel('振幅');
title(['幅频特性曲线 :Fs = ',num2str(fs),', N = ',num2str(N)])
运行结果如下图所示：
3. FFT程序分析
通过运行上面的代码，可以得到与上图一致的运行效果，该程序中，原函数是由直流分量、信号1、信号2、信号3、噪声分量叠加而成的（上图中暂时将直流分量与噪声分量设置为0），在实际工程使用中，将例程中生成的原函数替换成需要进行分析的序列就可以了。
上面的图片中，每一个信号的幅值就刚好等于所设置的振幅，此处有一个地方需要注意的就是，一般情况下，采样点N都不固定，都是随输入的因素影响，例如对一首采样率为44.1KHz的歌来分析，读取出来的数据长度N就不一定是2的n次方了，如果直接送入进行FFT计算，到N点结算数据，在FFT内由于使用蝶形运算，所以是需要N是2的n次方，因此会导致在N个采样点后方补零， 可能在某些频点出看到的赋值，会比预想值要小，这就是能量泄露。如果只需要考虑各频谱值得比例（即只看频率分布情况），不考虑实际FFT后的幅值的话，可以直接传入N进去，如果需要进行频谱分析，就要要考虑实际的幅值对分析结果的影响，此处就可以选择加窗，还有一种情况就是N足够大，并且是均匀分布，就可以考虑截取N的值，使得满足N等于2的n次方。
3.1 fs、N对FFT图像幅值的影响
上述程序中，fs = 128；N = 4096；（N是fs的整数倍，且N是2的n次方），下面来分析下不同的fs与N对实际FFT结果的影响。
- fs = 100，N = 100（这个情况下，信号3不满足fs ≥ 2f的情况，因此信号3不会显示出来）
- 增大N，fs = 100，N = 150，由于补零了，所以存在能量泄漏，计算出来的幅值会变小。
- 增大Fs，fs = 200，N=150，此时由于N与fs仍然是非整数倍关系，且样本数量不足够大，导致分辨率不足，使得不能够将所有的频率都能够对应上，造成频率泄漏，导致幅值变小。
- 继续加大N，直接扩大10倍，fs = 200，N = 1500，可以观察到随着样本数量增加的时候，幅频特性曲线会变窄，与实际的频率越接近，即精度越高。但是由于fs与N仍然非整数倍关系，因此仍然会有上一点存在的问题。
- 这一步，减少N，但是设置为fs的整数倍（往上滑，与第1点对比），fs = 100，N = 200，可以看出，与第1点对比的时候，幅值依然是实际值，但是精度更高了（由于N的数多了一倍）
- （对照组）扩大N为fs的2倍（与第3、5点对比）fs = 200，N= 400，因为同样N是fs的整数倍，幅值为实际值。
- 最后，设置fs、N均为2的n次方，这里的fs根据最高信号频率来选择即可，比如信号3是频率最高的信号（51Hz），因此fs选择128，N尽量的大，因此N选择4096，此时时域信号中的频率信息，都能够很好的表示出来。
3.2 直流分量、噪声分量对FFT图像的影响
在这一节考虑将直流分量，与噪声信号加入，通过改变以下两项来加入。其他值使用默认的参数(fs = 128，N = 1024)
- 设置直流分量为1，可以看到FFT图像x = 0处有一个直流分量。(实际使用中，直流分量对信号的分析用处不大，所以一般都会将此值设为0，或者将 用
x = x - mean(x)
的方法，将直流分量过滤掉) - 按照直流分量的计算公式，需要在程序里对该值进行修改（幅值 = 模值 / N ； 当 x = 0)
- 修改噪声分量的值为10，（即叠加一个幅度为±5的噪声信号），使用噪声前，需要对噪声进行消除直流处理。
- 可以看出来，在时域的图像上，已经看不清信号的内容了，但是在频域上，仍然能够分辨的出有3个主要信号，并且在采样范围内，均匀分布着噪声信号。后面可以通过设计带通滤波器，将3个有用信号的频段提取出来，其他频段的信号进行过滤，在利用IFFT还原出降噪后的波形。
3.3 总结
- 在能够自己选定fs的情况下，尽量选择fs为2的n次方，且N取fs的倍数，这样可以在幅频特性曲线上，得到与实际相符的幅度值。
- fs固定的情况下，（如音频文件，44.1K/48K/…），不能修改fs，则选择N为fs的整数倍，通常这种情况下N的数量都会足够大，因此N向前截取（即N值 ≤ 实际采样点）默认信息是随机分布，最后的一节信号不会对整体的频谱产生较大的影响。（若要考虑全部采样点，则可忽略对实际幅值的影响，实际幅值不一定需要）
- 叠加噪声的分析，可以转换到频域上进行处理。