快速傅里叶变换python_快速傅里叶变换(FFT) python实现

python代码

# coding=utf-8

from matplotlib.font_manager import FontProperties

import pylab as pl

import numpy as np

chinese_font = FontProperties(fname='/usr/share/fonts/truetype/wqy/wqy-microhei.ttc')

sampling_rate = 8000

fft_size = 512

t = np.arange(0, 1.0, 1.0/sampling_rate)

x = np.sin(2*np.pi*156.25*t) + 2*np.sin(2*np.pi*234.375*t)

xs = x[:fft_size]

xf = np.fft.rfft(xs)/fft_size

freqs = np.linspace(0, sampling_rate/2, fft_size/2+1)

xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e100))

pl.figure(figsize=(8,4))

pl.subplot(211)

pl.plot(t[:fft_size], xs)

pl.xlabel(u"时间(秒)",fontproperties=chinese_font)

pl.title(u"156.25Hz和234.375Hz的波形和频谱", fontproperties=chinese_font)

pl.subplot(212)

pl.plot(freqs, xfp)

pl.xlabel(u"频率(Hz)",fontproperties=chinese_font)

pl.subplots_adjust(hspace=0.4)

pl.show()

详解

因为绘图的时候如果图表中有中文乱码情况，可以设置成系统的中文字体：

chinese_font = FontProperties(fname='/usr/share/fonts/truetype/wqy/wqy-microhei.ttc')

定义了两个常数：sampling_rate, fft_size，分别表示数字信号的取样频率和FFT的长度：

sampling_rate = 8000

fft_size = 512

用取样时间数组t可以很方便地调用函数计算出波形数据，这里计算的是两个正弦波的叠加，一个频率是156.25Hz，一个是234.375Hz，为什么选择这两个奇怪的频率呢？因为这两个频率的正弦波在512个取样点中正好有整数个周期。满足这个条件波形的FFT结果能够精确地反映其频谱：

x = np.sin(2*np.pi*156.25*t) + 2*np.sin(2*np.pi*234.375*t)

N点FFT能精确计算的频率：

假设取样频率为fs, 取波形中的N个数据进行FFT变换。那么这N点数据包含整数个周期的波形时，FFT所计算的结果是精确的。于是能精确计算的波形的周期是: n*fs/N。对于8kHz取样，512点FFT来说，8000/512.0 = 15.625Hz，前面的156.25Hz和234.375Hz正好是其10倍和15倍。

从波形数据x中截取fft_size个点进行fft计算。np.fft库中提供了一个rfft函数，它方便我们对实数信号进行FFT计算。根据FFT计算公式，为了正确显示波形能量，还需要将rfft函数的结果除以fft_size：

xs = x[:fft_size]

xf = np.fft.rfft(xs)/fft_size

rfft函数的返回值是N/2+1个复数，分别表示从0(Hz)到sampling_rate/2(Hz)的N/2+1点频率的成分。于是可以通过下面的np.linspace计算出返回值中每个下标对应的真正的频率：

freqs = np.linspace(0, sampling_rate/2, fft_size/2+1)

最后我们计算每个频率分量的幅值，并通过 20*np.log10() 将其转换为以db单位的值。为了防止0幅值的成分造成log10无法计算，我们调用np.clip对xf的幅值进行上下限处理：

xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e100))

剩下的程序就是将时域波形和频域波形绘制出来，这里就不再详细叙述了。此程序的输出为：