FFT应用之Python实践

FFT的输入和输出：

输入2^n的numpy数组，输出2^n的numpy complex类型数组

>>> x = np.random.rand(8)
>>> x
array([ 0.15562099,  0.56862756,  0.54371949,  0.06354358,  0.60678158,
        0.78360968,  0.90116887,  0.1588846 ])
>>> xf = np.fft.fft(x)
>>> xf
array([ 3.78195634+0.j        , -0.53575962+0.57688097j,
       -0.68248579-1.12980906j, -0.36656155-0.13801778j,
        0.63262552+0.j        , -0.36656155+0.13801778j,
       -0.68248579+1.12980906j, -0.53575962-0.57688097j])

FFT的逆变换，变换回来

>>> np.fft.ifft(xf)
array([ 0.15562099 +0.00000000e+00j,  0.56862756 +1.91940002e-16j,
        0.54371949 +1.24900090e-16j,  0.06354358 -2.33573365e-16j,
        0.60678158 +0.00000000e+00j,  0.78360968 +2.75206729e-16j,
        0.90116887 -1.24900090e-16j,  0.15888460 -2.33573365e-16j])

 

FFT的输出特点

• 下标为0和N/2的两个复数的虚数部分为0，
• 下标为i和N-i的两个复数共轭，也就是其虚数部分数值相同、符号相反。

 

FFT的有用信息

由于虚数部共轭和虚数部为0等规律，真正有用的信息保存在下标从0到N/2的N/2+1个虚数中， 又由于下标为0和N/2的值虚数部分为0，因此只有N个有效的实数值。

 

FFT的物理含义

• 首先下标为0的实数表示了时域信号中的直流成分的多少
• 下标为i的复数a+b*j表示时域信号中周期为N/i个取样值的正弦波和余弦波的成分的多少， 其中a表示cos波形的成分，b表示sin波形的成分

 

例子：

让我们通过几个例子来说明一下，下面是对一个直流信号进行FFT变换：

>>> x = np.ones(8)
>>> x
array([ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.])
>>> np.fft.fft(x)/len(x)
array([ 1.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j])

所谓直流信号，就是其值不随时间变化，因此我们创建一个值全为1的数组x，我们看到它的FFT结果除了下标为0的数值不为0以外，其余的都为0。(为了计算各个成分的能量多少，需要将FFT的结果除以FFT的长度)，这表示我们的时域信号是直流的，并且其成分为1。

下面我们产生一个周期为8个取样的正弦波，然后观察其FFT的结果：

 

如何用linspace创建取样点x

要计算周期为8个取样的正弦波，就需要把0-2*pi的区间等分为8分，如果用np.linspace(0, 2*np.pi, 8)的话，产生的值为：

>>> np.linspace(0, 2*np.pi, 8)
array([ 0.        ,  0.8975979 ,  1.7951958 ,  2.6927937 ,  3.5903916 ,
    4.48798951,  5.38558741,  6.28318531])
>>> 2*np.pi / 0.8975979
7.0000000079986666

这没有等分8分，实际只等分7分，因为包括了最后一个2*pi的点。

如果要正确的等分8份

>>> np.linspace(0, 2*np.pi, 9, endpoint=False)
array([ 0.        ,  0.6981317 ,  1.3962634 ,  2.0943951 ,  2.7925268 ,
        3.4906585 ,  4.1887902 ,  4.88692191,  5.58505361])

FFT的结果

>>> np.fft.fft(y)/len(y)
array([  1.42979161e-18 +0.00000000e+00j,
        -4.44089210e-16 -5.00000000e-01j,
         1.53075794e-17 -1.38777878e-17j,
         3.87737802e-17 -1.11022302e-16j,
         2.91853672e-17 +0.00000000e+00j,
         0.00000000e+00 -9.71445147e-17j,
         1.53075794e-17 +1.38777878e-17j,   3.44085112e-16 +5.00000000e-01j])

让我们再来看看对正弦波的FFT的计算结果吧。可以看到下标为1的复数的虚数部分为-0.5，而我们产生的正弦波的放大系数(振幅)为1，它们之间的关系是-0.5*（-2）=1。再来看一下余弦波形:

>>> np.fft.fft(np.cos(x))/len(x)
array([ -4.30631550e-17 +0.00000000e+00j,
         5.00000000e-01 -2.52659764e-16j,
         1.53075794e-17 +0.00000000e+00j,
         1.11022302e-16 +1.97148613e-16j,
         1.24479962e-17 +0.00000000e+00j,
        -1.11022302e-16 +1.91429446e-16j,
         1.53075794e-17 +0.00000000e+00j,   5.00000000e-01 -1.35918295e-16j])

只有下标为1的复数的实数部分有有效数值0.5，和余弦波的放大系数1之间的关系是0.5*2=1。再来看2个例子：

>>> np.fft.fft(2*np.sin(2*x))/len(x)
array([  6.12303177e-17 +0.00000000e+00j,
         6.12303177e-17 +6.12303177e-17j,
        -1.83690953e-16 -1.00000000e+00j,
         6.12303177e-17 -6.12303177e-17j,
         6.12303177e-17 +0.00000000e+00j,
         6.12303177e-17 +6.12303177e-17j,
        -1.83690953e-16 +1.00000000e+00j,   6.12303177e-17 -6.12303177e-17j])
>>> np.fft.fft(0.8*np.cos(2*x))/len(x)
array([ -2.44921271e-17 +0.00000000e+00j,
        -3.46370983e-17 +2.46519033e-32j,
         4.00000000e-01 -9.79685083e-17j,
         3.46370983e-17 -3.08148791e-32j,
         2.44921271e-17 +0.00000000e+00j,
         3.46370983e-17 -2.46519033e-32j,
         4.00000000e-01 +9.79685083e-17j,  -3.46370983e-17 +3.08148791e-32j])

上面产生的是周期为4个取样(N/2)的正弦和余弦信号，其FFT的有效成分在下标为2的复数中，其中正弦波的放大系数为2,因此频域虚数部分的值为-1；余弦波的放大系数为0.8，因此其对应的值为0.4。

同频率的正弦波和余弦波通过不同的系数叠加，可以产生同频率的各种相位的余弦波，因此我们可以这样来理解频域中的复数:

• 复数的模（绝对值）代表了此频率的余弦波的振幅
• 复数的辐角代表了此频率的余弦波的相位

让我们来看最后一个例子：

>>> x = np.arange(0, 2*np.pi, 2*np.pi/128)
>>> y = 0.3*np.cos(x) + 0.5*np.cos(2*x+np.pi/4) + 0.8*np.cos(3*x-np.pi/3)
>>> yf = np.fft.fft(y)/len(y)
>>> yf[:4]
array([  1.00830802e-17 +0.00000000e+00j,
         1.50000000e-01 +6.27820821e-18j,
         1.76776695e-01 +1.76776695e-01j,   2.00000000e-01 -3.46410162e-01j])
>>> np.angle(yf[1])
4.1854721366992471e-017
>>> np.abs(yf[1]), np.rad2deg(np.angle(yf[1]))
(0.15000000000000008, 2.3980988870246962e-015)
>>> np.abs(yf[2]), np.rad2deg(np.angle(yf[2]))
(0.25000000000000011, 44.999999999999993)
>>> np.abs(yf[3]), np.rad2deg(np.angle(yf[3]))
(0.39999999999999991, -60.000000000000085)

在这个例子中我们产生了三个不同频率的余弦波，并且给他们不同的振幅和相位：

• 周期为128/1.0点的余弦波的相位为0， 振幅为0.3
• 周期为128/2.0点的余弦波的相位为45度， 振幅为0.5
• 周期为128/3.0点的余弦波的相位为-60度，振幅为0.8

对照yf[1], yf[2], yf[3]的复数振幅和辐角，我想你应该对FFT结果中的每个数值所表示的意思有很清楚的理解了吧。

 

总结: 在生成离散时数域,有两个关键的参数,一个是取的点数,一个是周期,或者解释为采样频率或者固有频率.

 

例子二: 合成三角波信号

方法是先生成三角波信号,然后对三角波信号做FFT,然后将得到的结果,用三角函数,再来表示三角波信号

# 产生size点取样的三角波，其周期为1
def triangle_wave(size):
    x = np.arange(0, 1, 1.0/size)
    y = np.where(x<0.5, x, 0)
    y = np.where(x>=0.5, 1-x, y)
    return x, y

fft_size = 256

 

# 计算三角波和其FFT
x, y = triangle_wave(fft_size)
fy = np.fft.fft(y) / fft_size

 

# 绘制三角波的FFT的前20项的振幅，由于不含下标为偶数的值均为0， 因此取
# log之后无穷小，无法绘图，用np.clip函数设置数组值的上下限，保证绘图正确
pl.figure()
pl.plot(np.clip(20*np.log10(np.abs(fy[:20])), -120, 120), "o")
pl.xlabel("frequency bin")
pl.ylabel("power(dB)")
pl.title("FFT result of triangle wave")

由频谱信号重新组合成时域信号

# 取FFT计算的结果freqs中的前n项进行合成，返回合成结果，计算loops个周期的波形
def fft_combine(freqs, n, loops=1):
    length = len(freqs) * loops
    data = np.zeros(length)
    index = loops * np.arange(0, length, 1.0) / length * (2 * np.pi)
    for k, p in enumerate(freqs[:n]):
        if k != 0: p *= 2 # 除去直流成分之外，其余的系数都*2
        data += np.real(p) * np.cos(k*index) # 余弦成分的系数为实数部
        data -= np.imag(p) * np.sin(k*index) # 正弦成分的系数为负的虚数部
    return index, data   

# 绘制原始的三角波和用正弦波逐级合成的结果，使用取样点为x轴坐标
pl.figure()
pl.plot(y, label="original triangle", linewidth=2)
for i in [0,1,3,5,7,9]:
    index, data = fft_combine(fy, i+1, 2)  # 计算两个周期的合成波形
    pl.plot(data, label = "N=%s" % i)
pl.legend()
pl.title("partial Fourier series of triangle wave")
pl.show()

 

同理,处理方波

 

总结

计算机只能处理离散的数据,看起来连续的结果实际上都是离散计算出来的。离散傅里叶变换很重要的概念是采样。理解离散傅里叶变换的物理意义更重要。FFT只是计算DFT的一种算法。