數字信号處理5：FIR濾波器設計

文章目錄

• 1. 濾波器初識
• 2. 最直覺的濾波方式：頻域濾波
• 3. 傅裡葉變換中的加窗
• 4. FIR濾波器設計
• 5. 總結

之前的一系列部落格中，詳細分解了從卷積到FFT的相關知識，不過那些屬于理論，是為了讓我們清楚認識到信号處理的本質。本篇部落格将會詳細講解數字信号處理最廣泛的應用——濾波器。

注意，本章所采用的dB标準為

d B = 20 l o g V V b a s e dB = 20log\frac{V}{V_{base}} dB=20logVbase​V​

其中， V b a s e V_{base} Vbase​表示某種作為标準的幅度值，它或者是電平，或者是信号振幅。在這裡，我們通常将它認為是原信号的“大小”，這個大小不光可以是信号在時域的幅度，也可以是在頻域的頻率分量的大小。

1. 濾波器初識

濾波器，如同字面意思是過濾信号的某些頻率分量。入門時書上最常說的幾種濾波器，比如低通、帶通、高通、帶阻等濾波器。如下是一個低通濾波器。

它将信号中的中高頻分量降低了60dB。按照dB的定義，也就是中高頻信号的幅度降低為原來的0.001倍。

實際上，更為廣泛的濾波器，并不隻有“濾”的作用，它同樣還可以對某些頻段進行增強，如下濾波器：

它将信号中的中高頻部分提高了60dB。按照dB的定義，也就是中高頻信号的幅度提高到原來的1000倍。

2. 最直覺的濾波方式：頻域濾波

回想一下《數字信号處理2：傅裡葉變換》部落格中講的傅裡葉變換的直覺了解那部分。傅裡葉變換之後的結果是一組虛數，這組虛數的相位反應的就是原信号中各頻率正弦波的相位，而這些虛數的絕對值某種程度上反應了對應頻率的信号在原始信号中的大小。是以，我們隻要對頻域上的這些虛數的絕對值進行調節，再作傅裡葉逆變換，即可以調整原信号中的某些頻率分量的大小，也就成為了濾波的一種手段。步驟通常如下：

1. 對信号加窗
2. 對加窗之後的結果作傅裡葉變換。
3. 根據你的目的對傅裡葉變換結果的某些頻率上的分量調節絕對值。
4. 将調節完畢後的結果作傅裡葉逆變換。

這便是頻域濾波。

3. 傅裡葉變換中的加窗

在上一部分中提到了對信号進行加窗，為什麼要進行加窗操作？加的又是什麼窗？

要了解加窗的意義，有兩種方式。

首先，回想一下部落格《數字信号處理2：傅裡葉變換》中對傅裡葉變換的推導過程，對于傅裡葉變換所做的所有推導，都是基于無限長周期信号進行的。我們假設用來做傅裡葉變換的信号是從周期信号中截取的一整個周期。是以，如果這段信号前後進行拼接成我們假想的周期信号後，中間沒有突變點（包括斜率突變、斷點等情況），那麼直接進行傅裡葉變換是沒有問題的。而如果前後進行拼接後，中間出現了突變點，那麼從直覺就可以知道，這個突變點中包含了相當多的頻率分量。但單從我們截取的這段信号來說，是不應該有這些多餘的頻率分量的。

下面舉個例子：

構造一個頻率為20的周期信号。首先，我們讓這段信号有整數個周期。使用matlab。

T = 200;

x = [0:1:T-1];
sinWave = sin(2 * pi / T * 20 * x);
plot(x, sinWave);
title("頻率為20的正弦波");
b = fft(sinWave);
plot(x, abs(b));
title("整數個周期正弦波的頻譜")
           

代碼中生成的sin波是完整的20個周期，采樣點200個，周期是10。繪制出來的圖像如下：

可以看到，頻譜上僅在周期為20的位置有值。

如果不給整數個周期。比如隻給19.5個周期。顯然，這樣的信号前後拼接成周期信号後，必然是會有突變點的。

x = [0:1:195];
sinWave = sin(2 * pi / T * 20 * x);
plot(x, sinWave);
title("19.5個周期的正弦波");
b = fft(sinWave);
%b = circshift(b', fix(length(sinWave) / 2))';
plot(x, abs(b));
title("非整數個周期正弦波的頻譜")
           

可以看到，這時在頻率20周圍出現了其他的頻率分量，這就是頻率洩漏。

實際運用中，一段信号通常有很多頻率的正弦波組成，我們不能保證總是能截取到整數個周期，是以，通過加窗的方式，我們使信号兩端加一個漸緩的過程，使得拼接成周期信号後，沒有突變點。

舉一個例子：漢甯窗。它的形狀如下：

對于上面提到的非整數個周期的信号，我們對其進行加窗操作：

x = [0:1:195];
window = hanning(length(x))';
plot(window);
title("漢甯窗");
sinWave = sin(2 * pi / T * 20 * x) .* window;
plot(x, sinWave);
title("加窗後的19.5個周期的正弦波");
b = fft(sinWave);
%b = circshift(b', fix(length(sinWave) / 2))';
plot(x, abs(b));
title("加窗後的非整數個周期正弦波的頻譜")
           

之後信号就變成這樣

其頻譜

可以看到，加窗之後，頻率洩漏的情況改善了很多。

從另一方面了解。之前的部落格《數字信号處理4：采樣定理》中證明了一個重要的結論：時域卷積等于頻域相乘，時域相乘等于頻域卷積。

是以，将信号做傅裡葉變換，其實隐含了一個過程，就是原信号的頻譜與窗函數的頻譜進行卷積。而如果不想改變原信号頻譜，最理想的情況就是讓它與一個沖激信号進行卷積，這樣原信号的頻譜其實隻是發生了位移。這樣一來，我們就希望能有一個窗函數，它的頻譜是一個沖激函數，或者盡量接近沖激函數。

對于一段信号，不加窗等價于加矩形窗。我們不妨看一下矩形窗的頻譜：

a=zeros(1, 400);
a(100:299) = ones(1, 200);
b = fft(a);
b = circshift(b', 200)';
plot(a);
title("矩形窗")

plot([-199:1:200], abs(b));
title("矩形窗頻譜")
           

顯然，矩形窗有很多旁瓣，與這樣的信号在頻域上做卷積，會改變原信号的頻譜。

而看一下漢甯窗：

a=zeros(1, 400);
a(100:299) = hanning(200)';

plot(a);
title("漢甯窗")
b = fft(a);
b = circshift(b', 200)';

plot([-199:1:200],abs(b));
title("漢甯窗頻譜")
           

可見漢甯窗的頻譜相當接近一個沖激函數。

至此，窗函數的意義已經明了：消除傅裡葉變換過程中的頻率洩漏。

4. FIR濾波器設計

濾波器是一個線性時不變系統，我們使用一個差分方程來表示該系統的信号傳輸特性：

y [ n ] = − ∑ k = 1 N a k y [ n − k ] + ∑ k = 0 M − 1 b k x [ n − k ] y[n] = -\sum_{k=1}^{N}a_ky[n-k] + \sum_{k=0}^{M-1}b_kx[n-k] y[n]=−k=1∑N​ak​y[n−k]+k=0∑M−1​bk​x[n−k]

如果 a k ≠ 0 a_k \neq 0 ak​​=0，那麼目前輸出不僅與目前輸入有關，還有過去的 M − 1 M-1 M−1個輸入，以及過去的 N N N個輸出有關。這樣的濾波器成為無限沖激響應濾波器（IIR Filter）。如何了解無限這個詞，目前輸出與過去的輸出有關，那麼可以預見，第一個輸入對後面所有的輸出都會有影響。

如果 a k = 0 a_k = 0 ak​=0，那麼目前輸出僅與目前輸入和過去的 M − 1 M-1 M−1個輸入有關。這樣的濾波器被成為有限沖激響應濾波器（FIR Filter）。

這兩種濾波器各有優缺點，IIR的優點在于計算快，能以較少的階數達到性能要求。但設計比較複雜，而且難以設計出具有準确頻率響應的濾波器，另外，IIR濾波器的相位不可能是線性的。

而FIR濾波器則反了過來，設計簡單，能快速設計出具有精确線性相位以及需要的頻率響應的濾波器。但是需要較高的階數才能達到濾波要求。

設 h [ n ] h[n] h[n]是濾波器， H [ k ] H[k] H[k]是 h [ n ] h[n] h[n]的傅裡葉變換。

我們從頻域濾波開始，既然從頻域調節就可以達到濾波器的目的，那麼幹脆将頻域上調節的系數看做 H [ k ] H[k] H[k]，然後用傅裡葉逆變換求出 h [ n ] h[n] h[n]，頻域相乘等于時域卷積，使用 h [ n ] h[n] h[n]與信号 x [ n ] x[n] x[n]進行卷積即可。

首先要明确的一個點，為什麼濾波器要是線性相位？按照傅裡葉變換公式

X [ k ] = ∣ X [ k ] ∣ e − j a k w 0 X[k] = |X[k]|e^{-ja_kw_0} X[k]=∣X[k]∣e−jak​w0​

然後有一濾波器 h [ n ] h[n] h[n]，其頻譜是

H [ k ] = ∣ H [ k ] ∣ e − j b k w 0 H[k] = |H[k]|e^{-jb_kw_0} H[k]=∣H[k]∣e−jbk​w0​

當 X [ k ] X[k] X[k]與 H [ k ] H[k] H[k]的每一項相乘。

Y [ k ] = H [ k ] X [ k ] = ∣ H [ k ] ∣ ∣ X [ k ] ∣ e − j ( a k + b k ) w 0 Y[k] = H[k]X[k] = |H[k]||X[k]|e^{-j(a_k + b_k)w_0} Y[k]=H[k]X[k]=∣H[k]∣∣X[k]∣e−j(ak​+bk​)w0​

我們已經知道傅裡葉變換是可逆的。是以對于 X [ k ] X[k] X[k]來說

x [ n ] = 1 N ∑ k = < N > ∣ X [ k ] ∣ e − j a k w 0 e j k w 0 n = 1 N ∑ k = < N > ∣ X [ k ] ∣ e j ( k n − a k ) w 0 x[n] = \frac{1}{N}\sum_{k=<N>}|X[k]|e^{-ja_kw_0}e^{jkw_0n} = \frac{1}{N}\sum_{k=<N>}|X[k]|e^{j(kn - a_k)w_0} x[n]=N1​k=<N>∑​∣X[k]∣e−jak​w0​ejkw0​n=N1​k=<N>∑​∣X[k]∣ej(kn−ak​)w0​

得出的結果必然是實數。

那麼對于 Y [ k ] Y[k] Y[k]：

y [ n ] = 1 N ∑ k = < N > ∣ X [ k ] H [ k ] ∣ e j ( k n − a k − b k ) w 0 y[n] = \frac{1}{N}\sum_{k=<N>}|X[k]H[k]|e^{j(kn - a_k - b_k)w_0} y[n]=N1​k=<N>∑​∣X[k]H[k]∣ej(kn−ak​−bk​)w0​

假設 ∣ H [ k ] ∣ = 1 |H[k]|=1 ∣H[k]∣=1，相當于沒有對 X [ k ] X[k] X[k]進行調節，此時隻看 b k b_k bk​。

當 b k = 0 b_k = 0 bk​=0時，就退化成 x [ n ] x[n] x[n]的傅裡葉逆變換，這個時候 x [ n ] x[n] x[n]自然是沒有問題，但是我們知道對于非全0的 h [ n ] h[n] h[n]來說， H [ k ] H[k] H[k]不可能是零相位的。

如果 b k b_k bk​是個常數，也就是 H [ k ] H[k] H[k]的相位是線性的，很顯然， y [ n ] y[n] y[n]隻是相當于将 x [ n ] x[n] x[n]延遲了 b k b_k bk​。

而如果相位是非線性的，則無法保證 y [ n ] y[n] y[n]是實數。

接下來的任務就是确定 b k b_k bk​的值。

注意，從這裡開始，就與大多數教科書講的不一樣了

設

H [ k ] = H a [ k ] e − j b k w 0 H[k] = H_a[k]e^{-jb_kw_0} H[k]=Ha​[k]e−jbk​w0​

其中， H a [ k ] H_a[k] Ha​[k]是其 H [ k ] H[k] H[k]的幅度響應，是以 ∣ H a [ k ] ∣ = ∣ H [ k ] ∣ |H_a[k]| = |H[k]| ∣Ha​[k]∣=∣H[k]∣。

由傅裡葉變換公式可知，對于實序列的傅裡葉變換，有

H [ k ] = H ∗ [ − k ] = H ∗ [ N − 1 − k ] H[k] = H^*[-k] = H^*[N-1-k] H[k]=H∗[−k]=H∗[N−1−k]

是以， ∣ H [ k ] ∣ = ∣ H [ − k ] ∣ |H[k]| = |H[-k]| ∣H[k]∣=∣H[−k]∣，意味着 H a [ k ] H_a[k] Ha​[k]既可以是偶函數也可以是奇函數。

很多教材将 H a [ k ] H_a[k] Ha​[k]是奇函數還是偶函數區分開來，并且由于濾波器長度N可以為奇數或偶數，是以衍生了四種情況，記起來頗為麻煩。

實際上，我們隻要考慮 H a [ k ] H_a[k] Ha​[k]是偶函數的情況，因為它已經可以構造任何形式的FIR濾波器。而奇函數的形式卻無法構造調節低頻的濾波器。

第一步：構造幅度響應

這一步是我們根據濾波器要求來構造的。舉例來說，信号采樣率是16kHz，我們要将4kHz以上的信号降低60dB。

這個時候，根據濾波器長度來區分不同的情況。

如果濾波器長度為20，那麼在0點處對應的是0Hz，每兩個點頻率差為 16000 / 20 = 800 16000/20=800 16000/20=800Hz。那麼第4點對應3200Hz，第5點對應4000Hz，第10點對應的是8000Hz，剛好是采樣率為16k能表示的最大信号頻率。

要将4kHz以上的信号降低60dB，則第5點到第10點，其幅度響應都為0.001（通過dB換算）。而第0點到第4點，幅度響應都為1。

現在确認了11個點，那剩下的9個點呢？别忘了此時幅度響應是偶函數，而且傅裡葉變換後的頻譜現在是周期為20的周期函數。是以第11點與第9點相同，一直往後，第19點與第1點相同。

是以得到的幅度響應為

如果濾波器長度為21，那麼每兩個點頻率差為 16000 / 21 ≈ 762 16000/21\approx762 16000/21≈762Hz。按照同樣的方式，可以得到下表：

注意其中3810Hz處将幅度響應設為0.01，因為過于陡峭的變化會導緻吉布斯現象出現抖動。其實以上兩個響應我沒有經過檢驗， 可能也會存在吉布斯現象。

第二步：構造相位響應

之前已經證明過FIR濾波器必須是線性相位的。但是我們不知道線性系數是多少。這裡我們先假設它是1。

使用上面的偶數長度濾波器的幅度響應。不過由于濾波器長度太短，為了防止吉布斯現象，我們降低6dB即可，這樣，對應的幅度響應約為0.5。

H_amp = [1, 1, 1, 1, 1, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1, 1, 1, 1];

N = length(H_amp);

w0 = 2 * pi / length(H_amp);

theta = (w0 * 1) .* [0 : 1 : N-1];

H_phase = exp(1i .* theta);

H = H_amp .* H_phase;

h = ifft(H);

plot(real(h));
title("h的實部");

plot(imag(h));
title("h的虛部");
           

顯示如下：

虛部是一些很小的數，這可能是由于計算誤差引起的，可以忽略掉。

但是實部是不太對的，因為畢竟之後我們要加窗的，這樣一加窗，這個濾波器肯定就變了，我們至少要保證起主要作用的那些值比較大的點在加窗之後要保留下來。是以，我們不妨将其循環移位 N / 2 N/2 N/2。再來看效果。

修改代碼如下，對h進行循環移位，之後加漢甯窗，并且使用freqz函數檢視濾波器的頻域響應。

H_amp = [1, 1, 1, 1, 1, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1, 1, 1, 1];

N = length(H_amp);

w0 = 2 * pi / length(H_amp);

theta = (w0 * 1) .* [0 : 1 : N-1];

H_phase = exp(1i .* theta);

H = H_amp .* H_phase;

h = ifft(H);

h = circshift(h',fix(N/2))';

plot(real(h));
title("h的實部");

plot(imag(h));
title("h的虛部");
           

其頻域響應為：

可以看到，非常完美的線性相位，并且濾波器名額也達到我們的要求。

其實對于循環移位這個操作，相當于把h[n]延遲了N/2。由傅裡葉變換的位移定理：

∑ n = < N > x [ n − a ] e − j k w 0 n = ∑ n = < N > x [ n − a ] e − j k w 0 ( n − a ) e − j k w 0 a = X [ k ] e − j k w 0 a \sum_{n=<N>}x[n-a]e^{-jkw_0n} = \sum_{n=<N>}x[n-a]e^{-jkw_0(n-a)} e^{-jkw_0a} = X[k]e^{-jkw_0a} n=<N>∑​x[n−a]e−jkw0​n=n=<N>∑​x[n−a]e−jkw0​(n−a)e−jkw0​a=X[k]e−jkw0​a

是以，對于FIR濾波器來說，線性相位的系數是N/2。若N為奇數，則為(N-1)/2。線性相位系數必須是正數。

之前我們實際上已經給h[n]延遲了1，因為相位的系數是1。這次我們直接将1改為N/2，并且去掉後面的循環移位操作。

H_amp = [1, 1, 1, 1, 1, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1, 1, 1, 1];

N = length(H_amp);

w0 = 2 * pi / length(H_amp);

theta = (w0 * fix(N/2)) .* [0 : 1 : N-1];

H_phase = exp(1i .* theta);

H = H_amp .* H_phase;

h = ifft(H);

%h = circshift(h',fix(N/2))';


plot(real(h));
title("h的實部");

plot(imag(h));
title("h的虛部");

window = hanning(N);

h = h .* window';

freqz(real(h));
title("h的頻率響應");
           

最終的頻率響應是

5. 總結

好了，至此FIR濾波器的構造方式已經講完了，實際上關于濾波器名額還有很多細節的東西需要學習，比如濾波器長度、截止頻率等。

總結一下構造FIR濾波器的步驟：

1. 根據要求構造幅度響應。
2. 構造相位響應，必須是線性相位，線性系數為：N為偶數時是N/2；N為奇數時是(N-1)/2。
3. 幅度響應和相位響應合起來成為濾波器頻譜 H [ k ] H[k] H[k]。
4. 将 H [ k ] H[k] H[k]進行傅裡葉逆變換并取其實部得到濾波器 h [ n ] h[n] h[n]。
5. 對 h [ n ] h[n] h[n]加窗。

具體的matlab代碼可以檢視我的github：FIR-filter。部落格中使用的示例代碼都可以從中找到。包含了一個我已經寫好的fir_filter.m用來直接生成FIR濾波器。