最近有朋友問我，他們被要求在單片機裡實現一個濾波器，參數等要求此處略去不表，他很憂愁怎麼在單片機實現濾波器進行濾波操作。
這個朋友真不是我自己
我告訴他，只需要Matlab就可，或者說只要你能確定濾波器抽頭的係數就可以了。為什麼濾波器只需要抽頭係數就可以了，這一切都要從濾波器的結構講起；為了方便講解，我們選用結構簡單的FIR濾波器講解。
本文通過講解FIR的系統結構，進而講解數字濾波器濾波原理和實現方法。文章涉及部分數字信號處理內容，實現方法設計Python語言和Matlab的使用以及部分線性代數的基礎知識。
查看FIR濾波器的系統結構
打開Matlab，在命令窗口輸入fdatool(新版MATLAB為filterDesigner)回車，Matlab的FIR相關基礎操作可以看（學習Verilog）6. FIR IP核的基礎功能使用總結
選擇FIR窗函數，72階，漢寧窗，fs:fc = 10:1
然後跟著下圖的操作順序來
按照操作來，接下來需要等待一段時間
等待一會後會打開Simulink，Simulink會出現下圖
Filter的實現框圖，點擊它
點擊Filter，會出現下圖
FIR濾波器實現框圖
各位電腦實際操作的時候可以縮小看，整個濾波器結構其實就是這樣子的。外部數據從上進入，然後經過延時到達下一級，我們稱每一級為抽頭；同一時刻下，抽頭的數據會跟抽頭係數相乘（圖中三角形），然後所有抽頭相乘的結構相加就是最後的結果。
為什麼這樣子就能實現濾波功能？
這裡首先告訴大家一個結論，所有濾波器實現的原理，無非是延時加權求和，具體的理論細節將在後續講解，FIR濾波器公式如下：
公式看著很簡單，也對應上了上述的系統結構，加權（係數相乘），延時，求和。不過，濾波器設計過程中，難點還是在於h(n)的設計實現，好在matlab幫助我們計算出了係數。
導出FIR濾波器係數
FIR濾波器實現過程中，最麻煩的係數已經獲取掉了，接下來我們就需要把係數導出，導出過程如下：
1. 直接導出數據到文本文件
導出數據
保存為txt文件
2. 導出C語言的頭文件，直接用在C/C++上
生成C 頭文件
由於家裡的電腦好久沒寫C了，Visual Studio2015不知道出什麼毛病，寫C一直報錯；Dev-C++也出現問題，一氣之下就用Python3寫完了後續過程。為了Python3方便調用，這裡我選擇了導出方法1。
3. 實現代碼講解
實驗所用的環境為Python3，需要預先安裝Numpy，Matplotlib。
參數設定
<code># FIR抽頭階數LEN = 73# FIR濾波器信號數據延時的數組
fir_data = np.zeros(LEN)# 製作一個信號源，對比效果# 採樣頻率
fs = 10000/<code>
首先設置FIR的階數，可能有人會有疑問，上面的圖片中FIR的階數是72，這裡怎麼是73；由於結構問題，FIR的數據進入的第一級只進行了延時，沒有係數相乘，或者說係數相乘0。而FIR的係數保持對稱關係，
為了保持這種關係，FIR係數最後會添加一個0，變成73個係數；但其實這裡可以去掉首尾的兩個0，係數就成了71個，不過這裡無傷大雅，我便沒有去除。
接著創建一個LEN長的一維零數組。設定採樣率為1000。
驗證FIR結構
<code>t = np.linspace(0, 1, fs)
f = np.cos(2*np.pi*10*t) + 2 * np.cos(2*np.pi*3000*t)# 獲取到濾波器，fs/fc = 10：1
coef = set_coef('coef.txt')# 分配一個fs長的零數組
ret = np.zeros(fs)# 獲取FIR延時加權求和的結果# 可以選擇FIR_INST函數，操作流程直觀# 可以選擇FIR_INST_ARRAY函數，選用了數組點乘切片等方式，去掉循環時間應該更快for i in range(fs):
ret[i] = fir_inst_array(f[i],coef)
# ret[i] = fir_inst(f[i],coef)/<code>
創造一個採樣率10000下，振幅為1的10hz正弦波和振幅為2的3000hz正弦波疊加的信號源。接著開始循環，每次循環過程中讓信號源數據經過FIR操作後，保存結果。
FIR結構實現
<code>def fir_inst(data, coef) :
# '''
# :param data: 原始信號 Orignal Signal
# :param coef: FIR抽頭係數數組
# :return: 濾波之後的信號
# '''
ret = 0
# 加權延時求和操作
# 第一個抽頭和最後的抽頭係數一定是0，所以循環只從LEN-1到1
# 第一個抽頭的操作有點不同，它接收進來的data
for i in range(LEN-1, 0, -1):
ret = ret + fir_data[i] * coef[i]
fir_data[i] = fir_data[i-1]
fir_data[0] = data
fir_sum_data[0] = fir_data[0] * coef[0]
return ret/<code>
這個函數的實現很直觀，可以照著這個函數寫出C/C++版本。
循環中，延時的FIR數據與抽頭係數相乘，然後結果進行累加，整個循環也完成了FIR整體數據延時一個循環週期的要求。最後將新進入的數據保存至FIR數據數組中。
FIR結構實現（數組實現方法）
<code>def fir_inst_array(data, coef):
# '''
# :param data: 同上 使用數組乘法等規則簡化和加快運算過程
# :param coef:
# :return:
# '''
# FIR_DATA數組使用切片索引等方式實現數組左移並在0位加上新的數據
fir_data[1:] = fir_data[0:LEN-1]
fir_data[0] = data
# 使用數組的點乘獲取加權求和的結果
ret = np.dot(fir_data, coef)
return ret/<code>
既然使用了Numpy，肯定有更適合數組運算的方法。在這個函數里面，FIR數據的延時使用切片的方式，使整體數據左移，然後將新加入的數據保存在0位。
FIR數據與抽頭係數相乘直接使用了Numpy的點乘，兩個數組點乘可以直接得到加權求和的結果。
結果對比
藍色部分為信號源，紅色的線就是濾波後的信號。結果表明，濾波器效果已經實現。不過，整個濾波器實現中，我們只是實現了濾波器結構，濾波器抽頭係數是由軟件確定。數字信號處理，研究的就是系統衝激響應序列h(n)的確定以及實現。
歡迎對代碼有更好的改進建議，Python用的不熟吐槽不要太狠。。。
對於IIR濾波器，其實也可以使用這種方法獲取IIR濾波器系統結構，根據參數寫出對應的IIR濾波器實現。
閱讀更多 BUG記錄 的文章
