1. 引言

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

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

圖-1 梳狀濾波器幅頻響應示意圖

2.關于線性移不變系統和線性相位系統

同樣的，關于線性移不變系統的特點，可以查看我們之前的文章《數字濾波器（1）——陷波濾波器》。這里為什么要提出另外一個術語：線性相位，這事關信號傳遞的保真。

在設計濾波器時，選擇具有線性相位(或近似線性相位)的系統非常重要。在許多應用中，如音頻處理或圖像處理，不僅需要關注濾波器的幅度響應，也需要關注濾波器的相位響應。在通過濾波器的信號的所有頻率得到相同的相對延遲（也就是線性相位）可以避免產生相位畸變。相位畸變會在信號中引入不想要的時間延遲，這對于音頻信號處理來說尤其重要，因為它可能會改變聲音的定位感知。

例如，對于一個音頻信號處理應用，如果不同頻率的聲音延遲不一致，那么聲音的定位就會受影響。同樣，對于圖像處理，如果圖像的各個頻率部分經過的延遲不同，可能會影響圖像的質量，造成圖像模糊等。

所以在設計濾波器時，控制相位響應以達到線性或近似線性相位是非常關鍵的一步。線性相位的濾波器可以確保過濾之后的信號不會出現相位畸變，從而確保原始信號的完整性。

在信號處理中，"線性時不變系統" 和 "線性相位系統" 是兩個具有不同含義的概念。他們的定義和特性如下：

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

線性相位系統（Linear Phase System）：系統的頻率響應或傳遞函數在相位上呈線性變化的系統。線性相位系統主要特點是各頻率成分通過該系統時具有相同的群延遲，也就是說，所有頻率的信號元素都同時到達。這點在很多應用中是很重要的，比如音頻信號處理和圖像處理，因為它可以防止相位畸變。

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

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

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

圖-2 一般濾波器的幅頻響應相頻響應示意圖

上圖所示的系統可以讓頻率低于1MHz的信號基本保持輸入時的響應輸出。然而，關于每個頻率對應的相位響應就不太確定了。如圖-2中在截止頻率范圍內的相頻響應部分。

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

所以這里對應的相位響應：

τ和β都為常數。要實現線性相位的特性，IIR是指望不上的，多為FIR系統。通過特別的單位脈沖響應序列的偶對稱或者奇對稱設計，讓頻率響應中的相頻響應和幅頻響應部分可以分離，從而我們在頻率響應的結果中看到了一個群延遲τ(e^jω)為所期待的定值結果。

回到時域中后，令:

則：

上式中，n和τ都是離散序列中表示時間的。這就意味著輸入的所有的有效信號在經過系統之后，都在原信號上增加了一個相同的延遲τ——這就是我們所說的群延遲的來歷。

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

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

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

（還在標題范圍內不？）在開始梳狀濾波器之前，我們再順便提一下窄帶選頻濾波器。和陷波濾波器相反，選頻濾波是要保留在所選頻帶內的信號。

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

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

圖-4 方波信號濾波前后的頻譜圖

圖-5 方波信號濾波前后的波形圖和濾波器的幅頻響應圖

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

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

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

濾波器加漢明窗測試代碼如下：

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


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


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


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


# 計算兩個濾波器的頻率響應
w1, h1 = freqz(taps_rectangular, worN=8000, fs=fs)
w2, h2 = freqz(taps_hamming, worN=8000, fs=fs)


# 創建兩個濾波器幅度響應圖
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)


# 創建兩個濾波器相位響應圖
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)

再回到梳狀濾波器這個話題。這要是寫作文，很擔心會走題啊。但是寫關于應用方面的東東，其實更多考慮的是面面俱到，因為單獨的話題往往會片面。

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

特性：

梳狀濾波器具有良好的二次諧波抑制效果。在應用了梳狀濾波器的通信系統中，可以有效避免二次諧波的產生，提高整個系統的信號質量。

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

在相同的工作頻率和帶寬下，梳狀濾波器相比于傳統濾波器會有更高的功率轉換效率。

應用場合：

無線通信：梳狀濾波器是無線通信系統中的關鍵部件。在射頻通信鏈路中，經常需要對發射機和接收機的信號進行濾波處理，以消除頻率上的雜波，提高通信的清晰度和穩定性。

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

射頻系統：在射頻系統中，梳狀濾波器主要用于選擇性地阻止某些頻率的信號，從而改善射頻設備的性能。

數字電視和無線電視接收設備中也經常使用到梳狀濾波器。

我們通過模擬信號以及濾波處理后的信號之間進行對比，來了解這種類型濾波器的特點。

如果需要保留基頻及其諧波對應的信號，那么就可以用到梳狀濾波器了。

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

圖-7-1 信號濾波前后對比（過濾后仍有50,200Hz頻率的信號疊加）

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

圖-7-2 信號濾波前后對比（去除另外一路信號，只有50Hz的信號過濾后）

圖-8 信號濾波前后的頻譜圖（留意50Hz和200Hz對應被保留的信號）

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

圖-9 該梳狀濾波器幅、相頻率響應圖

這是一個基于IIR系統設計的濾波器。所以在相頻響應部分可以看到并不是線性相位的。演示代碼中，我們可以輸出該濾波器的相應系數。對應的系統函數仍然可以簡化為（式-2）所示。

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

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

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

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

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

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

采樣頻率越高，就意味著同樣的時間內應用系統需要更多的緩存資源來保存和處理更多的數據。

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

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

雖然梳狀濾波器的信號傳遞架構仍然可以參考下圖，不過因為其規律性的梳齒特征的幅頻響應可以由（式-2）簡化到（式-3）所示，仍然是IIR系統，信號傳遞過程也可以相應簡化。

圖-11 濾波器的信號傳遞

（式-3）中分子對于零點的正負的選擇，將取決于我們對所選頻率的所需的位置。

如果是“-”號，那么梳狀濾波器必有零點對應于直流DC信號和高頻段；

如果是“+”號，那梳狀濾波器在DC和高頻段將沒有零點。

在以上兩種情況下，對于梳狀濾波器的零極點的分布將會有直接影響。實際應用需要根據需要進行調整。了解這些，對于調用現成的函數，也可以提供了解和參考。

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

4. 小結

在本文中，我們探討了包括梳狀濾波器在內的幾種濾波器的設計和應用，以及仿真，同時簡單探討了以下線性相位系統的重要性。在各種傳感器的應用過程中，總是會碰到各種狀況，我們需要使用者根據實際情況進行應對。

另外，實際應用過程中，還有哪些條件需要考慮的？比如，如果數據采樣的時鐘因為某種原因并不準有影響嗎？有無補償的方式呢？

安費諾傳感器（Amphenol Sensors）提供諸多滿足您應用的可靠、準確高精度的傳感器，但如果實際應用過程中碰到一些異常的信號，可以和我們一起來分析原因，解決問題，根據實際工況來改善應用。

峰值梳狀濾波器模擬測試代碼

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)
"""
#輸出零極點
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()


# 極坐標圖輸出
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')
# 參數(0,0)代表圓心坐標，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()