深入浅出解释 FFT（一）—— 用 fft 求频谱_频谱的幅度怎么计算

FFT 是离散傅立叶变换的快速算法，可以将一个信号变换
到频域。有些信号在时域上是很难看出什么特征的，但是如
果变换到频域之后，就很容易看出特征了。这就是很多信号
分析采用 FFT 变换的原因。另外，FFT 可以将一个信号的频谱
提取出来，这在频谱分析方面也是经常用的。 
虽然很多人都知道 FFT 是什么，可以用来做什么，怎么去
做，但是却不知道 FFT 之后的结果是什意思、如何决定要使用
多少点来做 FFT。

现在圈圈就根据实际经验来说说 FFT 结果的具体物理意义。
一个模拟信号，经过 ADC 采样之后，就变成了数字信号。采样
定理告诉我们，采样频率要大于信号频率的两倍，这些我就
不在此罗嗦了。

采样得到的数字信号，就可以做 FFT 变换了。N 个采样点，
经过 FFT 之后，就可以得到 N 个点的 FFT 结果。为了方便进行 FFT
运算，通常 N 取 2 的整数次方。

假设采样频率为 Fs，信号频率 F，采样点数为 N。那么 FFT
之后结果就是一个为 N 点的复数。每一个点就对应着一个频率
点。这个点的模值，就是该频率值下的幅度特性。具体跟原始
信号的幅度有什么关系呢？假设原始信号的峰值为 A，那么 FFT
的结果的每个点（除了第一个点直流分量之外）的模值就是 A
的 N/2 倍。而第一个点就是直流分量，它的模值就是直流分量
的 N 倍。而每个点的相位呢，就是在该频率下的信号的相位。
第一个点表示直流分量（即 0Hz），而最后一个点 N 的再下一个
点（实际上这个点是不存在的，这里是假设的第 N+1 个点，也
可以看做是将第一个点分做两半分，另一半移到最后）则表示
采样频率 Fs，这中间被 N-1 个点平均分成 N 等份，每个点的频率
依次增加。例如某点 n 所表示的频率为：Fn=(n-1)*Fs/N。
由上面的公式可以看出，Fn 所能分辨到频率为为 Fs/N，如果
采样频率 Fs 为 1024Hz，采样点数为 1024 点，则可以分辨到 1Hz。
1024Hz 的采样率采样 1024 点，刚好是 1 秒，也就是说，采样 1 秒
时间的信号并做 FFT，则结果可以分析到 1Hz，如果采样 2 秒时
间的信号并做 FFT，则结果可以分析到 0.5Hz。如果要提高频率
分辨力，则必须增加采样点数，也即采样时间。频率分辨率和
采样时间是倒数关系。
假设 FFT 之后某点 n 用复数 a+bi 表示，那么这个复数的模就是
An = 根号 a*a+b*b，相位就是 Pn=atan2 (b,a)。根据以上的结果，
就可以计算出 n 点（n≠1，且 n<=N/2）对应的信号的表达式为：
An/(N/2)*cos (2*pi*Fn*t+Pn)，即 2*An/N*cos (2*pi*Fn*t+Pn)。
对于 n=1 点的信号，是直流分量，幅度即为 A1/N。
由于 FFT 结果的对称性，通常我们只使用前半部分的结果，
即小于采样频率一半的结果。

好了，说了半天，看着公式也晕，下面圈圈以一个实际的
信号来做说明。

假设我们有一个信号，它含有 2V 的直流分量，频率为 50Hz、
相位为 - 30 度、幅度为 3V 的交流信号，以及一个频率为 75Hz、
相位为 90 度、幅度为 1.5V 的交流信号。用数学表达式就是如下：

S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)

式中 cos 参数为弧度，所以 - 30 度和 90 度要分别换算成弧度。
我们以 256Hz 的采样率对这个信号进行采样，总共采样 256 点。
按照我们上面的分析，Fn=(n-1)*Fs/N，我们可以知道，每两个
点之间的间距就是 1Hz，第 n 个点的频率就是 n-1。我们的信号
有 3 个频率：0Hz、50Hz、75Hz，应该分别在第 1 个点、第 51 个点、
第 76 个点上出现峰值，其它各点应该接近 0。实际情况如何呢？
我们来看看 FFT 的结果的模值如图所示。



图 1 FFT 结果
从图中我们可以看到，在第 1 点、第 51 点、和第 76 点附近有
比较大的值。我们分别将这三个点附近的数据拿上来细看：
1 点： 512+0i
2 点： -2.6195E-14 - 1.4162E-13i 
3 点： -2.8586E-14 - 1.1898E-13i

50 点：-6.2076E-13 - 2.1713E-12i
51 点：332.55 - 192i
52 点：-1.6707E-12 - 1.5241E-12i

75 点：-2.2199E-13 -1.0076E-12i
76 点：3.4315E-12 + 192i
77 点：-3.0263E-14 +7.5609E-13i

很明显，1 点、51 点、76 点的值都比较大，它附近的点值
都很小，可以认为是 0，即在那些频率点上的信号幅度为 0。
接着，我们来计算各点的幅度值。分别计算这三个点的模值，
结果如下：
1 点： 512
51 点：384
76 点：192
按照公式，可以计算出直流分量为：512/N=512/256=2；
50Hz 信号的幅度为：384/(N/2)=384/(256/2)=3；75Hz 信号的
幅度为 192/(N/2)=192/(256/2)=1.5。可见，从频谱分析出来
的幅度是正确的。
然后再来计算相位信息。直流信号没有相位可言，不用管
它。先计算 50Hz 信号的相位，atan2 (-192, 332.55)=-0.5236,
结果是弧度，换算为角度就是 180*(-0.5236)/pi=-30.0001。再
计算 75Hz 信号的相位，atan2 (192, 3.4315E-12)=1.5708 弧度，
换算成角度就是 180*1.5708/pi=90.0002。可见，相位也是对的。
根据 FFT 结果以及上面的分析计算，我们就可以写出信号的表达
式了，它就是我们开始提供的信号。

总结：假设采样频率为 Fs，采样点数为 N，做 FFT 之后，某
一点 n（n 从 1 开始）表示的频率为：Fn=(n-1)*Fs/N；该点的模值
除以 N/2 就是对应该频率下的信号的幅度（对于直流信号是除以
N）；该点的相位即是对应该频率下的信号的相位。相位的计算
可用函数 atan2 (b,a) 计算。atan2 (b,a) 是求坐标为 (a,b) 点的角
度值，范围从 - pi 到 pi。要精确到 xHz，则需要采样长度为 1/x 秒
的信号，并做 FFT。要提高频率分辨率，就需要增加采样点数，
这在一些实际的应用中是不现实的，需要在较短的时间内完成
分析。解决这个问题的方法有频率细分法，比较简单的方法是
采样比较短时间的信号，然后在后面补充一定数量的 0，使其长度
达到需要的点数，再做 FFT，这在一定程度上能够提高频率分辨力。
具体的频率细分法可参考相关文献。

[附录：本测试数据使用的 matlab 程序]
close all; % 先关闭所有图片
Adc=2; % 直流分量幅度
A1=3; % 频率 F1 信号的幅度
A2=1.5; % 频率 F2 信号的幅度
F1=50; % 信号 1 频率 (Hz)
F2=75; % 信号 2 频率 (Hz)
Fs=256; % 采样频率 (Hz)
P1=-30; % 信号 1 相位 (度)
P2=90; % 信号相位 (度)
N=256; % 采样点数
t=[0:1/Fs:N/Fs]; % 采样时刻

% 信号
S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);
% 显示原始信号
plot(S);
title (' 原始信号 ');

figure;
Y = fft (S,N); % 做 FFT 变换
Ayy = (abs (Y)); % 取模
plot (Ayy (1:N)); % 显示原始的 FFT 模值结果
title ('FFT 模值 ');

figure;
Ayy=Ayy/(N/2); % 换算成实际的幅度
Ayy(1)=Ayy(1)/2;
F=([1:N]-1)*Fs/N; % 换算成实际的频率值
plot (F (1:N/2),Ayy (1:N/2)); % 显示换算后的 FFT 模值结果
title (' 幅度 - 频率曲线图 ');

figure;
Pyy=[1:N/2];
for i="1:N/2"
Pyy (i)=phase (Y (i)); % 计算相位
Pyy (i)=Pyy (i)*180/pi; % 换算为角度
end;
plot (F (1:N/2),Pyy (1:N/2)); % 显示相位图
title (' 相位 - 频率曲线图 ');