1. 引言

請(qǐng)參考我們前一篇文章《數(shù)字濾波器（1）——陷波濾波器》的引言段落。我們還是圍繞著下面的兩個(gè)公式展開。公式中，（式-1）時(shí)時(shí)域的差分方程；（式-2）時(shí)z復(fù)頻域的系統(tǒng)函數(shù)形式。

前面提到的陷波濾波器，其功能主要是在相關(guān)信號(hào)中削弱或者消除特定頻率的信號(hào)。本文將圍繞但不限于梳狀濾波器進(jìn)行展開。其中，梳狀濾波器，一方面可以濾除特定的頻率（尤其是特定頻率的諧波），另一方面也可以在信號(hào)中對(duì)指定頻率及其倍頻的信號(hào)進(jìn)行揀選。例如圖-1所示的濾除指定諧波成分的梳狀濾波器的幅頻圖。

圖-1 梳狀濾波器幅頻響應(yīng)示意圖

2.關(guān)于線性移不變系統(tǒng)和線性相位系統(tǒng)

同樣的，關(guān)于線性移不變系統(tǒng)的特點(diǎn)，可以查看我們之前的文章《數(shù)字濾波器（1）——陷波濾波器》。這里為什么要提出另外一個(gè)術(shù)語(yǔ)：線性相位，這事關(guān)信號(hào)傳遞的保真。

在設(shè)計(jì)濾波器時(shí)，選擇具有線性相位(或近似線性相位)的系統(tǒng)非常重要。在許多應(yīng)用中，如音頻處理或圖像處理，不僅需要關(guān)注濾波器的幅度響應(yīng)，也需要關(guān)注濾波器的相位響應(yīng)。在通過濾波器的信號(hào)的所有頻率得到相同的相對(duì)延遲（也就是線性相位）可以避免產(chǎn)生相位畸變。相位畸變會(huì)在信號(hào)中引入不想要的時(shí)間延遲，這對(duì)于音頻信號(hào)處理來說尤其重要，因?yàn)樗赡軙?huì)改變聲音的定位感知。

例如，對(duì)于一個(gè)音頻信號(hào)處理應(yīng)用，如果不同頻率的聲音延遲不一致，那么聲音的定位就會(huì)受影響。同樣，對(duì)于圖像處理，如果圖像的各個(gè)頻率部分經(jīng)過的延遲不同，可能會(huì)影響圖像的質(zhì)量，造成圖像模糊等。

所以在設(shè)計(jì)濾波器時(shí)，控制相位響應(yīng)以達(dá)到線性或近似線性相位是非常關(guān)鍵的一步。線性相位的濾波器可以確保過濾之后的信號(hào)不會(huì)出現(xiàn)相位畸變，從而確保原始信號(hào)的完整性。

在信號(hào)處理中，"線性時(shí)不變系統(tǒng)" 和 "線性相位系統(tǒng)" 是兩個(gè)具有不同含義的概念。他們的定義和特性如下：

線性時(shí)不變系統(tǒng)（Linear Time-Invariant System，或稱LTI）：具有線性和時(shí)間不變性的系統(tǒng)。線性意味著系統(tǒng)對(duì)任何輸入信號(hào)的響應(yīng)，等于這個(gè)輸入信號(hào)的線性組合的相應(yīng)反應(yīng)的線性組合。時(shí)間不變性意味著如果輸入信號(hào)延后了一段時(shí)間，那么輸出信號(hào)也延后相同的時(shí)間。LTI系統(tǒng)的一個(gè)重要特性是它們的行為可以通過系統(tǒng)的脈沖響應(yīng)或者頻率響應(yīng)（即傳遞函數(shù)）完全描述。

線性相位系統(tǒng)（Linear Phase System）：系統(tǒng)的頻率響應(yīng)或傳遞函數(shù)在相位上呈線性變化的系統(tǒng)。線性相位系統(tǒng)主要特點(diǎn)是各頻率成分通過該系統(tǒng)時(shí)具有相同的群延遲，也就是說，所有頻率的信號(hào)元素都同時(shí)到達(dá)。這點(diǎn)在很多應(yīng)用中是很重要的，比如音頻信號(hào)處理和圖像處理，因?yàn)樗梢苑乐瓜辔换儭?/p>

一般的，非線性系統(tǒng)是很難實(shí)現(xiàn)線性相位的，要說線性相位系統(tǒng)需要基于FIR也不為過。由于信號(hào)可以分解成正余弦分量的方式，窺一斑可知全豹，這里不失一般性，我們看一個(gè)復(fù)指數(shù)信號(hào)通過一個(gè)單位脈沖響應(yīng)函數(shù)h(n)的線性系統(tǒng)會(huì)得到什么響應(yīng)輸出。

通過卷積我們可以得到響應(yīng)y(n):

響應(yīng)輸出是在原信號(hào)的基礎(chǔ)上有幅值|H(ω)|倍數(shù)和附加相位arg[H(ω)]=?(ω)的影響。在很多情況下，同一個(gè)系統(tǒng)函數(shù)，不同頻率的信號(hào)輸入后，輸出時(shí)得到的延遲是不相等的。通過設(shè)計(jì)，我們完全可以控制對(duì)輸入信號(hào)的響應(yīng)在不同頻率時(shí)對(duì)應(yīng)的幅值放大倍數(shù)。比如下圖的一個(gè)示例的低通濾波器的幅頻響應(yīng)和相頻響應(yīng)。

圖-2 一般濾波器的幅頻響應(yīng)相頻響應(yīng)示意圖

上圖所示的系統(tǒng)可以讓頻率低于1MHz的信號(hào)基本保持輸入時(shí)的響應(yīng)輸出。然而，關(guān)于每個(gè)頻率對(duì)應(yīng)的相位響應(yīng)就不太確定了。如圖-2中在截止頻率范圍內(nèi)的相頻響應(yīng)部分。

基本的線性相位要求是群延遲τ為常數(shù)：

所以這里對(duì)應(yīng)的相位響應(yīng)：

τ和β都為常數(shù)。要實(shí)現(xiàn)線性相位的特性，IIR是指望不上的，多為FIR系統(tǒng)。通過特別的單位脈沖響應(yīng)序列的偶對(duì)稱或者奇對(duì)稱設(shè)計(jì)，讓頻率響應(yīng)中的相頻響應(yīng)和幅頻響應(yīng)部分可以分離，從而我們?cè)陬l率響應(yīng)的結(jié)果中看到了一個(gè)群延遲τ(e^jω)為所期待的定值結(jié)果。

回到時(shí)域中后，令:

則：

上式中，n和τ都是離散序列中表示時(shí)間的。這就意味著輸入的所有的有效信號(hào)在經(jīng)過系統(tǒng)之后，都在原信號(hào)上增加了一個(gè)相同的延遲τ——這就是我們所說的群延遲的來歷。

加這一段關(guān)于線性相位的補(bǔ)充說明，是要各位留意在實(shí)際的濾波器設(shè)計(jì)過程中，需要根據(jù)應(yīng)用來選擇適用的設(shè)計(jì)，否則反而會(huì)引起信號(hào)的失真。如下圖所示的兩種濾波器的群延遲比較，線性相位的基于FIR，非線性相位的基于IIR。(跑題沒有？）

圖-3 線性相位和非線性相位的群延遲比較

2.窄帶選頻濾波器(Narrowband Selective Filter)

（還在標(biāo)題范圍內(nèi)不？）在開始梳狀濾波器之前，我們?cè)夙槺闾嵋幌抡瓗нx頻濾波器。和陷波濾波器相反，選頻濾波是要保留在所選頻帶內(nèi)的信號(hào)。

窄帶選頻濾波器屬于帶通濾波器的一個(gè)特例。帶通濾波器允許一定范圍內(nèi)的頻率通過，而阻止其它頻率的信號(hào)。如果這個(gè)范圍非常窄，這樣的帶通濾波器就可以被稱為窄帶選頻濾波器。它的主要作用是只允許非常特定的頻率的信號(hào)通過，而其它頻率的信號(hào)都被阻止。

比如我們打算通過過濾周期方波的形式，用軟件的方式過濾諧波，實(shí)現(xiàn)基波頻率正弦波的輸出。

圖-4 方波信號(hào)濾波前后的頻譜圖

圖-5 方波信號(hào)濾波前后的波形圖和濾波器的幅頻響應(yīng)圖

上面的示例中，濾波器的幅頻響應(yīng)并不是設(shè)想中的那樣在指定頻率處加一個(gè)非常窄的窗。我們比較一下理想窄帶選頻濾波器上加不同的窗口函數(shù)之后的濾波器幅頻響應(yīng)圖。

圖-6 濾波器加矩形窗和加漢明窗的濾波器幅頻響應(yīng)比較

理想在這里反而是骨感的（那些支棱起的旁瓣），而現(xiàn)實(shí)所使用的卻略顯豐滿，不過仍然是折衷的選擇。我們必須要在旁瓣大小和過渡區(qū)寬度之間做一個(gè)權(quán)衡選擇。

濾波器加漢明窗測(cè)試代碼如下：

from scipy.signal import firwin, freqz
import matplotlib.pyplot as plt
import numpy as np
# 設(shè)定采樣率、帶通頻率范圍
fs = 2000.0 # 采樣頻率(Hz)
f1, f2 = 590, 610 # 通帶頻率范圍(Hz)


# 設(shè)定濾波器長(zhǎng)度
numtaps = 61 # (濾波器系數(shù)數(shù)量，也就是“抽頭”)


# 用firwin設(shè)計(jì)帶通濾波器，先選用矩形窗
taps_rectangular = firwin(numtaps, [f1, f2], pass_zero=False, fs=fs, window='boxcar')


# 再設(shè)計(jì)使用漢明窗的濾波器
taps_hamming = firwin(numtaps, [f1, f2], pass_zero=False, fs=fs, window='hamming')


# 計(jì)算兩個(gè)濾波器的頻率響應(yīng)
w1, h1 = freqz(taps_rectangular, worN=8000, fs=fs)
w2, h2 = freqz(taps_hamming, worN=8000, fs=fs)


# 創(chuàng)建兩個(gè)濾波器幅度響應(yīng)圖
plt.figure(figsize=(10,6))
plt.plot(w1, abs(h1), label='Rectangular Window')
plt.plot(w2, abs(h2), label='Hamming Window')
plt.title('Frequency Response')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Gain')
plt.legend(loc='best')
plt.grid(True)


# 創(chuàng)建兩個(gè)濾波器相位響應(yīng)圖
plt.figure(figsize=(10,6))
plt.plot(w1, np.angle(h1), label='Rectangular Window')
plt.plot(w2, np.angle(h2), label='Hamming Window')
plt.title('Phase Response')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Phase (Radians)')
plt.legend(loc='best')
plt.grid(True)


plt.show()

3.梳狀濾波器(Comb Filter)

再回到梳狀濾波器這個(gè)話題。這要是寫作文，很擔(dān)心會(huì)走題啊。但是寫關(guān)于應(yīng)用方面的東東，其實(shí)更多考慮的是面面俱到，因?yàn)閱为?dú)的話題往往會(huì)片面。

梳狀濾波器，顧名思義，其頻率響應(yīng)曲線像梳子一樣，有著很多個(gè)頻率通道可以穿透或者過濾，中間隔離著獨(dú)立的頻帶。梳狀濾波器的主要特點(diǎn)是在一定頻率范圍內(nèi)有多個(gè)通帶或者阻帶，每個(gè)通帶/阻帶均勻分布。

特性：

梳狀濾波器具有良好的二次諧波抑制效果。在應(yīng)用了梳狀濾波器的通信系統(tǒng)中，可以有效避免二次諧波的產(chǎn)生，提高整個(gè)系統(tǒng)的信號(hào)質(zhì)量。

梳狀濾波器可以避免頻率疊加現(xiàn)象的出現(xiàn)，保證了信號(hào)在傳輸過程中的可靠性。

在相同的工作頻率和帶寬下，梳狀濾波器相比于傳統(tǒng)濾波器會(huì)有更高的功率轉(zhuǎn)換效率。

應(yīng)用場(chǎng)合：

無線通信：梳狀濾波器是無線通信系統(tǒng)中的關(guān)鍵部件。在射頻通信鏈路中，經(jīng)常需要對(duì)發(fā)射機(jī)和接收機(jī)的信號(hào)進(jìn)行濾波處理，以消除頻率上的雜波，提高通信的清晰度和穩(wěn)定性。

微波設(shè)備：梳狀濾波器在微波信號(hào)的處理和傳輸中也起著重要的作用。由于微波頻段的信號(hào)具有很大的帶寬，因此，梳狀濾波器可以針對(duì)特定的頻率進(jìn)行過濾，有效避免信號(hào)相互干擾。

射頻系統(tǒng)：在射頻系統(tǒng)中，梳狀濾波器主要用于選擇性地阻止某些頻率的信號(hào)，從而改善射頻設(shè)備的性能。

數(shù)字電視和無線電視接收設(shè)備中也經(jīng)常使用到梳狀濾波器。

我們通過模擬信號(hào)以及濾波處理后的信號(hào)之間進(jìn)行對(duì)比，來了解這種類型濾波器的特點(diǎn)。

如果需要保留基頻及其諧波對(duì)應(yīng)的信號(hào)，那么就可以用到梳狀濾波器了。

模擬信號(hào)中有50Hz及其諧波的部分需要保留，保留50Hz和200Hz的有用信號(hào)。事先聲明，這個(gè)濾波器是基于IIR生成的。先看結(jié)果，后查代碼[2]。

圖-7-1 信號(hào)濾波前后對(duì)比（過濾后仍有50,200Hz頻率的信號(hào)疊加）

為了更好地感覺梳狀選頻濾波的效果，我們留單個(gè)頻率的效果看看。

圖-7-2 信號(hào)濾波前后對(duì)比（去除另外一路信號(hào)，只有50Hz的信號(hào)過濾后）

圖-8 信號(hào)濾波前后的頻譜圖（留意50Hz和200Hz對(duì)應(yīng)被保留的信號(hào)）

圖-8中，50Hz及其諧波被保留，而其他頻率的信號(hào)則被削弱。

圖-9 該梳狀濾波器幅、相頻率響應(yīng)圖

這是一個(gè)基于IIR系統(tǒng)設(shè)計(jì)的濾波器。所以在相頻響應(yīng)部分可以看到并不是線性相位的。演示代碼中，我們可以輸出該濾波器的相應(yīng)系數(shù)。對(duì)應(yīng)的系統(tǒng)函數(shù)仍然可以簡(jiǎn)化為（式-2）所示。

圖-1所示的陷波梳狀濾波器和圖-9所示的峰值梳狀濾波器的幅頻響應(yīng)來看，兩種濾波器的功能剛好是相反的。所有的極點(diǎn)和零點(diǎn)都在演示代碼中可以輸出。我們可以看一下各自的零點(diǎn)、極點(diǎn)在z平面的分布情況。

圖-10 陷波（左）和峰值（右）梳狀濾波器的零極點(diǎn)分布圖

Z平面零極點(diǎn)分布特點(diǎn)：零點(diǎn)在單位圓上，極點(diǎn)離單位圓的距離小于1。其中：

陷波梳狀濾波器的目標(biāo)是對(duì)一種或幾種特定頻率的信號(hào)進(jìn)行衰減或者濾除。所以，在對(duì)應(yīng)于這些頻率的單位圓上，放置了零點(diǎn)來最大限度地衰減掉這些頻率的信號(hào)。同時(shí)，為了使得在非濾波頻率處信號(hào)可以通過，所以在零點(diǎn)附近的單位圓內(nèi)部，放置了極點(diǎn)來抵消在其他頻率處零點(diǎn)引起的衰減。至于極點(diǎn)到單位圓的距離，即半徑ρ<1，決定了濾波器在濾波頻率處的衰減深度以及過渡帶的寬度。 ? ?

峰值梳狀濾波器剛好相反，它的目標(biāo)是放大一些特定頻率的信號(hào)，或者在某些頻率處衰減。所以，它在濾波頻率對(duì)應(yīng)的位置放置了極點(diǎn)來放大這些頻率的信號(hào)，同時(shí)在極點(diǎn)附近的單位圓上放置零點(diǎn)來抵消在其他頻率處極點(diǎn)引起的放大。極點(diǎn)離單位圓的距離，即半徑ρ<1，決定了濾波器在濾波頻率處的放大程度以及過渡帶的寬度。

在后面的演示代碼中，我們可以看到該處理函數(shù)的采樣頻率是1000Hz。所以這個(gè)頻譜圖也只到500Hz。那為什么不用更高的采樣頻率？有以下喜憂參半的幾點(diǎn)要考慮：

采樣頻率越高，就意味著同樣的時(shí)間內(nèi)應(yīng)用系統(tǒng)需要更多的緩存資源來保存和處理更多的數(shù)據(jù)。

其次，按照奈奎斯特采樣定理，你的采樣頻率如果為Fs，則能夠有效采樣和還原的信號(hào)頻率上限是Fs/2。所以，當(dāng)你提高采樣頻率時(shí)，能夠處理的信號(hào)頻率范圍也變得更廣，可以過濾的基頻倍數(shù)也越多和越高。

再有需要考慮，是否真的需要更高的采樣頻率？這完全取決于應(yīng)用需要。很多情況下，諧波頻率越高，對(duì)應(yīng)的能量也越小，在多次諧波之后，對(duì)應(yīng)諧波所附能量幾乎可以忽略的情況下，就沒有必要再額外增加這些資源。

雖然梳狀濾波器的信號(hào)傳遞架構(gòu)仍然可以參考下圖，不過因?yàn)槠湟?guī)律性的梳齒特征的幅頻響應(yīng)可以由（式-2）簡(jiǎn)化到（式-3）所示，仍然是IIR系統(tǒng)，信號(hào)傳遞過程也可以相應(yīng)簡(jiǎn)化。

圖-11 濾波器的信號(hào)傳遞

（式-3）中分子對(duì)于零點(diǎn)的正負(fù)的選擇，將取決于我們對(duì)所選頻率的所需的位置。

如果是“-”號(hào)，那么梳狀濾波器必有零點(diǎn)對(duì)應(yīng)于直流DC信號(hào)和高頻段；

如果是“+”號(hào)，那梳狀濾波器在DC和高頻段將沒有零點(diǎn)。

在以上兩種情況下，對(duì)于梳狀濾波器的零極點(diǎn)的分布將會(huì)有直接影響。實(shí)際應(yīng)用需要根據(jù)需要進(jìn)行調(diào)整。了解這些，對(duì)于調(diào)用現(xiàn)成的函數(shù)，也可以提供了解和參考。

理論上，我們可以用陷波濾波器或窄帶選頻濾波器串并聯(lián)的方式來搭建梳狀濾波器，但是一定要了解這種實(shí)現(xiàn)方式對(duì)于資源和時(shí)延，以及是否線性相位等方面對(duì)應(yīng)用的影響。

4. 小結(jié)

在本文中，我們探討了包括梳狀濾波器在內(nèi)的幾種濾波器的設(shè)計(jì)和應(yīng)用，以及仿真，同時(shí)簡(jiǎn)單探討了以下線性相位系統(tǒng)的重要性。在各種傳感器的應(yīng)用過程中，總是會(huì)碰到各種狀況，我們需要使用者根據(jù)實(shí)際情況進(jìn)行應(yīng)對(duì)。

另外，實(shí)際應(yīng)用過程中，還有哪些條件需要考慮的？比如，如果數(shù)據(jù)采樣的時(shí)鐘因?yàn)槟撤N原因并不準(zhǔn)有影響嗎？有無補(bǔ)償?shù)姆绞侥兀?/p>

安費(fèi)諾傳感器（Amphenol Sensors）提供諸多滿足您應(yīng)用的可靠、準(zhǔn)確高精度的傳感器，但如果實(shí)際應(yīng)用過程中碰到一些異常的信號(hào)，可以和我們一起來分析原因，解決問題，根據(jù)實(shí)際工況來改善應(yīng)用。

峰值梳狀濾波器模擬測(cè)試代碼

import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as signal
from scipy.fft import fft


# Sampling rate and time vector
fs = 1000.0
t = np.arange(0, 0.5, 1/fs)


# Input signal
freqs = np.array([50.0, 120.0, 270.0, 322.0])
x = np.sin(2.0*np.pi*freqs[:,None]*t).sum(axis=0)


# Design filter
w0 = 50.0 # Frequency to be removed from signal (Hz)
Q = 30.0 # Quality factor
b, a = signal.iircomb(w0, Q, ftype='peak', fs=fs, pass_zero=False)


# Calculate Zeros, Poles & gain of the filter
zeros, poles, gain = signal.tf2zpk(b, a)
print('Zeros: ', zeros)
print('Poles: ', poles)
"""
#輸出零極點(diǎn)
print('Zeros: ', zeros)
for zero in zeros:
    print('Abs of Zero:',abs(zero))
    
print('Poles: ', poles)
for pole in poles:
    print('Abs of Pole:',abs(pole))
"""
print('Gain: ', gain)


# Frequency response
freq, h = signal.freqz(b, a, fs=fs)
amplitude = np.abs(h)
phase = np.angle(h)


# Apply filter
y = signal.lfilter(b, a, x)


# Normalized FFT
N = len(t)
fftx = (fft(x)/ len(t))[:N//2]
ffty = (fft(y)/ len(t))[:N//2]


frequencies = np.linspace(0, fs/2, len(fftx))


# Output to figures
plt.figure(figsize=(8, 6))


# Input signal
plt.subplot(2, 1, 1)
plt.plot(t, x, color='blue')
plt.grid()
plt.title('Input signal')


# Filtered signal
plt.subplot(2, 1, 2)
plt.plot(t, y, color='red')
plt.title('Filtered signal')
plt.grid()
plt.tight_layout()
plt.show()


# FFT of input signal
plt.figure(figsize=(8, 6))


plt.subplot(2, 1, 1)
#plt.semilogx(t*fs, np.abs(fftx))
plt.semilogx(frequencies, np.abs(fftx),color='blue')
plt.grid()
plt.title('FFT of input signal')


# FFT of filtered signal
plt.subplot(2, 1, 2)
#plt.semilogx(t*fs, np.abs(ffty))
plt.semilogx(frequencies, np.abs(ffty),color='red')
plt.title('FFT of filtered signal')
plt.grid()
plt.tight_layout()
plt.show()


# Plot frequency and phase response
plt.figure(figsize=(8, 4))


plt.subplot(2, 1, 1)
plt.plot(freq, amplitude, color='blue')
plt.grid()
plt.title('Frequency response')


plt.subplot(2, 1, 2)
plt.plot(freq, phase, color='red')
plt.grid()
plt.title('Phase response')


plt.tight_layout()
plt.show()


# 極坐標(biāo)圖輸出
plt.figure(figsize=(8,8))
plt.scatter(0, 0, marker='+', color='green')
plt.scatter(np.real(poles), np.imag(poles), marker='x', color='blue', label='Poles')
plt.scatter(np.real(zeros), np.imag(zeros), marker='o', color='red', label='Zeros')
# 參數(shù)(0,0)代表圓心坐標(biāo)，1代表半徑
circle = plt.Circle((0, 0), 1, fill = False, linestyle='dashed') 
plt.gca().add_patch(circle)


plt.title('Pole / Zero Plot')
plt.xlabel('Real')
plt.ylabel('Imaginary')
plt.legend()
plt.grid()
plt.axis('equal') # equal scaling in both x and y
plt.show()