Matlab 计算 FFT 的方法及幅值问题

欢迎转载，但请一定要给出原文链接，标注出处，支持原创！谢谢~http目录1、Matlab FFT 函数介绍2、Matlab FFT 程序3、FFT 程序分析3.1 fs 、N 对 FFT 图像幅值的影响。3.2 直流分量、噪声分量对 FFT 图像的影响3.3 总结1、Matlab FFT 函数介绍 FFT (Fast Fourier Transform) 中文为快速傅里叶变换，作用是将离...

Mr-Ma Technology

43649人浏览 · 2020-04-12 16:20:33

Mr-Ma Technology · 2020-04-12 16:20:33 发布

欢迎转载，但请一定要给出原文链接，标注出处，支持原创！谢谢~
https://blog.csdn.net/qq_29225913/article/details/105467006

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。

FFT 计算公式如下：
Xk 长度与 N 相同。
根据奈科斯特定律，只有 f=fs/2 范围内的信号才是被采样到的有效信号，因此得到的频谱肯定是关于 N/2 对称的（就是只看前一半的波形就好）。
第k点的实际频率的计算为 f(k) = k * (fs / n) — — (横轴的频率范围为：f = n * fs / N;)
X[0] 为直流分量，幅值 = 模值(X[0]) / N
X[k] 为个点的频率分量（除X[0]外），幅值 = 模值(X[k]) / (N / 2)

2、Matlab FFT 程序

先摆上代码，运行一下看效果。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   功能：FFT 运用方法（例程）                  % 
%   作者：Mr-Ma Technology(马健维)             %
%   时间：2020.04.12                           %
%   转载请注明出处                             %
%   https://blog.csdn.net/qq_29225913/article/details/105467006                
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 描述
% 使用 Matlab 的 FFT，有一个注意事项，采样数量是采样频率的整数倍，且最好是 2 的 n 次幂。
% 否则会导致采样分辨率匹配不上，表现为在对应的频点处没有采样到，出现偏移，因此会导致
% 描绘图像时，发现所需要的频点处显示出来的幅度值小于预期的值。
% 可以通过修改参数定义的值得定义来直观的看FFT的效果

%% 参数定义
    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 结果的影响。

1、 fs = 100，N = 100（这个情况下，信号 3 不满足 fs ≥ 2f 的情况，因此信号 3 不会显示出来）
2、增大 N ，fs = 100，N = 150，由于补零了，所以存在能量泄漏，计算出来的幅值会变小。
3、增大Fs，fs = 200，N=150，此时由于 N 与 fs 仍然是非整数倍关系，且样本数量不足够大，导致分辨率不足，使得不能够将所有的频率都能够对应上，造成频率泄漏，导致幅值变小。
4、继续加大 N ，直接扩大10倍，fs = 200，N = 1500，可以观察到随着样本数量增加的时候，幅频特性曲线会变窄，与实际的频率越接近，即精度越高。但是由于 fs 与 N 仍然非整数倍关系，因此仍然会有上一点存在的问题。
5、这一步，减少 N ，但是设置为 fs 的整数倍（往上滑，与第1点对比），fs = 100， N = 200，可以看出，与第 1 点对比的时候，幅值依然是实际值，但是精度更高了（由于N 的数多了一倍）
6、（对照组）扩大 N 为 fs 的 2 倍（与第 3 、5 点对比）fs = 200，N= 400，因为同样 N 是 fs 的整数倍，幅值为实际值。
7、最后，设置 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值 ≤ 实际采样点）默认信息是随机分布，最后的一节信号不会对整体的频谱产生较大的影响。（若要考虑全部采样点，则可忽略对实际幅值的影响，实际幅值不一定需要）
叠加噪声的分析，可以转换到频域上进行处理。