FFT 变换的 MATLAB 实现
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 个点、第 50 个点、
第 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 结果以及
上面的分析计算,我们就可以写出信号的表达式了,它就是我们开始
提供的信号。
本测试数据使用的 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/18
0);
%显示原始信号
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('相位-频率曲线图');