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()