A、C計權數字濾波器的設計實現

8 人贊了文章

最近在做聲學實驗，需要將測得的聲壓數據做三分之一倍頻程（A計權），其實直接三分之一倍頻程是很好實現的（後續更新會講到實現方法），但是如何進行A計權呢？有的辦法是在做完倍頻程後直接在響應頻帶加減相應中心頻率的修正分貝值，但是如果我想在做完頻域計權後進行時域計權處理（如求任何瞬時時間上的 A 計權和時間計權聲級公式1），這種方法是沒辦法直接將上述結果通過IFFT變換回時域信號繼續處理的。

公式1 瞬時時間上的 A 計權和時間計權聲級

在解決這個問題前，需要看看A計權的頻域圖形，如圖1所示，不同計權的頻譜其實都可以看成一個加權濾波器，只要直接將時域信號通過該濾波器進行濾波，就可以得到計權後的時域信號，這不僅可以在時域上進行後處理，還可以對時域信號進行三分之一倍頻程分析（而不用在三分之一倍頻程圖上進行A計權修正），因此該問題關鍵就轉化為如何設計各個計權下的濾波器。

圖1 不同計權下的頻率圖（摘自模態空間）

通過查閱標準標準 JJG 188-2002聲級計檢定規程，我們可以得到A計權幅頻響應曲線（公式2）。

公式2 A計權幅頻響應

其中，f1 =20.6Hz , f2 =107.7Hz , f3 = 737.9Hz , f4 =12194Hz , A1000 =-2.000dB.我們先轉化為未取對數的公式，並且將自然頻率轉化為模擬角頻率：

公式3 轉化後的A計權幅頻響應曲線

為了設計數字濾波器，我們可以先通過上式推出S複數域的傳遞函數，然後利用雙線性變換法設計對應的IIR數字濾波器。我們首先來觀察一個一階系統：

公式4 一階系統傳遞函數

由信號與系統的知識，令可以得到系統的頻域表示，化簡併進一步表示系統的幅頻響應可以得到：

公式5 一階系統幅頻響應

這個結果是不是很像公式3中的分母一部分，而一個系統可以看成多個子系統級聯起來的，級聯的幅頻響應是乘積關係，因此我們很容易聯想到公式3是由6個一階系統、1個四階微分系統和常數增益環節構成，因此級聯後的A系統的傳遞函數為：

公式6 A計權的傳遞函數

接下來就很容易得到該模擬濾波器的分子分母係數，然後只需要利用matlab的雙線性變換函數（bilinear）就得到A計權的數字濾波器。

圖2 matlab設計後的結果

當然C計權的濾波器也可以仿照上述的方法設計，C計權的公式如下：

公式7 C計權幅頻響應

其中f1與f4和上邊h取值相同， $C_{1000}=-0.062$ .

下面將設計兩種計權的matlab代碼貼出來：

%A計權function [B,A] = adsgn(Fs); f1 = 20.598997; f2 = 107.65265;f3 = 737.86223;f4 = 12194.217;A1000 = 1.9997;pi = 3.14159265358979;NUMs = [ (2*pi*f4)^2*(10^(A1000/20)) 0 0 0 0 ];DENs = conv([1 +4*pi*f4 (2*pi*f4)^2],[1 +4*pi*f1 (2*pi*f1)^2]); DENs = conv(conv(DENs,[1 2*pi*f3]),[1 2*pi*f2]);[B,A] = bilinear(NUMs,DENs,Fs);

%C計權function [B,A] = adsgn2(Fs); f1 = 20.598997; f4 = 12194.217;C1000 = 0.062;pi = 3.14159265358979;NUMs = [ (2*pi*f4)^2*(10^(C1000/20)) 0 0 ];DENs = conv([1 +4*pi*f4 (2*pi*f4)^2],[1 +4*pi*f1 (2*pi*f1)^2]); [B,A] = bilinear(NUMs,DENs,Fs);