嵌入式开发日记（8）——用python实现FIR滤波（未完待续）
第一階段的方法是根據(jù)單位時間內(nèi)的加速度絕對值差值來判斷震顫程度，存在很多問題。因此設(shè)想采用更加高級的算法來加以改進。
這部分的主要工作有：? 1 學(xué)習(xí)數(shù)字信號處理的濾波算法，重點學(xué)習(xí)python下使用FIR濾波算法
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?2 利用python的數(shù)據(jù)分析工具對傳感器信號加以分析，利用好更多的信號數(shù)據(jù)
FIR濾波器很簡單，它實際上是一個全零點模型（MA滑動平均模型），濾波器系數(shù)只包含滑動平均的B，而自回歸的A=1。濾波過程就是X和B的卷積，再聯(lián)系濾波器的原理，容易知道，FIR濾波器的系數(shù)B就是濾波器的沖激響應(yīng)h。
因此它的設(shè)計就簡單了，只要設(shè)定濾波器的頻率響應(yīng)H，進行ifft后就得到?jīng)_激響應(yīng)h，把h直接作為濾波器系數(shù)B對原信號濾波即可。
經(jīng)查找，python下可以通過scipy.signal這個科學(xué)計算模塊實現(xiàn)對信號的處理濾波。python在科學(xué)計算領(lǐng)域有三個非常受歡迎庫，numpy、SciPy、matplotlib。numpy是一個高性能的多維數(shù)組的計算庫，SciPy是構(gòu)建在numpy的基礎(chǔ)之上的，它提供了許多的操作numpy的數(shù)組的函數(shù)。SciPy是一款方便、易于使用、專為科學(xué)和工程設(shè)計的python工具包，它包括了統(tǒng)計、優(yōu)化、整合以及線性代數(shù)模塊、傅里葉變換、信號和圖像圖例，常微分方差的求解等。
?
基礎(chǔ)知識：
傅立葉變換的形象講解：https://zhuanlan.zhihu.com/p/19763358
圖片說明（摘自上面鏈接）：
numpy教程：快速傅里葉變換模塊numpy.fft??https://blog.csdn.net/pipisorry/article/details/51050297
下面主要是對官方文檔的復(fù)現(xiàn)與翻譯：
巴特沃斯濾波器https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html：
scipy.signal.butter（N，Wn，btype ='低'，模擬=假，輸出='ba'，fs =無）
| N?：?int 過濾器的順序。 Wn?：?array_like 給出臨界頻率的標(biāo)量或長度為2的序列。對于巴特沃斯濾波器，這是增益下降到通帶的1 / sqrt（2）的點（“-3 dB點”）。 對于數(shù)字濾波器，Wn與fs的單位相同。默認(rèn)情況下，?fs為2個半周/采樣，因此這些從0到1歸一化，其中1是奈奎斯特頻率。（因此，Wn處于半周期/樣本中。） 對于模擬濾波器，Wn是角頻率（例如rad / s）。 btype?：?{'?lowpass?'，'?highpass?'，'bandpass'，'bandstop'}，可選 過濾器的類型。默認(rèn)為'低通'。 模擬?：?bool，可選 如果為True，則返回模擬濾波器，否則返回數(shù)字濾波器。 輸出?：?{'ba'，'zpk'，'sos'}，可選 輸出類型：分子/分母（'ba'），極點零（'zpk'）或二階段（'sos'）。默認(rèn)為'ba'。 fs?：?float，可選 數(shù)字系統(tǒng)的采樣頻率。 版本1.2.0中的新功能。 |
| b，a?：?ndarray，ndarray 分子（b）和分母（a）IIR濾波器的多項式。只有返回output='ba'。 z，p，k?：?ndarray，ndarray，float IIR濾波器傳遞函數(shù)的零點，極點和系統(tǒng)增益。只有返回output='zpk'。 sos?：?ndarray IIR濾波器的二階截面表示。只有返回output=='sos'。 |
例子：左圖高通，右圖低通
#設(shè)計模擬濾波器并繪制其頻率響應(yīng)，顯示關(guān)鍵點： from scipy import signal import matplotlib.pyplot as pltb, a = signal.butter(4, 100, 'low', analog=True) w, h = signal.freqs(b, a) plt.semilogx(w, 20 * np.log10(abs(h))) plt.title('Butterworth filter frequency response') plt.xlabel('Frequency [radians / second]') plt.ylabel('Amplitude [dB]') plt.margins(0, 0.1) plt.grid(which='both', axis='both') plt.axvline(100, color='green') # cutoff frequency plt.show()#產(chǎn)生由10 Hz和20 Hz組成的信號，以1 kHz采樣 t = np.linspace(0, 1, 1000, False) # 1 second sig = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*20*t) fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True) ax1.plot(t, sig,color='red') ax1.set_title('10 Hz and 20 Hz sinusoids') ax1.axis([0, 1, -2, 2])#設(shè)計15 Hz的數(shù)字高通濾波器以消除10 Hz音調(diào)，并將其應(yīng)用于信號。（建議在過濾時使用二階段格式，以避免傳遞函數(shù)（ba）格式的數(shù)值錯誤）： sos = signal.butter(10, 15, 'hp', fs=1000, output='sos') filtered = signal.sosfilt(sos, sig) ax2.plot(t, filtered,color='green') ax2.set_title('After 15 Hz high-pass filter') ax2.axis([0, 1, -2, 2]) ax2.set_xlabel('Time [seconds]') plt.tight_layout() plt.show()?
?使用IIR或FIR濾波器沿一維過濾數(shù)據(jù):https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lfilter.html
scipy.signal.lfilter?
scipy.signal.lfilter（b，a，x，軸= -1，zi =無）
使用IIR或FIR濾波器沿一維過濾數(shù)據(jù)。
使用數(shù)字濾波器過濾數(shù)據(jù)序列x。這適用于許多基本數(shù)據(jù)類型（包括對象類型）。濾波器是標(biāo)準(zhǔn)差分方程的直接形式II轉(zhuǎn)置實現(xiàn)（見注釋）。
| b?：?array_like 1-D序列中的分子系數(shù)向量。 a?：?array_like 1-D序列中的分母系數(shù)向量。如果a[0]?不是1，則a和b都?xì)w一化a[0]。 x?：?array_like N維輸入數(shù)組。 axis?：?int，可選 輸入數(shù)據(jù)陣列的軸，沿著該軸應(yīng)用線性濾波器。過濾器沿該軸應(yīng)用于每個子陣列。默認(rèn)值為-1。 zi?：?array_like，可選 過濾器延遲的初始條件。它是長度的向量（或用于N維輸入的向量的數(shù)組）?。如果zi為None或未給出，則假定初始休息。有關(guān)更多信息，請參閱max(len(a),?len(b))?-?1lfiltic |
| y?：?數(shù)組 數(shù)字濾波器的輸出。 zf?：?array，可選 如果zi為None，則不返回，否則，zf保存最終的濾波器延遲值。 |
例子：右圖為只有filt:
from scipy import signal import matplotlib.pyplot as plt import numpy as np#生成要過濾的噪聲信號 t = np.linspace(-1, 1, 201) x = (np.sin(2*np.pi*0.75*t*(1-t) + 2.1) +0.1*np.sin(2*np.pi*1.25*t + 1) +0.18*np.cos(2*np.pi*3.85*t)) xn = x + np.random.randn(len(t)) * 0.08#創(chuàng)建order3低通道butterworth過濾器： b, a = signal.butter(3, 0.05) #將過濾器應(yīng)用于xn。使用lfilter_zi選擇過濾器的初始條件： #zi = signal.lfilter_zi(b, a) #z, _ = signal.lfilter(b, a, xn, zi=zi*xn[0]) #再次應(yīng)用過濾器，以與filtfilt相同的順序過濾結(jié)果： #z2, _ = signal.lfilter(b, a, z, zi=zi*z[0]) #使用filtfilt應(yīng)用過濾器： y = signal.filtfilt(b, a, xn) #繪制原始信號和各種過濾版本： plt.figure plt.plot(t, xn, 'b', alpha=0.75) plt.plot(t, z, 'r--', t, z2, 'r', t, y, 'k') plt.legend(('noisy signal', 'lfilter, once', 'lfilter, twice','filtfilt'), loc='best') plt.grid(True) plt.show()lfilter還是filtfilt的區(qū)別：https://dsp.stackexchange.com/questions/19084/applying-filter-in-scipy-signal-use-lfilter-or-filtfilt
-
filtfilt是零相位濾波，在濾波時不會移位信號。由于在所有頻率下相位均為零，因此它也是線性相位。及時過濾需要您預(yù)測未來，因此不能在“在線”實際應(yīng)用中使用，僅用于信號記錄的離線處理。
-
lfilter是一種因果時間過濾，類似于現(xiàn)實生活中的電子過濾器。它不能是零相位。它可以是線性相位（對稱FIR），但通常不是。通常它會在不同的頻率上增加不同的延遲量。
一個例子和圖像應(yīng)該使它顯而易見。雖然濾波器的頻率響應(yīng)幅度相同（左上角和右上角），但零相位低通與原始信號對齊，只是沒有高頻內(nèi)容，而最小相位濾波以因果方式延遲信號：
from __future__ import division, print_function import numpy as np from numpy.random import randn from numpy.fft import rfft from scipy import signal import matplotlib.pyplot as pltb, a = signal.butter(4, 0.03, analog=False)# Show that frequency response is the same impulse = np.zeros(1000) impulse[500] = 1# Applies filter forward and backward in time imp_ff = signal.filtfilt(b, a, impulse)# Applies filter forward in time twice (for same frequency response) imp_lf = signal.lfilter(b, a, signal.lfilter(b, a, impulse))plt.subplot(2, 2, 1) plt.semilogx(20*np.log10(np.abs(rfft(imp_lf)))) plt.ylim(-100, 20) plt.grid(True, which='both') plt.title('lfilter')plt.subplot(2, 2, 2) plt.semilogx(20*np.log10(np.abs(rfft(imp_ff)))) plt.ylim(-100, 20) plt.grid(True, which='both') plt.title('filtfilt')sig = np.cumsum(randn(800)) # Brownian noise sig_ff = signal.filtfilt(b, a, sig) sig_lf = signal.lfilter(b, a, signal.lfilter(b, a, sig)) plt.subplot(2, 1, 2) plt.plot(sig, color='silver', label='Original') plt.plot(sig_ff, color='#3465a4', label='filtfilt') plt.plot(sig_lf, color='#cc0000', label='lfilter') plt.grid(True, which='both') plt.legend(loc="best")scipy.signal.filtfilt模塊：https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.filtfilt.html#scipy.signal.filtfilt
scipy.signal.filtfilt（b，a，x，axis = -1，padtype ='odd'，padlen = None，method ='pad'，irlen = None?）[source]
向數(shù)字濾波器前后應(yīng)用信號。
此功能應(yīng)用線性數(shù)字濾波器兩次，一次向前，一次向后。組合濾波器的零相位和濾波器順序是原始濾波器的兩倍。
該功能提供了處理信號邊緣的選項。
| b?：?（N，）array_like 濾波器的分子系數(shù)向量。 a?：?（N，）array_like 濾波器的分母系數(shù)向量。如果a[0]?不是1，則a和b都?xì)w一化a[0]。 x?：?array_like 要過濾的數(shù)據(jù)數(shù)組。 axis?：?int，可選 應(yīng)用濾波器的x軸。默認(rèn)值為-1。 padtype?：?str或None，可選 必須是'奇數(shù)'，'偶數(shù)'，'常數(shù)'或無。這決定了用于應(yīng)用濾波器的填充信號的擴展類型。如果padtype為None，則不使用填充。默認(rèn)值為“奇數(shù)”。 padlen?：?int或None，可選 在應(yīng)用過濾器之前在軸的兩端?延伸x的元素數(shù)。該值必須小于?。?意味著沒有填充。默認(rèn)值為。x.shape[axis]?-?1padlen=03?*?max(len(a),?len(b)) 方法?：?str，可選 確定處理信號邊緣的方法，“pad”或“gust”。當(dāng)方法是“pad”時，信號被填充;?填充的類型由padtype和padlen決定，irlen?被忽略。當(dāng)方法是“陣風(fēng)”時，使用Gustafsson的方法，并忽略padtype和padlen。 irlen?：?int或None，可選 當(dāng)方法是“陣風(fēng)”時，irlen指定濾波器的脈沖響應(yīng)的長度。如果irlen為None，則不會忽略脈沖響應(yīng)的任何部分。對于長信號，指定?irlen可以顯著提高濾波器的性能。 |
| y?：?ndarray 濾波后的輸出與x的形狀相同。 |
?
參考資料：
SciPy完整的教程：https://docs.scipy.org/doc/scipy/reference/index.html。
numpy的使用：https://blog.csdn.net/cxmscb/article/details/54583415
SciPy上手：?https://blog.csdn.net/q583501947/article/details/76735870
利用Python scipy.signal.filtfilt() 實現(xiàn)信號濾波：https://blog.csdn.net/weixin_37996604/article/details/82864680
python實現(xiàn)低通濾波器：https://blog.csdn.net/kkkxiong1/article/details/84941992
??https://blog.csdn.net/daiyinger/article/details/48289587
python實現(xiàn)高通濾波器：https://blog.csdn.net/thoughts_storms/article/details/32318921
python實現(xiàn)帶通濾波器：https://blog.csdn.net/weixin_39739342/article/details/84948308
python實現(xiàn)數(shù)字濾波器FIR，IIR：https://blog.csdn.net/qq_21970857/article/details/46685923
數(shù)據(jù)分析scipy常用方法?：https://www.cnblogs.com/why957/p/9308916.html
總結(jié)
以上是生活随笔為你收集整理的嵌入式开发日记（8）——用python实现FIR滤波（未完待续）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 【论文笔记】Enabling techn
- 下一篇： 探索技术论坛 GhostXP SP2 精