00 載入python擴(kuò)展庫

import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

01 定義信號

t=np.arange(0,10,1/256)
xt=1.5*np.sin(2*np.pi*40*t)+3*np.sin(2*np.pi*60*t)+2*np.sin(2*np.pi*80*t)

#40Hz，60Hz，80Hz的三個諧波疊加

02 傅里葉變換

xf=np.fft.fft(xt) #傅里葉變換
af=np.abs(xf)/len(xf) #求復(fù)數(shù)的模
afn=np.short(af.size/2) #整型
f=np.arange(0,256/2,256/len(t)) #定義頻率點(diǎn)
plt.subplot(211) #子圖1
plt.plot(t,xt) #時域信號
plt.subplot(212) #子圖2
plt.plot(f,2*af[:afn]) #頻域信號
plt.xlabel('Hz')

03 設(shè)計(jì)butterworth低通濾波器

N, Wn = sig.buttord(45, 55, 3, 30,fs=256) #數(shù)字濾波器
b, a = sig.butter(N, Wn, btype='lowpass',fs=256) #數(shù)字濾波器
w, h = sig.freqz(b,a,fs=256)
plt.plot(w, 20 * np.log10(abs(h)))
plt.xlabel('fre [Hz]')
plt.ylabel('response [dB]')

04 使用低通濾波器進(jìn)行濾波

x2t=sig.lfilter(b,a,xt)
x2f=np.fft.fft(x2t) #傅里葉變換
a2f=np.abs(x2f)/len(x2f) #求復(fù)數(shù)的模
a2fn=np.short(a2f.size/2) #整型
f=np.arange(0,256/2,256/len(t)) #定義頻率點(diǎn)
plt.subplot(211) #子圖1
plt.plot(t,x2t)
plt.subplot(212) #子圖2
plt.plot(f,2*a2f[:a2fn]) #頻域信號

可見，經(jīng)過butter低通濾波，60Hz，80Hz的頻率消失了；

05 設(shè)計(jì)butterworth 高通濾波器

N, Wn = sig.buttord(55, 45, 3, 30,fs=256) #數(shù)字濾波器
b, a = sig.butter(N, Wn, btype='highpass',fs=256) #數(shù)字濾波器
w, h = sig.freqz(b,a,fs=256)
plt.plot(w, 20 * np.log10(abs(h)))
plt.xlabel('fre [Hz]')
plt.ylabel('response [dB]')

06 使用高通濾波器進(jìn)行濾波

可見，經(jīng)過butter高通濾波，40Hz的頻率消失了；

07 設(shè)計(jì)butterworth 帶通濾波器

N, Wn = sig.buttord([45,75], [40,80], 3, 30,fs=256) #數(shù)字濾波器
b, a = sig.butter(N, Wn, btype='bandpass',fs=256) #數(shù)字濾波器
w, h = sig.freqz(b,a,fs=256)
plt.plot(w, 20 * np.log10(abs(h)))
plt.xlabel('fre [Hz]')
plt.ylabel('response [dB]')

08 使用帶通濾波器進(jìn)行濾波

可見，經(jīng)過butter濾波，只有60Hz留下了；