查看“︁格策爾演算法”︁的源代码

'''格策爾演算法'''或'''格茲爾演算法'''( 英語:'''Goertzel algorithm''' )是[[数字信号处理|數位訊號處理]]的一種運算技巧，此運算技巧提供一個有效率的方式來估計部分區域的[[离散傅里叶变换|離散傅立葉轉換]]，廣泛的運用在數字電話中的的[[双音多频|雙音多頻信號]](每個撥號的數字鍵由兩個頻率的音所組成，一個低頻，一個高頻)，此演算法在1958年被{{tsl|en|Gerald_Goertzel|傑拉德 · 格策爾}}所提出<ref>Goertzel, G. (January 1958),&nbsp;[http://www.jstor.org/discover/10.2307/2310304?uid=3739664&uid=2134&uid=2&uid=70&uid=4&uid=3739256&sid=55978747833 "An Algorithm for the Evaluation of Finite Trigonometric Series"],&nbsp;''American Mathematical Monthly'',&nbsp;'''65'''&nbsp;(1): 34–35,&nbsp;[[Digital object identifier|doi]]:[https://doi.org/10.2307%2F2310304 10.2307/2310304]</ref>。

格策爾演算法與離散傅立葉轉換的相似處在於他們都可以分析某個特定頻段的離散訊號<ref>Mock, P. (March 21, 1985), "[http://www.ti.com/lit/an/spra168/spra168.pdf Add DTMF Generation and Decoding to DSP-μP Designs] {{Wayback|url=http://www.ti.com/lit/an/spra168/spra168.pdf |date=20201021042635 }}" (PDF), EDN, [https://www.worldcat.org/title/edn-with-eee-exclusively-for-designers-and-design-managers-in-electronics/oclc/310949239 ISSN 0012-7515] {{Wayback|url=https://www.worldcat.org/title/edn-with-eee-exclusively-for-designers-and-design-managers-in-electronics/oclc/310949239 |date=20181106004702 }}; also found in DSP Applications with the TMS320 Family, Vol. 1, Texas Instruments, 1989.</ref><ref>Chen, Chiouguey J. (June 1996),&nbsp;[http://focus.ti.com/lit/an/spra066/spra066.pdf ''Modified Goertzel Algorithm in DTMF Detection Using the TMS320C80 DSP''] (PDF), Application Report, Texas Instruments, SPRA066</ref><ref>Schmer, Gunter (May 2000),&nbsp;[http://focus.ti.com/lit/an/spra096a/spra096a.pdf ''DTMF Tone Generation and Detection: An Implementation Using the TMS320C54x''] {{Wayback|url=http://focus.ti.com/lit/an/spra096a/spra096a.pdf |date=20110604204843 }} (PDF), Application Report, Texas Instruments, SPRA096a</ref>；相反的，它們的不同處在於，格策爾演算法每次疊代的運算都是使用實數的乘法。雖然說在全頻域的計算上，格策爾演算法會比其他的傅立葉轉換快速演算法的複雜度來的高，但是它能區段式的分析每個小區段的頻率組成，因此可以編寫成較簡單的運算架構，實際應用在處理器內的數值計算會更有效率。

格策爾演算法逆向操作生成出弦波，而這個過程只需花費一個乘法和一個加法運算<ref><nowiki>http://haskell.cs.yale.edu/wp-content/uploads/2011/01/AudioProc-TR.pdf</nowiki>.</ref>。

== 演算法 ==
格策爾演算法把[[离散傅里叶变换|離散傅立葉轉換]]看成是一組濾波器，將輸入的訊號與[[濾波器]]中的[[冲激响应|脈衝響應]]做卷積運算，求得[[濾波器]]的輸出，即得到頻率域其中一點的頻率<math>X(k)</math>。此演算法利用[[旋轉因子]]<math>{\omega}^k_N</math>的週期性，將[[离散傅里叶变换|離散傅立葉轉換]]轉換為線性的濾波運算。

因為旋轉因子
{{NumBlk|:|<math>{\omega}^{-kN}_N=e^{-j(2{\pi/N})(-kN)}=1</math>|1}}
可得轉換後第k點的頻率為
{{NumBlk|:|<math>X(k)={\omega}^{-kN}_N\sum_{m=0}^{N-1}x(m){\omega}^{km}_N=\sum_{m=0}^{N-1}x(m){\omega}^{-k(N-m)}_N \qquad\quad ,k=0,1,2,...,N-1</math>|2}}
定義<math>y_{k}(n)</math>為
{{NumBlk|:|<math> y_{k}(n)=\sum_{m=0}^{N-1}x(m){\omega}^{-k(n-m)}_N </math>|3.A}}
可將<math>y_{k}(n)</math>理解為由兩個訊號的[[卷积|卷積]]運算得出的結果
{{NumBlk|:|<math> y_{k}(n)=x(n)\otimes h_k(n) </math>|3.B}}
其中<math>x(n)</math>式輸入的N點訊號，另外一個<math>h_k(n)</math>則被看作是[[无限冲激响应|IIR]]濾波器的[[冲激响应|脈衝頻率響應]]
{{NumBlk|:|<math> h_k(n)={\omega}^{-kn}_N \  u(n) </math>|4}}
對比(2)和(3)式，可推知(3.A)進行卷積運算，當n=N時，濾波器的輸出<math>y_{k}(N)</math>即為<math>X(k)</math>:
{{NumBlk|:|<math>X(k)=y_k(n)\lfloor_{n=N}</math>|5}}
對(4)進行Z轉換，可得一階[[无限冲激响应|IIR]][[传递函数|轉移函數]]
{{NumBlk|:|<math>H_k(z)=\frac {1}{1-{{\omega}^{-k} \ z^{-1}}}</math>|6}}
圖一為此系統的[[流程图|流程圖]]，其對應的[[遞迴關係式|差分方程式]]為:
{{NumBlk|:|<math>y_{k}(n)={\omega}^{-k}_N \ y_k(n-1)+x(n), \qquad \ y(-1)=0</math>|7}}
[[File:格策爾一階系統圖.PNG|缩略图|550x550像素|圖一、格策爾一階濾波器系統示意圖]]
依照此差分方程進行[[迭代法|疊代運算]]，疊代到<math>n=N </math>時即可依據(5)式得到<math>X(k)</math>。而依照轉移函數(6)式進行運算時，可以先將[[旋轉因子]]<math>{\omega}^k_N</math>儲存起來，每次疊代包含一次複數乘法，則按照(1)式計算N點[[离散傅里叶变换|離散傅立葉轉換]]時則需要<math>4N^2</math>次實數乘法運算和<math>N(4N-2)</math>次加/減法<ref>http://www.docin.com/p-577391532.html</ref>，加/減法與乘法運算皆為<math>4N^2</math>次，當N不大時運算效率不佳，若改為接下來改進的的格策爾演算法(二階)，所需的實數乘法次數約為原本的一半。

將式(6)上下同乘以<math>1-{\omega}^k_N \ z^{-1}</math>，可得第k點的頻率響應轉移函數為{{NumBlk|:|<math>\begin{alignat}{2} H_k(z) & = \frac {1-{{\omega}^{k}_N \ z^{-1}}}{(1-{{\omega}^{-k}_N \ z^{-1}})(1-{{\omega}^{k}_N \ z^{-1})}}
\\ & = \frac {1-{{\omega}^{k}_N \ z^{-1}}}{1-2 \ \cos((2{\pi}/N)k) z^{-1}+z^{-2})} \\ \end{alignat}</math>
|8}}
[[File:格策爾二階系統圖.PNG|缩略图|566x566像素|圖二、格策爾二階濾波器系統圖]]
此轉移函數所對應的系統流程圖如圖二所示，複數分析(8)式，可得知此二階濾波器有一對共軛的極點與一個零點。圖二中在計算<math>x(n)</math>的轉換結果<math>X(k)</math>時，會有兩個步驟:
# 共軛極點疊代計算   依序將輸入訊號<math>x(0),x(1),x(2),...,x(n-1)</math>放入濾波器做疊代運算，共作N次疊代，計算量是<math>2N</math>次實數乘法與<math>4N</math>次實數加/減法
# 零點疊代計算   輸入訊號<math>x(n)</math>是N點的訊號從<math>n=0,1,2,3,...,N-1</math>。加入<math>x(N)=0</math>的邊界條件，可以按照圖二的流程圖計算出<math>y_{k}(N)</math>，此即為所求的<math>x(n)</math>離散傅立葉轉換<math>X(k)</math>，此步驟的計算量為4次實數乘法與4次實數加/減法。
綜合以上步驟，總共的計算量為<math>2N+4</math>次實數乘法運算以及<math>4N+4</math>次實數加法運算，而使用此計算演算法只需儲存<math>{\omega}^k_N</math>與<math>\cos((2{\pi}/N)k)</math>兩個參數<ref>{{Cite web|url=http://en.dsplib.org/content/goertzel.html|title=格策爾介紹網站(英文)|accessdate=2017-06-22|author=|date=|publisher=|archive-url=https://web.archive.org/web/20170622095124/http://en.dsplib.org/content/goertzel.html|archive-date=2017-06-22|dead-url=yes}}</ref>。

== 相關條目 ==
* [[双音多频|雙音多頻]]
* [[Chirp-Z轉換]]
* [[頻率偏移調變]](FSK)
* [[相位偏移調變]](PSK)

== 參考資料 ==
<references />

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

== 外部連結 ==
* [https://web.archive.org/web/20170619103127/http://en.dsplib.org/content/goertzel.html 格策爾演算法網站介紹](英文)
* [http://www.embedded.com/design/configurable-systems/4006427/A-DSP-algorithm-for-frequency-analysis 頻域分析數位訊號數理演算法] {{Wayback|url=http://www.embedded.com/design/configurable-systems/4006427/A-DSP-algorithm-for-frequency-analysis |date=20170801100259 }}(英文)
* [http://www.embedded.com/print/4024443 格策爾演算法網站介紹] {{Wayback|url=http://www.embedded.com/print/4024443 |date=20170801035205 }}(英文) 作者: Kevin Bank,2002

{{Authority control}}
[[Category:快速傅立葉轉換演算法]]
[[Category:线性滤波器]]
[[Category:離散轉換]]