Python 是目前的熱門語(yǔ)言，一直覺得掌握一門編程語(yǔ)言對(duì)作為搞技術(shù)的來(lái)說(shuō)還是很有必要的，結(jié)合工作中能用到的一些數(shù)據(jù)處理和分析的內(nèi)容，覺得從數(shù)據(jù)分析入手，爭(zhēng)取能夠掌握Python在數(shù)據(jù)處理領(lǐng)域的一些應(yīng)用。下面是基于Python的numpy進(jìn)行的數(shù)字信號(hào)的頻譜分析介紹

一、傅里葉變換

傅里葉變換是信號(hào)領(lǐng)域溝通時(shí)域和頻域的橋梁，在頻域里可以更方便的進(jìn)行一些分析。傅里葉主要針對(duì)的是平穩(wěn)信號(hào)的頻率特性分析，簡(jiǎn)單說(shuō)就是具有一定周期性的信號(hào)，因?yàn)楦道锶~變換采取的是有限取樣的方式，所以對(duì)于取樣長(zhǎng)度和取樣對(duì)象有著一定的要求。

二、基于Python的頻譜分析

# _*_ coding:utf-8 _*_

import numpy as np #導(dǎo)入一個(gè)數(shù)據(jù)處理的模塊

import pylab as pl #導(dǎo)入一個(gè)繪圖模塊，matplotlib下的模塊

sampling_rate = 8000 ##取樣頻率

fft_size =512 #FFT處理的取樣長(zhǎng)度

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

#np.arange(起點(diǎn)，終點(diǎn)，間隔)產(chǎn)生1s長(zhǎng)的取樣時(shí)間

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

#兩個(gè)正弦波疊加，156.25HZ和234.375HZ，因此如上面簡(jiǎn)單

#的介紹FFT對(duì)于取樣時(shí)間有要求，

#N點(diǎn)FFT進(jìn)行精確頻譜分析的要求是N個(gè)取樣點(diǎn)包含整數(shù)個(gè)

#取樣對(duì)象的波形。

#因此N點(diǎn)FFT能夠完美計(jì)算頻譜對(duì)取樣對(duì)象的要求

#是n*Fs/N（n*采樣頻率/FFT長(zhǎng)度），

#因此對(duì)8KHZ和512點(diǎn)而言，

#完美采樣對(duì)象的周期最小要求是8000/512=15.625HZ,

#所以156.25的n為10,234.375的n為15。

xs = x[:fft_size]# 從波形數(shù)據(jù)中取樣fft_size個(gè)點(diǎn)進(jìn)行運(yùn)算

xf = np.fft.rfft(xs)/fft_size # 利用np.fft.rfft()進(jìn)行FFT計(jì)算，rfft()是為了更方便

#對(duì)實(shí)數(shù)信號(hào)進(jìn)行變換，由公式可知/fft_size為了正確顯示波形能量

# rfft函數(shù)的返回值是N/2+1個(gè)復(fù)數(shù)，分別表示從0(Hz)

#到sampling_rate/2(Hz)的分。

#于是可以通過(guò)下面的np.linspace計(jì)算出返回值中每個(gè)下標(biāo)對(duì)應(yīng)的真正的頻率：

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

# np.linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None)

#在指定的間隔內(nèi)返回均勻間隔的數(shù)字

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

#最后我們計(jì)算每個(gè)頻率分量的幅值，并通過(guò) 20*np.log10()

#將其轉(zhuǎn)換為以db單位的值。為了防止0幅值的成分造成log10無(wú)法計(jì)算，

#我們調(diào)用np.clip對(duì)xf的幅值進(jìn)行上下限處理

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

pl.subplot(211)

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

pl.xlabel(u"時(shí)間(秒)")

pl.title(u"The Wave and Spectrum 156.25Hz234.375Hz")

pl.subplot(212)

pl.plot(freqs, xfp)

pl.xlabel(u"Hz")

pl.subplots_adjust(hspace=0.4)

pl.show()

#繪圖顯示結(jié)果

現(xiàn)在來(lái)看看頻譜泄露，將采樣對(duì)象的頻率改變

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

我們明顯看出，第一個(gè)對(duì)象的頻譜分析出現(xiàn)“泄露”，能量分散到其他頻率上，

沒法準(zhǔn)確計(jì)算到計(jì)算對(duì)象的頻譜特性。

窗函數(shù)

上面我們可以看出可以通過(guò)加“窗”函數(shù)的方法來(lái)處理，盡量保證FFT長(zhǎng)度內(nèi)

的取樣對(duì)象是對(duì)稱的。

import pylab as pl

import scipy.signal as signal

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

pl.plot(signal.hann(512))#漢明窗函數(shù)

pl.show()

對(duì)上述出現(xiàn)頻譜泄露的函數(shù)進(jìn)行加窗處理，后面會(huì)介紹一下各種加窗函數(shù)的原理和效果。