使用Python进行FFT频谱分析
2015-12-29 02:21阅读：
FFT频谱分析原理
采样定理：采样频率要大于信号频率的两倍。
N个采样点经过FFT变换后得到N个点的以复数形式记录的FFT结果。
假设采样频率为Fs，采样点数为N。那么FFT运算的结果就是N个复数（或N个点），每一个复数就对应着一个频率值以及该频率信号的幅值和相位。第一个点对应的频率为0Hz（即直流分量），最后一个点N的下一个点对应采样频率Fs。其中任意一个采样点n所代表的信号频率：
Fn=(n-1)*Fs/N。
这表明，频谱分析得到的信号频率最大为(N-1)*Fs/N,对频率的分辨能力是Fs/N。采样频率和采样时间制约着通过FFT运算能分析得到的信号频率上限，同时也限定了分析得到的信号频率的分辨率。
每一个复数的模值对应该点所对应的频率值的幅度特性，具体的定量关系如下：
假设信号由以下周期的原始信号叠加而成：
Y = A1 + A2 Cos
(2*PI*ω2*t + φ2 * PI/180)
+ A3 Cos (2*PI*ω3*t + φ3
* PI/180)
那么，在经过FFT分析后得到的第一个点的模值是A1的N倍，而且只有在FFT结果点对应的频率在ω2,ω3时，其模值才明显放大，在其他频率点，模值接近于0。在这些模值明显放大的点中，除第一个点之外的其它点模值是相应信号幅值的N/2倍。
每个复数的相位就是在该频率值下信号的相位：φ2，φ3。
FFT结果有对称性，通常我们只是用前半部分的结果，也就是小于采样频率一半的结果。同时也只有采样频率一半以内、具有一定幅值的信号频率才是真正的信号频率。
Python实践FFT频谱分析
假如信号S是有1个直流信号和4个周期信号叠加而成，如下公式所列（t为自变量，pi为圆周率值）现要对其进行FFT分析并绘制频谱图。
S = 2.0 + 3.0 *
cos(2.0 * p
i
*
50
*
t
-
pi
*
30/180)
+
1.5
*
cos(2.0
*
pi
*
75
*
t
+
pi
*
90/180)
+
1.0
*
cos(2.0
*
pi
*
150
*
t
+
pi
*
120/180)
+
2.0
*
cos(2.0
*
pi
*
220
*
t
+
pi
*
30/180)
我们先使用
Python绘制其1秒内的波形图：
import
numpy
as
np
import
pylab
as
pl
import
math
#
采样步长
t
=
[x/1048.0
for
x
in
range(1048)]
#
设计的采样值
y
=
[2.0
+
3.0
*
math.cos(2.0
*
math.pi
*
50
*
t0
-
math.pi
*
30/180)
+
1.5
*
math.cos(2.0
*
math.pi
*
75
*
t0
+
math.pi
*
90/180)
+
1.0
*
math.cos(2.0
*
math.pi
*
150
*
t0
+
math.pi
*
120/180)
+
2.0
*
math.cos(2.0
*
math.pi
*
220
*
t0
+
math.pi
*
30/180)
for
t0
in
t
]
pl.plot(t,y)
pl.show()
波形图如图1所示。
图
1
信号
S在1秒内的波形
现在我们要对该信号在
0-1秒时间内进行频谱分析，我们在这1秒时间内采样1048点，则我们的采样频率为1048Hz。这两个设定决定了我们频谱分析所能得到的最高信号频率是1047Hz
(1048Hz*(1048-1)/1048)，但只有小于524Hz的频率才可能是真实的信号频率。我们能分辨的最小频率差为1Hz（1048Hz/1048)。输入如下的python代码：
N=len(t)
#
采样点数
fs=1048.0
#
采样频率
df
=
fs/(N-1)
#
分辨率
f
=
[df*n
for
n
in
range(0,N)]
#
构建频率数组
使用下面的代码进行快速傅里叶变换并对结果模值进行缩放：
Y
=
np.fft.fft(y)*2/N
#*2/N
反映了FFT变换的结果与实际信号幅值之间的关系
absY
=
[np.abs(x)
for
x
in
Y]
#求傅里叶变换结果的模
此时我们得到了信号的FFT分析结果，它存储在列表Y内，同时我们也获取了Y内每一个元素（是复数）的模组成的列表absY。根据理论分析，absY内的元素大部分都极其接近0，只有极少数峰值提示的是信号频率的幅度。我们写一小段代码显示幅度较大的信号的信息。代码如下：
i=0
while
i
<
len(absY):
if
absY[i]>0.01:
print('freq:M,
Y:
%5.2f
+
%5.2f
j,
A:%3.2f,
phi:
%6.1f'
%(i,
Y[i].real,
Y[i].imag,
absY[i],math.atan2(Y[i].imag,Y[i].real)*180/math.pi))
i+=1
该段代码的执行结果如下：
freq:
0,
Y:
4.00
+
0.00
j,
A:4.00,
phi:
0.0
freq:
50,
Y:
2.60
+
-1.50
j,
A:3.00,
phi:
-30.0
freq:
75,
Y:
0.00
+
1.50
j,
A:1.50,
phi:
90.0
freq:
150,
Y:
-0.50
+
0.87
j,
A:1.00,
phi:
120.0
freq:
220,
Y:
1.73
+
1.00
j,
A:2.00,
phi:
30.0
freq:
828,
Y:
1.73
+
-1.00
j,
A:2.00,
phi:
-30.0
freq:
898,
Y:
-0.50
+
-0.87
j,
A:1.00,
phi:
-120.0
freq:
973,
Y:
0.00
+
-1.50
j,
A:1.50,
phi:
-90.0
freq:
998,
Y:
2.60
+
1.50
j,
A:3.00,
phi:
30.0
该段结果很好的反映了原始信号的频率、幅度和相位，我们在代码中已经对FFT结果与信号幅值做了一个换算2/N，但在0Hz处，实际的信号应该是FFT结果的1/N，也就是A：4.00再除以2。
使用如下代码绘制频谱图，得到图2：
pl.plot(f,absY)
pl.show()
图
2
信号
S的频谱分析图
由图2可见，频谱分析的结果具有对称性，同220Hz对称的频率为828Hz，可以计算出对称轴为524Hz。但只有小于524Hz的频率才是信号频率。我们通过修改采样频率和采样点数，我们发现作为对称轴的频率会变化，同样的对称轴右侧的信号频率也会发生变化。提示右侧的频率没有实际意义，这也从一个方面印证了采样定理。
本文大量参考了博文：
利用matlab怎样进行频谱分析
在此对博主表示感谢，原文链接：http://blog.sina.cn/dpool/blog/s/blog_a07f4fe301013gj3.html?vt=4.
