引言

數字濾波器在語音與圖像處理、模式識別、雷達信號處理、頻譜分析等應用中都具有重要作用。它能避免模擬濾波器所無法克服的溫漂和噪聲等問題， 同時比模擬濾波器精度高、穩定性好、體積小、更加靈活， 因而得到廣泛應用。在聲波測井中， 通常需要對信號進行精確的濾波，并且對濾波器具有嚴格的實時性要求。本文利用輔助Matlab設計工具， 設計了一種基于FPGA且可滿足測井需求的高階快速數字濾波器。

1 線性相位FIR濾波器結構

數字濾波器的種類很多， 分類的方法也不盡相同。從數字濾波器的單位沖擊響應來看， 數字濾波器分為有限沖激響應數字濾波器（FIR） 和無限沖激響應數字濾波器（IIR）。相對于IIR濾波器，FIR濾波器可以進行準確的線性相位設計， 而且其結構具有穩定的量化濾波器系數。針對于聲波測井處理具有線性相位要求的聲波信號， FIR濾波器是首選。

在時域中， FIR濾波器的輸入輸出過程是一個輸入信號與單位脈沖響應進行線性卷積的過程， 其差分方程表達式為：

其中， y （n） 為濾波輸出， x （n） 為采樣數據， h （n） 為濾波器抽頭系數。其結構如圖1 （a）所示， 圖中， N-1階FIR濾波器要用N個系數描述， 通常， 需要N個乘法器和N-1個兩輸入加法器才能實現。不難發現， 乘法器的系數正好是傳輸函數的系數， 因此， 該結構且稱為直接型結構。

對于系數對稱的FIR線性相位濾波器， 可將式（1） 寫成如下形式：

系數對稱的改進型FIR濾波器的結構如圖1（b） 所示。該結構把系數對稱（相同或相反） 的抽頭合并之后再作乘法， 這樣可使乘法器數量降為原有的一半， 但也增加了額外的加法器。

圖1 FIR濾波器結構

2 設計方法與指標

FDATool是Matlab信號處理工具箱里專用的濾波器設計分析工具， 該工具的主要作用是按照設計指標提取濾波器系數。用FDATool設計數字濾波器的關鍵在于濾波器類型、窗函數、濾波器階數、截止頻率等參數的選擇。其中窗函數用于決定阻帶衰減和過渡帶帶寬， 常用的窗函數有矩形窗、漢寧窗、海明窗和布萊克曼窗。矩形窗和漢寧窗阻帶衰減較小， 而布萊克曼窗過渡帶較大， 相對來說， 海明窗更符合設計要求， 它的最小阻帶可以達到54.5dB， 歸一化過渡帶帶寬為3.11π/M （濾波器階數N=2M+1）。針對聲波測井信號， 設計時應進行表1所列的參數設置。

表1 濾波器參數選擇

圖2所示是該濾波器的幅頻和相頻響應曲線，該曲線在通帶內保持線性相位， 阻帶衰減大于52dB， 過渡帶帶寬為1.65kHz。抽頭系數可以在工具箱中量化為定點整型數據， 以便在FPGA實現階數為127的濾波器， 該濾波器一共有128個系數。對于階數較大的濾波器來說， 其量化對阻帶衰減和過渡帶的影響極小。

圖2 濾波器幅頻和相頻響應特性曲線

3 基于FPGA的濾波器設計

用FPGA設計FIR濾波器的關鍵在于如何處理占用大量資源的乘法單元。分布式算法（DA） 的提出可將乘法運算轉換為移位相加運算， 從而節約硬件資源。若令Hk為濾波器系數， xk （n） 為n時刻的采樣輸入， y （n） 為n時刻的系統響應， 那么， 式（1） 就可以等效于下式：

把數據的源數據格式規定為2的補碼形式，則有：

式中， xkb （n） 為二進制數， 可取值為0或1；xk0 （n） 為符號位， 為1表示數據為負， 為0表示數據為正。因此， 將（4） 式代入（3） 式可得：

式（5） 的形式被稱為分布式算法。可以看出， 方括號內表示輸入變量的一個數據位和所有濾波器抽頭系數H0～HN的每一位進行“與” 運算并求和。而指數部分則說明了求和結果的位權，整數乘以2b就是左移b位， 對此可以通過硬件連線來實現， 而不占用邏輯資源。這樣就可以通過建立查找表來實現方括號中的運算。查找表可用所有輸入變量的同一位進行尋址， 這便是基于查找表的分布式算法（LUT-DA）。

LUT-DA算法的查找表大小為B·2N bits， 其中B為輸入數據的位寬， N為濾波器階數。隨著濾波器階數的增加， 查找表大小是2的指數增長；當B為16， N為128時， 查找表的大小已經不可想象。故將查找表分割成多個子表， 可以有效解決這個問題， 這也衍生了比較有效的串行LUT-DA算法和并行LUT-DA算法， 但兩者都有不足的地方。對于串行結構， 要完成一次輸出， 需要大于B的多個時鐘周期； 而對于并行結構， 雖然可以一個時鐘周期完成一次輸出， 但需要復制B個完全相同的LUT表， 而這會增加硬件資源的開銷。

為了兼顧速度和面積， 本文設計了一種基于DA算法原理的CSD-DA算法。首先， 將系數式（3） 中的固定系數Hk按2的冪展開后可得：

然后交換移位和累加順序， 則可得到下式：

式中， Hkb是值為0或者1的權重系數； Sk為1表示Hk為正， 為-1則表示Hk為負； s′kb的值可取0、-1或者1。經過（4） 式的展開， 乘法運算將被全部轉換為移位相加運算， 其中權重為0的部分可以剔除而不進行計算。為了更進一步減少Hkb陣列中的非零項， 可將Hk編碼為CSD碼， 即從二進制編碼最低的有效位開始， 用10···01來取代所有大于或等于2的1序列， 1表示該位為-1。由于CSD表示其中任何相鄰的兩位中， 必包含一個0，故1的數量最多不會超過N/2。平均來說， CSD表示其中大約有1/3的位為非零值， 這比補碼表示少大約1/3的非零位。假設h= （15） 10= （01111） 2，y=hx=x （23+22+21+20）， 而如果將（15） 10編碼為（10001） csd， 那么， Y=x· （24-20）。采用二進制編碼方式， 將用到3個加法器， 而用CSD編碼， 則只用了一個減法器， 可見， CSD編碼可以從本質上減少硬件資源開銷。經過CSD編碼優化后， s′kb非零值的個數會遠小于Hkb的非零值個數。

對于線性相位系數對稱的FIR濾波器， 為了減少乘法單元， 可選擇圖3所示的結構。由于所有的乘法運算都可轉化為大量的加法和減法運算， 故將導致關鍵路徑過長， 系統運行速度較低。而加入流水線寄存器， 則可減少關鍵路徑長度， 從而提高系統的最大工作頻率。在b為定值時， s′kb的非零值個數存在不確定性， 故在進行流水線設計的時候， 可根據s′kb進行靈活的分割， 路徑越長， 加入的流水線寄存器越多。為了防止中間結果的溢出， 寄存器的位寬要有冗余設計， 對于有符號的數， 其位寬取值為M+log2N-1， M為上級累加器位寬， N為濾波器階數。

圖3 流水線CSD-DA算法局部結構

從圖3的流水優化CSD-DA算法結構可見， 所有乘法都會轉換為移位加法， 移位運算可用硬件連線實現， 整個結構經過了合理的流水線分割。

表2所列是不同結構濾波器的綜合結果。其中并行結構是最差的一種結構， 它占用資源多、速度慢； 串行LUT-DA結構， 雖然占用資源少，最大工作頻率高， 但畢竟是串行結構， 不能在一個時鐘節拍完成對一個采樣點的濾波運算； 而流水線CSD-DA結構則無論在速度， 還是面積上，都具有比較明顯的優勢。如果工作時鐘為75MHz， 那么， 一個時鐘節拍便可以完成一次輸出， 處理330個采樣點的單通道信號僅需4.4μs，可以滿足測井實時性要求。

表2 濾波器的綜合結果

4 結果分析

為了驗證濾波器的功能是否正確， 可對本設計在Modelsim中進行仿真。若原始波形為帶噪聲的聲波信號， 那么， 其濾波的結果如圖4所示。

圖4 濾波器在Modelsim中的仿真結果

圖5所示是濾波器在Matlab中的仿真結果， 可以看出， Modelsim和Matlab的仿真結果一致。在頻域上， 對比圖5 （a） 和圖5 （b） 可以看到， 其濾波后的波形只保留了5kHz~18kHz的頻譜部分， 這說明流水線CSD-DA結構的數字濾波器設計正確。

圖5 濾波器在Matlab中的仿真結果

5 結束語

本文詳細講述了通過Matlab工具設計FIR線性相位濾波器的方法， 并針對聲波信號設計了優于傳統結構的流水線CSD-DA結構， 該結構具有較明顯的速度和面積優勢。文中也通過仿真實驗證實了設計的合理性和正確性。但值得指出的是， 該結構只適合固定濾波器系數的場合， 而如果要進行修改， 則需要重新對系數進行CSD編碼和流水線分割。