姓名:韓宜真
學號:17020120095
轉載自:https://mp.weixin.qq.com/s/uTVhAWR5N1ZgoxmWL4Zj_g
【嵌牛導讀】本文介紹了一種用python進行數(shù)字信號處理的方法
【嵌牛鼻子】CPU計算能力基于python
【嵌牛提問】如何用python做信號處理?
【嵌牛正文】
Python 是目前的熱門語言,一直覺得掌握一門編程語言對作為搞技術的來說還是很有必要的喘沿,結合工作中能用到的一些數(shù)據(jù)處理和分析的內(nèi)容,覺得從數(shù)據(jù)分析入手艺挪,爭取能夠掌握Python在數(shù)據(jù)處理領域的一些應用被丧。下面是基于Python的numpy進行的數(shù)字信號的頻譜分析介紹
一儿礼、傅里葉變換
傅里葉變換是信號領域溝通時域和頻域的橋梁乍钻,在頻域里可以更方便的進行一些分析肛循。傅里葉主要針對的是平穩(wěn)信號的頻率特性分析蛛株,簡單說就是具有一定周期性的信號,因為傅里葉變換采取的是有限取樣的方式育拨,所以對于取樣長度和取樣對象有著一定的要求。
二欢摄、基于Python的頻譜分析
# _*_ coding:utf-8 _*_
import numpy as np #導入一個數(shù)據(jù)處理的模塊
import pylab as pl #導入一個繪圖模塊熬丧,matplotlib下的模塊
sampling_rate = 8000 ##取樣頻率
fft_size ?=512 ? #FFT處理的取樣長度
t = np.arange(0,1.1,1.0/sampling_rate)
#np.arange(起點,終點怀挠,間隔)產(chǎn)生1s長的取樣時間
x = np.sin(2*np.pi*156.25*t)+2*np.sin(2*np.pi*234.375*t)
#兩個正弦波疊加析蝴,156.25HZ和234.375HZ,因此如上面簡單
#的介紹FFT對于取樣時間有要求绿淋,
#N點FFT進行精確頻譜分析的要求是N個取樣點包含整數(shù)個
#取樣對象的波形闷畸。
#因此N點FFT能夠完美計算頻譜對取樣對象的要求
#是n*Fs/N(n*采樣頻率/FFT長度),
#因此對8KHZ和512點而言吞滞,
#完美采樣對象的周期最小要求是8000/512=15.625HZ,
#所以156.25的n為10,234.375的n為15佑菩。
xs = x[:fft_size]# 從波形數(shù)據(jù)中取樣fft_size個點進行運算
xf = np.fft.rfft(xs)/fft_size # 利用np.fft.rfft()進行FFT計算,rfft()是為了更方便
#對實數(shù)信號進行變換裁赠,由公式可知/fft_size為了正確顯示波形能量
# rfft函數(shù)的返回值是N/2+1個復數(shù)殿漠,分別表示從0(Hz)
#到sampling_rate/2(Hz)的分。
#于是可以通過下面的np.linspace計算出返回值中每個下標對應的真正的頻率:
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))
#最后我們計算每個頻率分量的幅值佩捞,并通過 20*np.log10()?
#將其轉換為以db單位的值绞幌。為了防止0幅值的成分造成log10無法計算,
#我們調(diào)用np.clip對xf的幅值進行上下限處理
pl.figure(figsize=(8,4))
pl.subplot(211)
pl.plot(t[:fft_size], xs)
pl.xlabel(u"時間(秒)")
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()
#繪圖顯示結果
現(xiàn)在來看看頻譜泄露一忱,將采樣對象的頻率改變
x = np.sin(2*np.pi*100*t)+2*np.sin(2*np.pi*234.375*t)
我們明顯看出莲蜘,第一個對象的頻譜分析出現(xiàn)“泄露”,能量分散到其他頻率上帘营,
沒法準確計算到計算對象的頻譜特性票渠。
窗函數(shù)
上面我們可以看出可以通過加“窗”函數(shù)的方法來處理,盡量保證FFT長度內(nèi)
的取樣對象是對稱的仪吧。
import pylab as pl
import scipy.signal as signal
pl.figure(figsize=(8,3))
pl.plot(signal.hann(512))#漢明窗函數(shù)?
pl.show()