帶通濾波器設計

相關文檔： SciPy-數值計算庫

用scipy.signal庫設計帶通濾波器，并用頻率掃描波測試其頻率響應。

# -*- coding: utf-8 -*-
import scipy.signal as signal
import pylab as pl
import numpy as np

#取樣頻率為8kHz
sampling_rate = 8000.0

# 設計一個帶通濾波器：
# 通帶為0.2*4000 - 0.5*4000
# 阻帶為<0.1*4000, >0.6*4000
# 通帶增益的最大衰減值為2dB
# 阻帶的最小衰減值為40dB
b, a = signal.iirdesign([0.2, 0.5], [0.1, 0.6], 2, 40)

# 使用freq計算濾波器的頻率響應
w, h = signal.freqz(b, a)

# 計算增益
power = 20*np.log10(np.clip(np.abs(h), 1e-8, 1e100))

# 繪制增益
pl.subplot(211)
pl.plot(w/np.pi*sampling_rate/2, power)
pl.title(u"freqz計算的濾波器頻譜")
pl.ylim(-100,20)
pl.ylabel(u"增益(dB)")

# 產生2秒鐘的取樣頻率為sampling_rate Hz的頻率掃描信號
# 開始頻率為0， 結束頻率為sampling_rate/2
t = np.arange(0, 2, 1/sampling_rate)
sweep = signal.chirp(t, f0=0, t1 = 2, f1=sampling_rate/2)

# 將頻率掃描信號進行濾波
out = signal.lfilter(b, a, sweep)

# 將波形轉換為能量
out = 20*np.log10(np.abs(out))

# 找到所有局部最大值的下標
index = np.where(np.logical_and(out[1:-1] > out[:-2], out[1:-1] > out[2:]))[0] + 1

# 繪制濾波之后的波形的增益
pl.subplot(212)
pl.plot(t[index]/2.0*4000, out[index] )
pl.title(u"頻率掃描波測量的濾波器頻譜")
pl.ylim(-100,20)
pl.ylabel(u"增益(dB)")
pl.xlabel(u"頻率(Hz)")

pl.subplots_adjust(hspace=0.3)
pl.show()