python频谱分析

import numpy as np
import matplotlib.pyplot as pl


sampling_rate = 8000

fft_size = 512

# 首先定义了两个常数：sampling_rate, fft_size，分别表示数字信号的取样频率和FFT的长度.由于快速离散傅里叶算法的影响这里的N必须是2的整数次幂，也称为FFT_SIZE。如果你不是取的整数次，小于则信号补0，大于则信号截断.

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)

#取样时间数组t可以很方便地调用函数计算出波形数据，这里计算的是两个正弦波的叠加，一个频率是156.25Hz，一个是234.375Hz.

用取样时间数组t可以很方便地调用函数计算出波形数据，这里计算的是两个正弦波的叠加，一个频率是156.25Hz，一个是234.375Hz：
对于8KHZ取样，512点FFT来说，8000/512=15.625Hz,156.25Hz和234.375Hz正好是其10倍和15倍。

xs = x[:fft_size]

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

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

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

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

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

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

pl.figure(figsize=(8,4))
pl.subplot(211)
pl.plot(t[:fft_size], xs)
pl.xlabel(u"时间(秒)")
pl.title(u"156.25Hz和234.375Hz的波形和频谱")
pl.subplot(212)
pl.plot(freqs, xfp)
pl.xlabel(u"频率(Hz)")
pl.subplots_adjust(hspace=0.4)
pl.show()