6. 頻率響應

DSP (ex: 音響、通訊設備) 系統的目的都是希望達到無失真的傳遞，在輸出端重現原始訊號。因此，理想的訊號處理系統，在不同頻率範圍，需表現均勻的強度響應。頻率響應 Frequency Response 是用來衡量系統在特定頻率範圍的操作特性。

6.1 頻率響應定義

頻率響應 (Frequency Response) ：系統對於輸入訊號在特定頻率範圍的響應情形，可定義為頻率的函式，響應則包含強度 (Magnitude) 與相位角 (Phase Angle) 等等量化數據，通常是以頻譜的方式呈現，用來表示系統的操作特性。

6.2 DSP系統的頻率響應定義

DSP系統的頻率響應 (Frequency Response) 可以定義為：
$
\begin{align}
H(\omega)=\sum_{n=-\infty}^\infty h\lbrack n\rbrack\cdot e^{-j\omega n}
\end{align}
$

參考上圖

假設有一個DSP系統 $h\lbrack n\rbrack$，當輸入訊號 $x\lbrack n\rbrack$，產生 $y\lbrack n\rbrack$ 的訊號輸出

可寫成卷積 (Convolution) 公式如下：

$
\begin{align}
y\lbrack n\rbrack &= h\lbrack n\rbrack\ast x\lbrack n\rbrack\\
&=\sum_{k=-\infty}^\infty h\lbrack k\rbrack\cdot x\lbrack n-k\rbrack
\end{align}
$

此時，若輸入的訊號為一個複指數訊號，即 $x\lbrack n\rbrack=e^{j\omega n}$

⇒ $
\begin{align}
y\lbrack n\rbrack &= h\lbrack n\rbrack\ast x\lbrack n\rbrack\\
&=\sum_{k=-\infty}^\infty h\lbrack k\rbrack\cdot x\lbrack n-k\rbrack\\
&=\sum_{k=-\infty}^\infty h\lbrack k\rbrack\cdot e^{j\omega(n-k)}\\
&=\left(\sum_{k=-\infty}^\infty h\lbrack k\rbrack\cdot e^{-j\omega k}\right)\cdot e^{j\omega(n)}
\end{align}
$

則 $H(\omega)=\sum_{k=-\infty}^\infty h\lbrack k\rbrack\cdot e^{-j\omega k}$

因此觀察上述推導，可以發現所謂的頻率響應，就只是DSP系統 $h\lbrack n\rbrack$ 求取 離散時間傅立葉轉換 (Discrete-Time Fourier Transform)的結果。

問題討論：

為什麼要使用複指數訊號，即 $x\lbrack n\rbrack=e^{j\omega n}$作為頻率響應的輸入訊號??

• 頻率響應的目的是希望可以分析DSP系統在各種不同頻率的訊號的工作狀況
• 其實最好的方式應該是灌入弦波訊號較為直觀，而且 $x\lbrack n\rbrack=e^{j\omega n}$ 在真實的世界中是不存在的
• 但是在數學推導上面如果輸入的是一個複指數訊號，可以將頻率響應簡化成對DSP系統的傅立葉轉換分析，整個過程會簡化很多
• 因此在訊號分析中探討頻率響應的系統特性都是假設灌入的訊號為複指數訊號
• 由上述的數學式可以導出當 $x\lbrack n\rbrack=e^{j\omega n}$，則 $y\lbrack n\rbrack=H(\omega)\cdot x\lbrack n\rbrack$，可以看出複指數訊號對於DSP系統是一個特徵函數 (Eigenfunction)，因此在DSP的理論中都喜歡用$x\lbrack n\rbrack=e^{j\omega n}$輸入系統分析系統的特性。

6.2 Moving Average Filter頻率響應 (數學式推導)

移動平均濾波器定義

$h\lbrack n\rbrack=\frac1M\{1,\;1,\;\dots,\;1\},\;n=0,\;1,\;\dots,\;(M-1)$

Moving Average Filter頻率響應推導

$
\begin{align}
H(\omega) &= \sum_{k=-\infty}^\infty h\lbrack k\rbrack\cdot e^{-j\omega k}\\
&= \frac1M\sum_{k=0}^{M-1}e^{-j\omega k}\\
&= \frac1M\left(\sum_{k=0}^\infty e^{-j\omega k}-\sum_{k=M}^\infty e^{-j\omega k}\right)\\
&= \frac1M\left(\sum_{k=0}^\infty e^{-j\omega k}\right)\cdot\left(1-e^{-j\omega M}\right)\\
&= \frac1M\cdot\frac{1-e^{-j\omega M}}{1-e^{-j\omega}}\\
&= \frac1M\cdot\frac{\sin\left(\frac{M\omega}2\right)}{\sin\left({\displaystyle\frac\omega2}\right)}\cdot e^\frac{-j(M-1)\omega}2
\end{align}
$

所以

Moving Average Filter的強度頻率響應為

$\left|H(\omega)\right|=\left|\frac{\sin\left(\frac{M\omega}2\right)}{\sin\left({\displaystyle\frac\omega2}\right)}\right|$

6.3 Moving Average Filter頻率響應 (Python Code)

如果依照上述推導Moving Average Filter的方式，一旦設計出了新的濾波器，就要從新推導一次頻率響應的數學式，太辛苦了，而且搞不好會有推導不出的可能性。

在Python的Scipy函式庫中已經有幫忙做出濾波器的頻率響應的Library可以使用。

• 程式碼

    import numpy as np
    import scipy.signal as signal
    import matplotlib.pyplot as plt
  
    # 輸入濾波器的大小
    filter_size = eval( input( "Please enter filter size: " ) )
    # 產生濾波器的係數
    h = np.ones( filter_size ) / filter_size
  
    # 產生濾波器頻率響應點
    # w: x座標，單位: *pi rad/sample，數值範圍: [0, pi)
    # H: y座標，頻率響應點，為複數，可計算出Magnitude與Phase
    w, H = signal.freqz( h )
    # 計算 Magnitude
    mag = abs( H )
    # 計算 Phase
    phase = np.angle( H )
  
    plt.plot( w, mag )
    plt.xlabel( r'$\omega$ (*pi rad/sample)' )
    plt.ylabel( 'Magnitude' )
    plt.title("Magnitude Response (Default 512 points)")
    plt.show( )
  
    plt.plot( w, phase )
    plt.xlabel( r'$\omega$ (*pi rad/sample)' )
    plt.ylabel( 'Phase (radians)' )
    plt.title("Phase Response (Default 512 points)")
    plt.show( )
  

上述程式碼主要的核心程式兩個

• 決定出濾波器的係數
  • 這部分在下一個章節 7. 濾波器設計 會提到
• 將濾波器係數透過 scipy.signal.freqz 函式庫就可以算出頻率響應的點

• M=5 輸出結果
  
  

• M=15 輸出結果
  
  

透過上述兩個結果可以得出幾個結論

• Moving Average濾波器是一個低頻率波器 (Low-Pass Filter)
• Moving Average濾波器長度影響低頻可通過的頻帶寬度
• Moving Average濾波器長度越長，高頻部可通過的頻帶Magnitude越低，但是ripple較多
• 濾波器伴隨著相位偏移，以Moving Average濾波器來說，當濾除的頻率越高，相位出現落後的現象也越嚴重

注意!!，在時域訊號處理中不存在完美濾波的濾波器，無法完全排除不需要的頻帶的訊號

頻率響應結果做圖說明

• 點數量的說明
  • 預設計算出的點是512點，而且只有正頻率 (不同於傅立葉分析會有正負頻率)
  • 要調整點數scipy.signal.freqz有參數可以修正，點的數量改變就是修改頻率響應圖的解析度
    • 回想一下， 5. 頻譜分析 章節也有提到這樣的觀念

• $x$座標單位說明 (很重要)

  • 函式庫算出的點 $x$座標單位都是 “$rad/sample$”，但這是甚麼意思呢?

  • 角頻率與 $Hz$ 的轉換公式為 $\omega=2\pi f$，把上述Moving Average濾波器 $M=15$ 的$x$座標全部除以 $2\pi$，可以得到下圖

    • 做圖一
      

  • 此時頻率響應圖的 $x$座標單位是 “$rad/sample$”，而且可以看到最大值只有 $0.5$ ???

  • 因為做頻率響應分析時需要參考數位訊號系統的取樣頻率 ($f_s$) ，但是在設計濾波器時可能不一定知道系統的$f_s$ 是多少，因此Python或是Matlab這類的工具在畫出頻率響應圖的時候 $x$座標單位給的是 "Normalized Frequency”。
    • 回想一下在 5. 頻譜分析 章節有提到這跟Nyquist-Shannon 取樣定理有關
  • Normalized Frequency公式 $f_n=\frac f{f_s}=\frac{\displaystyle\frac{rad}{sec}}{\displaystyle\frac{sample}{sec}}=\frac{rad}{sample}$，因此想要把頻率響應的 $x$座標改成較為直覺的 $Hz$ 做單位就是把$x$座標乘上 $\frac{f_s}{2\pi}$ (公式推導: 1. 訊號的描述 )
    • 做圖二 (假設 $f_s=1000\;Hz$)
      

6.4 頻率響應 v.s. 頻譜分析

Ref:

• 頻率響應
• Frequency Response of the Running Average Filter
• What is actually Normalized Frequency?
• [SciPy-user] normalized frequencies to pi radians / sample
• Normalized Frequency（x PI rad/sample）含义
• 数字角频率ω与模拟角频率Ω的理解
• 浅谈模拟频率，模拟角频率，数字角频率之间的关系
• 在离散傅里叶变换（DFT）中如何理解采样频率、模拟频率和数字频率三者之间的关系？

• 数字信号处理中，模拟角频率和数字角频率是什么关系?