Goertzelアルゴリズム – 日本の辞書

Goertzel_algorithm
Goertzelアルゴリズムは、離散フーリエ変換（DFT）の個々の項を効率的に評価するためのデジタル信号処理（DSP）の手法です。これは、従来のアナログ電話のキーパッドのプッシュボタンによって生成されるデュアルトーン多重周波数信号（DTMF）トーンの認識など、特定の実用的なアプリケーションで役立ちます。このアルゴリズムは、1958年にGeraldGoertzelによって最初に記述されました。
DFTと同様に、Goertzelアルゴリズムは、離散信号から選択可能な1つの周波数成分を分析します。直接DFT計算とは異なり、Goertzelアルゴリズムは、実数値入力シーケンスに実数値演算を使用して、各反復で単一の実数値係数を適用します。スペクトル全体をカバーする場合、Goertzelアルゴリズムは、高速フーリエ変換（FFT）アルゴリズムよりも複雑度が高くなりますが、選択した少数の周波数成分を計算する場合は、数値的に効率的です。Goertzelアルゴリズムの単純な構造により、小型プロセッサや組み込みアプリケーションに最適です。
Goertzelアルゴリズムは、正弦波合成関数として「逆に」使用することもできます。これは、生成されたサンプルごとに1回の乗算と1回の減算のみを必要とします。

コンテンツ
1 アルゴリズム
2 数値安定性
3 DFT計算
4 アプリケーション
4.1 パワースペクトル用語 4.2 実数値演算を使用した単一のDFT項 4.3 位相検出 4.4 実際の算術演算における複雑な信号
5 計算の複雑さ
6 も参照してください
7 参考文献
8 参考文献
9 外部リンク

アルゴリズム
では、の
です。方程式とラベルの間の一貫性のないスタイル。可能であれば、このセクション
Goertzelアルゴリズムの主な計算は、デジタルフィルターの形式であるため、このアルゴリズムは、Goertzelフィルターと呼ばれることがよくフィルタは入力シーケンスで動作しますX
[ n ] { x }
$x""$
パラメータを使用して2段階のカスケードでω 0
{ omega _ {0}}

、分析する頻度を示し、サンプルあたりのラジアンに正規化します。
最初の段階では、中間シーケンスを計算します。 s [ n ] { s }

： s [ n] =X
[ n] + 2 cos（（ ω
0）。 s [ n− 1 ] − s
[ n− 2
] { s = x +2 cos（ omega _ {0}）s -s。}

（1）第2段階では、次のフィルターをに適用します s [ n ] { s }

、出力シーケンスを生成します y [ n ] { y }

： y [ n] = s
[ n] − e − j ω 0 s [ n− 1
] { y = s -e ^ {-j omega _ {0}}s。}

（2）最初のフィルターステージは、直接型構造の2次IIRフィルターであることがわかります。この特定の構造には、その内部状態変数がそのステージからの過去の出力値と等しいという特性が入力値X
[ n ] { x }

ためにn < 0
{ n <0}

すべて0に等しいと推定されます。評価をサンプルで開始できるように、初期フィルター状態を確立します。X
[ 0 ] { x }

、フィルターの状態には初期値が割り当てられます s [ −2 ] = s
[ −1 ] = 0
{ s = s = 0}

。エイリアシングの危険を回避するために、頻度ω 0
{ omega _ {0}}

多くの場合、0からπの範囲に制限されます（ナイキスト-シャノンのサンプリング定理を参照）。この範囲外の値を使用することは無意味ではありませんが、指数関数は2πの周期で周期的であるため、この範囲内のエイリアス周波数を使用することと同じです。ω 0
{ omega _ {0}}

。
第2ステージのフィルターは、計算で過去の出力を使用しないため、 FIRフィルターであることがわかります。
Z変換法は、フィルターカスケードの特性を研究するために適用できます。式（1）で与えられる最初のフィルターステージのZ変換は次のとおりです。 S （（ z
）。 X （（ z
）。 1 − 2 cos（（ ω 0 ）。z − 1 + z − 2= 1（（ 1− e + j ω 0 z− 1 ）。（（ 1− e − j ω 0 z− 1
）。 { { begin {aligned} { frac {S（z）} {X（z）}}＆= { frac {1} {1-2 cos（ omega _ {0}）z ^ { -1} + z ^ {-2}}} \＆= { frac {1} {（1-e ^ {+ j omega _ {0}} z ^ {-1}）（1-e ^ {-j omega _ {0}} z ^ {-1}）}}。 end {aligned}}}

（3）式（2）で与えられる2番目のフィルターステージのZ変換は次のようになります。 Y （（ z
）。 S （（ z
）。= 1 − e − j
ω0 − 1 { { frac {Y（z）} {S（z）}} = 1-e ^ {-j omega _ {0}} z^{-1}。}

（4）次に、2つのフィルターステージのカスケードの結合伝達関数は次のようになります。 S （（ z
）。 X （（ z
）。 Y （（ z
）。 S （（ z
）。= Y（（ z
）。 X （（ z
）。 = （（ 1− e − j ω 0 z− 1 ）。（（ 1− e + j ω 0 z− 1 ）。（（ 1− e − j ω 0 z− 1 ）。 =1 − e + j ω 0 z − 1 { { begin {aligned} { frac {S（z）} {X（z）}} { frac {Y（z）} {S（z）}} = { frac {Y（z） } {X（z）}}＆= { frac {（1-e ^ {-j omega _ {0}} z ^ {-1}）} {（1-e ^ {+ j omega _ { 0}} z ^ {-1}）（1-e ^ {-j omega _ {0}} z ^ {-1}）}} \＆= { frac {1} {1-e ^ { + j omega _ {0}} z^{-1}}}。end{aligned}}}

（5）これは、同等の時間領域シーケンスに変換して戻すことができ、用語はインデックスの最初の入力用語にロールバックされます。n = 0
{ n = 0}

： y [ n] =X
[ n] + e + j 0 y
[ n− 1 ] = ∑
k= − ∞ X
[ k] e + j 0（（ n− k
）。= e
j 0 ∑ k= 0 nX
[ k] e −
j 0
以来 ∀ k < 0 X
[ k] = 0。
{ { begin {aligned} y ＆= x + e ^ {+ j omega _ {0}} y \＆= sum _ {k =- infty} ^ {n} x e ^ {+ j omega _ {0}（nk）} \＆= e ^ {j omega _ {0} n} sum _ {k = 0} ^ {n} x e ^ {-j omega _ {0} k} qquad { text {since}} forall k <0、x =0。end{aligned}}}

（6）

数値安定性
フィルタのZ変換の極が次の場所にあることがわかります。e + j ω 0
{ e ^ {+ j omega _ {0}}}

とe −
jω 0
{ e ^ {-j omega _ {0}}}

、複素Z変換平面の原点を中心とする単位半径の円上。このプロパティは、低精度の算術演算と長い入力シーケンスを使用して計算した場合、フィルタープロセスがわずかに安定しており、数値エラーの累積に対して脆弱であることを示しています。数値的に安定したバージョンがChristianReinschによって提案されました。

DFT計算
DFT項を計算する重要なケースでは、次の特別な制限が適用されます。
フィルタリングはインデックスで終了しますn = N
{ n = N}

、どこ N { N}

DFTの入力シーケンスの項の数です。
Goertzel分析用に選択された周波数は、特別な形式に制限されていますω 0 = 2 π k
N{ omega _ {0} = 2 pi { frac{k}{N}}。}
（7）
インデックス番号 k { k}

DFTの「周波数ビン」がインデックス番号のセットから選択されていることを示しますk ∈ {{ 0 1 2
、。、N − 1
}{ k in {0，1，2、…、N-1}。}
（8）
これらの置換を式（6）に変換し、次の項を観察します。e + j 2
π 1
{ e ^ {+ j2 pi k} = 1}

、式（6）は、次の形式になります。 y [ N] = ∑ n = 0 N X [ n] e − j 2 π n k N { y = sum _ {n = 0} ^ {N} x e ^ {-j2 pi { frac{nk}{N}}}。}

（9）式（9）の右辺は、DFT項の定義式と非常に似ていることがわかります。X
[ k ] { X }

、インデックス番号のDFT用語 k { k}

、しかし完全に同じではありません。式（9）に示されている合計には、N + 1
{ N + 1}

用語を入力しますが、 N { N}

DFTを評価するときに、入力項を使用できます。単純ですがエレガントではない方法は、入力シーケンスを拡張することですX
[ n ] { x }

もう1つの人工的な価値を持つX
[ N] = 0
{ x = 0}

。式（9）から、最終結果に対する数学的効果は項を削除することと同じであることがわかります。X
[ N ] { x }

合計から、意図したDFT値を提供します。
ただし、余分なフィルターパスを回避するより洗練されたアプローチが式（1）から、拡張入力項がX
[ N] = 0
{ x = 0}

最終ステップで使用されます、 s [ N] = 2 cos（（ ω
0）。 s [ N− 1 ] − s
[ N− 2
] { s = 2 cos（ omega _ {0}）s -s。}

（10）したがって、アルゴリズムは次のように完了することができます。
入力項の処理後にIIRフィルターを終了しますX
[ N− 1 ]
{ x }

、
式（10）を適用して構築する s [ N ] { s }

以前の出力から s [ N− 1 ]
{ s }
と s
[ N− 2 ]
{ s }

、
計算された式（2）を適用します s [ N ] { s }

値と s [ N− 1 ]
{ s }

フィルタの最終的な直接計算によって生成されます。
最後の2つの数学演算は、代数的に組み合わせることで簡略化されます。 y [ N] = s
[ N] − e − j 2 πk N s
[ N− 1 ] =（（ 2 cos （（ ω 0 ）。 s [ N− 1 ] − s
[ N− 2 ]
）。− e − j 2 π kN s
[ N− 1 ] = e j 2π k N s
[ N− 1 ] − s
[ N− 2
] { { begin {aligned} y ＆= s -e ^ {-j2 pi { frac {k} {N}}} s \＆=（2 cos（ omega _ {0}）s -s ）-e ^ {-j2 pi { frac {k} {N}}} s ＆= e ^ {j2 pi { frac {k} {N}}} s -s。end{aligned}}}

（11）期間中にフィルターの更新を停止することに注意してくださいN − 1
{ N-1}

そして、式（11）ではなく式（2）をすぐに適用すると、最終的なフィルター状態の更新が失われ、誤った位相の結果が得られます。
Goertzelアルゴリズム用に選択された特定のフィルタリング構造は、その効率的なDFT計算の鍵です。1つの出力値のみが観察できます y [ N ] { y }

はDFTの計算に使用されるため、他のすべての出力項の計算は省略されます。FIRフィルターが計算されないため、IIRステージの計算 s [ 0 ] s [ 1 ] { s 、s }

、などは、第1ステージの内部状態を更新した直後に破棄できます。
これはパラドックスを残しているようです。アルゴリズムを完了するには、FIRフィルターステージをIIRフィルターステージからの最後の2つの出力を使用して一度評価する必要がありますが、計算効率のために、IIRフィルターの反復はその出力値を破棄します。これは、直接形式のフィルター構造のプロパティが適用される場所です。IIRフィルターの2つの内部状態変数は、IIRフィルター出力の最後の2つの値を提供します。これらは、FIRフィルターステージを評価するために必要な項です。

アプリケーション

パワースペクトル用語
式（6）を調べると、項を計算するための最後のIIRフィルターパス y [ N ] { y }

補足入力値を使用するX
[ N] = 0
{ x = 0}

前の項に大きさ1の複素乗数を適用します y [ N− 1 ]
{ y }

。その結果、 y [ N ] { y }

と y [ N− 1 ]
{ y }

同等の信号電力を表します。式（11）を適用し、項から信号電力を計算することも同様に有効です。 y [ N ] { y }

または式（2）を適用し、項から信号電力を計算します y [ N− 1 ]
{ y }

。どちらの場合も、DFT項で表される信号電力について次の式になります。X
[ k ] { X }

：X
[ k]X ′
[ k] = y
[ N] y ′
[ N] = y
[ N− 1 ] y ′
[ N− 1 ] = s 2
[ N− 1 ] + s 2
[ N− 2 ] − 2 cos（（ 2π k N
）。 s [ N− 1 ] s
[ N− 2
] { { begin {aligned} X 、X’＆= y 、y’ = y 、y’ ＆= s ^ {2} + s ^ {2} -2 cos left（2 pi { frac {k} {N}} right）、s 、s。end {aligned}}}

（12）以下の擬似コードでは、変数sprevとsprev2はIIRフィルターからの出力履歴を一時的に保存しますが、は入力を保存する配列xのインデックス付き要素です。 x
ここで定義されたNtermここで選択されたKtermω=2×π×Kterm/Nterms;cr：= cos（ω）ci：= sin（ω）係数：=2×crsprev：= 0sprev2：= 00からNterms-1の範囲の各インデックスnに対して s：= x +coeff×sprev–sprev2 sprev2：= sprev sprev：= s終わりパワー：=sprev2×sprev2+sprev×sprev-係数×sprev×sprev2
計算を整理して、入力サンプルが更新間でフィルター状態を維持するソフトウェアオブジェクトに単独で配信され、他の処理が行われた後に最終的な電力結果にアクセスできるようにすることができます。

実数値演算を使用した単一のDFT項
実数値の入力データの場合は、特に入力ストリームが物理プロセスの直接測定から生じる組み込みシステムで頻繁に発生します。前のセクションの図と比較すると、入力データが実数値の場合、フィルターの内部状態変数sprevとsprev2実数値も観察できるため、最初のIIRステージで複雑な演算は必要ありません。実数値演算の最適化は、通常、変数に適切な実数値データ型を適用するのと同じくらい簡単です。
入力項を使用して計算した後X
[ N− 1 ]
{ x }

、およびフィルターの反復が終了した場合、式（11）を適用してDFT項を評価する必要が最終的な計算では複素数値の算術を使用しますが、これは実数と虚数の項を分離することで実数値の算術に変換できます。c r= cos （（ 2π k N
）。 c I =
sin （（ 2π k N
）。 y
[ N] = c r s
[ N− 1 ] − s
[ N− 2 ] + j c I s [ N− 1
] { { begin {aligned} c_ {r}＆= cos（2 pi { tfrac {k} {N}}）、\ c_ {i}＆= sin（2 pi { tfrac {k} {N}}）、\ y ＆= c_ {r} s -s + jc_ {i}s。end{aligned }}}

（13）パワースペクトルアプリケーションと比較すると、唯一の違いは、終了に使用される計算です。（信号電力の実装と同じIIRフィルター計算）XKreal = sprev * cr –sprev2;XKimag = sprev * ci;

位相検出
このアプリケーションでは、DFT用語の同じ評価が必要ですX
[ k ] { X }

、前のセクションで説明したように、実数値または複素数値の入力ストリームを使用します。次に、信号位相は次のように評価できます。ϕ =
日焼け− 1 ℑ（（ X
[ k ] ）。 ℜ （（ X
[ k ] ）。 { phi = tan ^ {-1} { frac { Im（X ）} { Re（X ）}}、}

（14）逆正接関数を計算するときは、特異点、象限などに適切な予防策を講じます。

実際の算術演算における複雑な信号
複素数信号は実数部と虚数部に線形に分解されるため、Goertzelアルゴリズムは実数部のシーケンスに対して個別に実数演算で計算でき、次のようになります。y r n ]
{ y _ { text {r}} }

、および架空の部分のシーケンスを超えて、y I n ]
{ y _ { text {i}} }

。その後、2つの複素数値の部分的な結果を再結合できます。 y [ n] = y r n ] +j y I n
] { y = y _ { text {r}} + jy _ { text{i}}。}

（15）

計算の複雑さ
で
計算の複雑さの理論によると、 M { M}
M
{ M}

データセットへのGoertzelアルゴリズムの適用 N { N}

「操作あたりのコスト」が K { K}

複雑さがある O （（ KN M ）。 { O（KNM）}

。
単一の
DFTビン
を計算するにはX（（ f ）。 { X（f）}

長さの複雑な入力シーケンスの場合 N { N}

、Goertzelアルゴリズムには2 N
{ 2N}

掛け算と4 N
{ 4 N}

ループ内の加算/減算、および4つの乗算と4つの最終的な加算/減算、合計2 N + 4
{ 2N + 4}

掛け算と4 N + 4
{ 4N + 4}

足し算/引き算。これは、それぞれに対して繰り返されます M { M}

周波数。
対照的に、データセットでFFTを使用する N { N}

値には複雑さがあります O （（ K N ログ N ）。
{ O（KN log N）}

。
これは、使用するFFTアルゴリズムに依存するため、直接適用するのは困難ですが、典型的な例は、基数2のFFTです。 2 ログ 2 （（ N ）。 { 2 log _ {2}（N）}

掛け算と 3 ログ 2 （（ N ）。 { 3 log _ {2}（N）}

DFTビン
ごとの加算/減算 N
{ N}

ビン。
複雑度の順序式で、計算された項の数が M { M}

より小さい
ログ N { log N}

、Goertzelアルゴリズムの利点は明らかです。ただし、FFTコードは比較的複雑であるため、「作業単位あたりのコスト」の要素 K { K}

多くの場合、FFTの方が大きく、実用上の利点は、 M { M}

の数倍
ログ 2 （（ N ）。 { log _ {2}（N）}

。
基数2FFTまたはGoertzelアルゴリズムのどちらがより効率的かを判断するための経験則として、項の数を調整します N { N}

データセットで、最も近い正確な2の累乗まで上向きに、これを呼び出しますN 2
{ N_ {2}}

、およびGoertzelアルゴリズムは、次の場合により高速になる可能性がM ≤ 5 N2 N
ログ 2 （（ N 2 ）。
{ M leq { frac {5N_ {2}} {6N}} log _ {2}（N_ {2}）}
FFTの実装と処理プラットフォームは、相対的なパフォーマンスに大きな影響を与えます。一部のFFT実装は、内部複素数計算を実行してオンザフライで係数を生成し、「作業単位あたりのコストK」を大幅に増加させます。FFTおよびDFTアルゴリズムは、数値効率を高めるために事前に計算された係数値のテーブルを使用できますが、これには外部メモリにバッファリングされた係数値へのアクセスが多く必要であり、数値の利点の一部に対抗するキャッシュ競合の増加につながる可能性が
どちらのアルゴリズムも、複素数値ではなく実数値の入力データを使用すると、約2倍の効率が得られます。ただし、これらのゲインはGoertzelアルゴリズムでは当然ですが、特定のアルゴリズムバリアントを使用しないとFFTでは達成されません実数値データの変換に特化しています。

も参照してください
ブルースタインのFFTアルゴリズム（chirp-Z）
周波数偏移変調（FSK）
位相偏移変調（PSK）

参考文献
^ Goertzel、G.（1958年1月）、「有限三角級数の評価のためのアルゴリズム」、American Mathematical Monthly、65（1）：34–35、doi：10.2307 / 2310304、JSTOR 2310304
^ Mock、P.（1985年3月21日）、「DSP-μP設計へのDTMF生成とデコードの追加」（PDF）、EDN、ISSN 0012-7515 ; TMS320ファミリを使用したDSPアプリケーション、Vol。1、テキサスインスツルメンツ、1989年。
^ Chen、Chiouguey J.（1996年6月）、TMS320C80 DSPを使用したDTMF検出の修正Goertzelアルゴリズム（PDF）、アプリケーションレポート、Texas Instruments、SPRA066
^ Schmer、Gunter、DTMFトーンの生成と検出：TMS320C54xを使用した実装（PDF）、アプリケーションレポート、Texas Instruments、SPRA096a
^ チェン、エリック; Hudak、Paul、Haskellでのオーディオ処理とサウンド合成（ PDF）、2017-03-28にオリジナル（PDF）からアーカイブ
^ 紳士、WM（1969年2月1日）。「フーリエ係数を計算するためのGoertzel（ワット）の方法のエラー分析」。コンピュータジャーナル。12（2）：160–164。土井：10.1093 / comjnl/12.2.160。
^ Stoer、J .; Bulirsch、R.（2002）、数値解析入門、Springer、ISBN 9780387954523 ^ 「Goertzelのアルゴリズム」。Cnx.org。2006-09-12 。
^ 「ElectronicEngineeringTimes|グローバルエレクトロニクスコミュニティの接続」。EETimes 。
^ Elmenreich、Wilfried（2011年8月25日）。「Goertzelフィルターを使用して周波数を効率的に検出する」。
^ 押す; フラナリー; Teukolsky; Vetterling（2007）、「Chapter 12」、数値レシピ、The Art of Scientific Computing、ケンブリッジ大学出版局

参考文献
プロアキス、JG; Manolakis、DG（1996）、Digital Signal Processing：Principles、Algorithms、and Applications、Upper Saddle River、NJ：Prentice Hall、pp。480–481、Bibcode：1996dspp.book ….. P

外部リンク
ウェイバックマシンでのGoertzelアルゴリズム（2018年6月28日アーカイブ）
周波数分析のためのDSPアルゴリズム
KevinBanksによるGoertzelアルゴリズム”