格策爾演算法

格策爾演算法把離散傅立葉轉換看成是一組濾波器，將輸入的訊號與濾波器中的脈衝響應做卷積運算，求得濾波器的輸出，即得到頻率域其中一點的頻率${\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)

1. Goertzel, G. (January 1958), "An Algorithm for the Evaluation of Finite Trigonometric Series", American Mathematical Monthly, 65 (1): 34–35, doi:10.2307/2310304
2. Mock, P. (March 21, 1985), "Add DTMF Generation and Decoding to DSP-μP Designs （页面存档备份，存于）" (PDF), EDN, ISSN 0012-7515 （页面存档备份，存于）; also found in DSP Applications with the TMS320 Family, Vol. 1, Texas Instruments, 1989.
3. Chen, Chiouguey J. (June 1996), Modified Goertzel Algorithm in DTMF Detection Using the TMS320C80 DSP (PDF), Application Report, Texas Instruments, SPRA066
4. Schmer, Gunter (May 2000), DTMF Tone Generation and Detection: An Implementation Using the TMS320C54x （页面存档备份，存于） (PDF), Application Report, Texas Instruments, SPRA096a
5. http://haskell.cs.yale.edu/wp-content/uploads/2011/01/AudioProc-TR.pdf.
6. http://www.docin.com/p-577391532.html
7. . [2017-06-22]. （原始内容存档于2017-06-22）.

• Proakis, J. G.; Manolakis, D. G. (1996), Digital Signal Processing: Principles, Algorithms, and Applications, Upper Saddle River, NJ: Prentice Hall, pp. 480–481

• 格策爾演算法網站介紹(英文)
• 頻域分析數位訊號數理演算法 （页面存档备份，存于）(英文)
• 格策爾演算法網站介紹 （页面存档备份，存于）(英文) 作者: Kevin Bank,2002