格策爾演算法

格策爾演算法或格茲爾演算法( 英語:Goertzel algorithm )是數位訊號處理的一種運算技巧，此運算技巧提供一個有效率的方式來估計部分區域的離散傅立葉轉換，廣泛的運用在數字電話中的的雙音多頻信號(每個撥號的數字鍵由兩個頻率的音所組成，一個低頻，一個高頻)，此演算法在1958年被傑拉德 · 格策爾（英语：Gerald_Goertzel）所提出^[1]。

格策爾演算法與離散傅立葉轉換的相似處在於他們都可以分析某個特定頻段的離散訊號^[2][3][4]；相反的，它們的不同處在於，格策爾演算法每次疊代的運算都是使用實數的乘法。雖然說在全頻域的計算上，格策爾演算法會比其他的傅立葉轉換快速演算法的複雜度來的高，但是它能區段式的分析每個小區段的頻率組成，因此可以編寫成較簡單的運算架構，實際應用在處理器內的數值計算會更有效率。

格策爾演算法逆向操作生成出弦波，而這個過程只需花費一個乘法和一個加法運算^[5]。

演算法

格策爾演算法把離散傅立葉轉換看成是一組濾波器，將輸入的訊號與濾波器中的脈衝響應做卷積運算，求得濾波器的輸出，即得到頻率域其中一點的頻率${\displaystyle X(k)}$。此演算法利用旋轉因子${\displaystyle {\omega }_{N}^{k}}$的週期性，將離散傅立葉轉換轉換為線性的濾波運算。

因為旋轉因子

${\displaystyle {\omega }_{N}^{-kN}=e^{-j(2{\pi /N})(-kN)}=1}$

1

可得轉換後第k點的頻率為

${\displaystyle X(k)={\omega }_{N}^{-kN}\sum _{m=0}^{N-1}x(m){\omega }_{N}^{km}=\sum _{m=0}^{N-1}x(m){\omega }_{N}^{-k(N-m)}\qquad \quad ,k=0,1,2,...,N-1}$

2

定義${\displaystyle y_{k}(n)}$為

${\displaystyle y_{k}(n)=\sum _{m=0}^{N-1}x(m){\omega }_{N}^{-k(n-m)}}$

3.A

可將${\displaystyle y_{k}(n)}$理解為由兩個訊號的卷積運算得出的結果

${\displaystyle y_{k}(n)=x(n)\otimes h_{k}(n)}$

3.B

其中${\displaystyle x(n)}$式輸入的N點訊號，另外一個${\displaystyle h_{k}(n)}$則被看作是IIR濾波器的脈衝頻率響應

${\displaystyle h_{k}(n)={\omega }_{N}^{-kn}\ u(n)}$

4

對比(2)和(3)式，可推知(3.A)進行卷積運算，當n=N時，濾波器的輸出${\displaystyle y_{k}(N)}$即為${\displaystyle X(k)}$:

${\displaystyle X(k)=y_{k}(n)\lfloor _{n=N}}$

5

對(4)進行Z轉換，可得一階IIR轉移函數

${\displaystyle H_{k}(z)={\frac {1}{1-{{\omega }^{-k}\ z^{-1}}}}}$

6

圖一為此系統的流程圖，其對應的差分方程式為:

${\displaystyle y_{k}(n)={\omega }_{N}^{-k}\ y_{k}(n-1)+x(n),\qquad \ y(-1)=0}$

7

圖一、格策爾一階濾波器系統示意圖

依照此差分方程進行疊代運算，疊代到${\displaystyle n=N}$時即可依據(5)式得到${\displaystyle X(k)}$。而依照轉移函數(6)式進行運算時，可以先將旋轉因子${\displaystyle {\omega }_{N}^{k}}$儲存起來，每次疊代包含一次複數乘法，則按照(1)式計算N點離散傅立葉轉換時則需要${\displaystyle 4N^{2}}$次實數乘法運算和${\displaystyle N(4N-2)}$次加/減法^[6]，加/減法與乘法運算皆為${\displaystyle 4N^{2}}$次，當N不大時運算效率不佳，若改為接下來改進的的格策爾演算法(二階)，所需的實數乘法次數約為原本的一半。

將式(6)上下同乘以${\displaystyle 1-{\omega }_{N}^{k}\ z^{-1}}$，可得第k點的頻率響應轉移函數為

${\displaystyle {\begin{alignedat}{2}H_{k}(z)&={\frac {1-{{\omega }_{N}^{k}\ z^{-1}}}{(1-{{\omega }_{N}^{-k}\ z^{-1}})(1-{{\omega }_{N}^{k}\ z^{-1})}}}\\&={\frac {1-{{\omega }_{N}^{k}\ z^{-1}}}{1-2\ \cos((2{\pi }/N)k)z^{-1}+z^{-2})}}\\\end{alignedat}}}$

8

圖二、格策爾二階濾波器系統圖

此轉移函數所對應的系統流程圖如圖二所示，複數分析(8)式，可得知此二階濾波器有一對共軛的極點與一個零點。圖二中在計算${\displaystyle x(n)}$的轉換結果${\displaystyle X(k)}$時，會有兩個步驟:

1. 共軛極點疊代計算 依序將輸入訊號${\displaystyle x(0),x(1),x(2),...,x(n-1)}$放入濾波器做疊代運算，共作N次疊代，計算量是${\displaystyle 2N}$次實數乘法與${\displaystyle 4N}$次實數加/減法
2. 零點疊代計算 輸入訊號${\displaystyle x(n)}$是N點的訊號從${\displaystyle n=0,1,2,3,...,N-1}$。加入${\displaystyle x(N)=0}$的邊界條件，可以按照圖二的流程圖計算出${\displaystyle y_{k}(N)}$，此即為所求的${\displaystyle x(n)}$離散傅立葉轉換${\displaystyle X(k)}$，此步驟的計算量為4次實數乘法與4次實數加/減法。

綜合以上步驟，總共的計算量為${\displaystyle 2N+4}$次實數乘法運算以及${\displaystyle 4N+4}$次實數加法運算，而使用此計算演算法只需儲存${\displaystyle {\omega }_{N}^{k}}$與${\displaystyle \cos((2{\pi }/N)k)}$兩個參數^[7]。

相關條目

• 雙音多頻
• Chirp-Z轉換
• 頻率偏移調變(FSK)
• 相位偏移調變(PSK)