Python: FFT的输入与输出分析

目的

本文简述如何使用numpy的fft lib进行快速傅里叶变换，以及对fft变换后结果的分析。由于水平有限，不当之处望指正。

数据生成

使用如下代码生成仿吉他C和弦声音片段，如下代码所示。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

import librosa
import IPython.display as ipd

sr = 8000
l = 3

t = sr * l
t = np.array(range(0, t)) / sr
x = np.sin(2 * np.pi * 130.81 * t) + 0.8 * np.sin(2 * np.pi * 164.81 * t) + 1.2 * np.sin(2 * np.pi * 196.00 * t)
x = x / max(x)

ipd.Audio(x, rate=sr)

上述代码中：

• 和声的C和弦由C3(130.81HZ), E3(164.81HZ), G3(196.00HZ)组成。
• G3的音量(1.2)>C3的音量(1.0)>E3(0.8)的音量。
• 声音的总长度为3秒钟(l变量)。
• 声音的采样率为8000HZ(sr变量)。
• 合成的声音进行归一化处理。

FFT

numpy的FFT变换非常简单，使用fft()函数即可。

y = np.fft.fft(x)

所以，如何读懂FFT的返回结果y呢？使用如下代码将FFT的返回结果y映射到频率上即可：

fr = np.array(range(0, len(x))) / l
region = (int)(len(fr) / 2)
plt.plot(fr[0:region], abs(y)[0:region])

• 首先，type(y)的结果是numpy.darray，由于fft输出的对称性，y[0:1/2] 与y[1/2:1]的数据是冗余的。因此只取结果的前半部分即可。
• 其次，type(y[0])的结果是numpy.complex128的复数。实部(real)和虚部(imaginary)通过计算可以得到y[0]频率下的振幅(Magnitude)和相位(Phase)[2]。
• y[0]频率是频率为[ 0, 8000 / len(fr) )的频率范围。以此类推，y[1]频率为[ 8000 / len(y), 2 * 8000 / len(y) )的频率范围。
• 如果将上图放大看，或者只打印fr, abs(y)中的一部分数据，则上图中的突起部分正好落在C3, E3以及G3频率处。
• abs(y)只打印复数数组y中的实部。

分辨率

上述使用3秒的数据(24000个点)分析频率特征比较消耗计算时间，也可以通过更少的数据来分析，如下述代码使用1024个采样点。由于采样点变少，而分析的频率范围(0 ~ 8000HZ)不变，因此输出结果的精度降低了。

xslice = x[0:1024]
frslice = np.array(range(0, len(xslice))) / (1024 / 8000)

yslice = np.fft.fft(xslice)
cslice = abs(yslice)
plt.plot(frslice[0:100], cslice[0:100])

总结

• numpy的fft使用很方便，输入一维数组，输出同样大小的一维复数数组。
• t数组是信号x的时域范围，上例中是0 ~ 3秒。
• fr数组是变换结果y的频率范围，上例中是0 ~ 8000HZ。

引用

[1] https://www.youtube.com/watch?v=HEcdmDSVHsE&list=PL6rD7wR_ngfnpLFdhbjCTTwRY8PNk6qph&index=3
[2] https://stackoverflow.com/questions/25624548/fft-real-imaginary-abs-parts-interpretation