最近要把巴特沃斯濾波算法移植到嵌入式端，之前一直是在服務(wù)器端用python搞的递宅，簡(jiǎn)單粗暴娘香，幾行代碼就解決問(wèn)題了苍狰，如下：
import numpy as np
from scipy import signal
def data_filt(array):
array = array-np.mean(array, axis=1, keepdims=True)
b, a = signal.butter(4, [0.5, 30], 'bandpass',fs=Fs)
print('b = ', b)
print('a = ', a)
data = signal.filtfilt(b, a, array, axis=1)
return data
b = [ 0.00842863 0. -0.03371451 0. 0.05057177 0.
-0.03371451 0. 0.00842863]
a = [ 1. -6.05023914 16.06728061 -24.56190002 23.72091821
-14.84696022 5.88026366 -1.34582827 0.13646518]
采用scipy.signal里面的butter函數(shù)，先計(jì)算濾波系數(shù)b,a,然后利用signal.filtfilt進(jìn)行iir濾波烘绽。
效果展示如下淋昭，輸入一個(gè)采樣率為250Hz的信號(hào)，進(jìn)行[0.5 30]Hz范圍內(nèi)的帶通濾波安接。
可以看出翔忽，對(duì)高頻分量和低頻分量都進(jìn)行了有效濾除，效果還不錯(cuò)盏檐。
總結(jié)下總共兩步：
- 計(jì)算butter系數(shù)b歇式，a
- 將b，a輸入到iir濾波器胡野，進(jìn)行濾波
其中iir濾波的公式如下：
現(xiàn)在要做的是材失，把這兩步用c實(shí)現(xiàn)。
1. 計(jì)算butter系數(shù)
找了很多教程硫豆，最后發(fā)現(xiàn)了一個(gè)神奇的網(wǎng)站龙巨，https://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html，可以直接幫助計(jì)算參數(shù)b熊响，a,且給出濾波代碼旨别。
如上圖，把參數(shù)填好后耘眨，直接點(diǎn)擊網(wǎng)頁(yè)下方的提交昼榛，濾波系數(shù)就生成了,如下：
其中參數(shù)b,a就隱藏在Recurrence relation里面境肾。根據(jù)公式
y[n] = ( 1 * x[n- 8])
+ ( 0 * x[n- 7])
+ ( -4 * x[n- 6])
+ ( 0 * x[n- 5])
+ ( 6 * x[n- 4])
+ ( 0 * x[n- 3])
+ ( -4 * x[n- 2])
+ ( 0 * x[n- 1])
+ ( 1 * x[n- 0])
+ ( -0.1364651827 * y[n- 8])
+ ( 1.3458282682 * y[n- 7])
+ ( -5.8802636583 * y[n- 6])
+ ( 14.8469602185 * y[n- 5])
+ (-23.7209182114 * y[n- 4])
+ ( 24.5619000249 * y[n- 3])
+ (-16.0672806050 * y[n- 2])
+ ( 6.0502391423 * y[n- 1])
這里面的對(duì)應(yīng)的系數(shù)就是
。即
奥喻。同理可得
跟上面python的結(jié)果對(duì)比下偶宫，發(fā)現(xiàn)a是一致的，但是b不一致环鲤，原因是因?yàn)樗麄冎g差了一個(gè)gain纯趋。就是Results那副圖里面的 gain at centre: mag = 1.185137516e+02, b的每個(gè)值都除以1.185137516e+02，就和python的結(jié)果一致了冷离。
2. 將參數(shù)b吵冒，a用于iir濾波
主要就是對(duì)公式的實(shí)現(xiàn)。
比較牛逼的點(diǎn)就是西剥，上面推薦的那個(gè)網(wǎng)站把這個(gè)也給做好了痹栖，我們基本上只需要復(fù)制。
這是網(wǎng)站生成的Ansi C的代碼瞭空。
#define NZEROS 8
#define NPOLES 8
#define GAIN 1.185137516e+02
static float xv[NZEROS+1], yv[NPOLES+1];
static void filterloop()
{ for (;;)
{ xv[0] = xv[1]; xv[1] = xv[2]; xv[2] = xv[3]; xv[3] = xv[4]; xv[4] = xv[5]; xv[5] = xv[6]; xv[6] = xv[7]; xv[7] = xv[8];
xv[8] = next input value / GAIN;
yv[0] = yv[1]; yv[1] = yv[2]; yv[2] = yv[3]; yv[3] = yv[4]; yv[4] = yv[5]; yv[5] = yv[6]; yv[6] = yv[7]; yv[7] = yv[8];
yv[8] = (xv[0] + xv[8]) - 4 * (xv[2] + xv[6]) + 6 * xv[4]
+ ( -0.1364651827 * yv[0]) + ( 1.3458282682 * yv[1])
+ ( -5.8802636583 * yv[2]) + ( 14.8469602180 * yv[3])
+ (-23.7209182110 * yv[4]) + ( 24.5619000250 * yv[5])
+ (-16.0672806050 * yv[6]) + ( 6.0502391423 * yv[7]);
next output value = yv[8];
}
}
利用這個(gè)代碼揪阿，對(duì)信號(hào)進(jìn)行濾波疗我，并同python的結(jié)果對(duì)比，如下：
發(fā)現(xiàn)效果基本上都有了溺健，但是還是差一點(diǎn)點(diǎn)麦牺，仔細(xì)看還是沒有python濾波的效果好。
那么問(wèn)題出在了那里呢鞭缭？
第一步計(jì)算參數(shù)b,a枕面，python和c的結(jié)果是一樣的。所以只能是出在了第二步：iir濾波缚去。
python用到的是signal.filtfilt()函數(shù)潮秘。我們看下它的源碼的核心部分：
原因找到了，在filtfilt函數(shù)里易结，進(jìn)行完Forward filter濾波——lfilter以后枕荞，把結(jié)果y做了翻轉(zhuǎn)，然后又進(jìn)行了一次lfilter濾波——即Backword濾波搞动。最后再把y翻轉(zhuǎn)回來(lái)躏精，進(jìn)行return。
相比較的話鹦肿，我們的代碼里只進(jìn)行了Forward濾波矗烛，沒有Backword。所以我們仿照這個(gè)方法箩溃，把生成的結(jié)果y進(jìn)行翻轉(zhuǎn)瞭吃，再做一次相同的濾波，最后再翻轉(zhuǎn)回來(lái)涣旨，作為最終的濾波結(jié)果歪架。對(duì)照?qǐng)D如下：
YES!! 這一次python和c的結(jié)果就一模一樣了。大功告成霹陡，完美和蚪！
