固定電話のボタンを押すと「ピポパ」という音が鳴ります。あの音には、押されたボタンを相手の交換機が聞き分けるための情報が、2つの音の組み合わせとして埋め込まれています。このとき交換機がやりたいのは「全周波数のスペクトルを描くこと」ではなく、「あらかじめ決まった8個の周波数のうち、どれとどれが鳴っているか」だけを判定することです。スペクトル全体を計算するのは明らかに過剰で、もったいない計算です。
このように「特定の数個の周波数の強さだけを知りたい」場面で、FFT（高速フーリエ変換）よりも少ない計算量で答えを出すのがゲーツェルアルゴリズム（Goertzel algorithm）です。離散フーリエ変換（DFT）のうち、欲しい1本のビン $X[k]$ だけを2次のIIRフィルタとして計算する巧妙な方法で、メモリも演算量も小さいため、マイコンやDSPなど資源の限られた組込み機器で広く使われています。
ゲーツェルアルゴリズムが活きる応用先は数多くあります。
- DTMF復号: 電話のプッシュ音（Dual-Tone Multi-Frequency）を、低演算量のマイコンで8周波だけ監視して鍵を判定する
- トーン検出・FSK受信: モデムやテレメトリで特定周波数の有無を検出する、産業用の制御トーン検出
- 楽器のチューナー / 周波数推定: 注目する音階の周波数成分の強さだけを安価に測る
本記事では、DFTの単一ビンの式から出発して、それをz変換し、共役な極を持つ2次IIRフィルタの形に変形する過程を一行ずつ追います。そのうえで「N点・1周波数あたり」の演算量を数え、FFTより少なくなる条件を明らかにし、最後にPythonでゲーツェルフィルタ・複数トーンのパワー検出・DTMF復号器・FFTとの演算量比較を実装します。
本記事の内容
- DFTの単一ビン $X[k]$ の意味と、なぜ全ビンを計算したくないのか
- ゲーツェルアルゴリズムの2次IIRフィルタとしての導出（z変換・共役極の整理を省略せず）
- パワー（マグニチュード二乗）を最終1ステップだけで求める式の導出
- N点・1周波数あたりの演算量と、FFTより有利になる条件
- Pythonでのゲーツェルフィルタ実装、複数トーンのパワー検出、DTMF復号器、演算量比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ゲーツェルアルゴリズムとは
スペクトルアナライザを思い浮かべてください。横軸に周波数、縦軸に強さが並ぶグラフを描くには、すべての周波数ビンを計算する必要があります。FFTはこれを $O(N\log N)$ という驚くほど速い計算量で実現する魔法のような道具です。しかし、もし知りたいのが「440Hzの音は鳴っているか？」という1点だけだったら、グラフ全体を描くのは無駄な仕事です。聞きたいのは1音だけなのに、オーケストラ全体を採譜しているようなものです。
ゲーツェルアルゴリズムは、この「1点だけ知りたい」という要求に正面から応えます。アイデアの核心は、DFTの定義式 $X[k]=\sum_n x[n]e^{-j2\pi kn/N}$ を「フィルタへの入力を1サンプルずつ流し込む再帰計算」として書き直すことです。フィルタとして実装すると、係数は $\cos(2\pi k/N)$ という実数1個だけで済み、しかも各サンプルでの計算が乗算1回・加算2回という非常に軽い操作になります。サンプルを流し終わった最後の1回だけ、複素数の仕上げ計算をして $X[k]$ を取り出します。
直感的には、ゲーツェルフィルタは「周波数 $\omega_k = 2\pi k/N$ にぴったり共鳴する振り子」です。入力信号にその周波数の成分が含まれていれば、振り子はだんだん大きく揺れ、含まれていなければ揺れません。N個のサンプルを流し終えたときの揺れの大きさが、その周波数の成分の強さ（DFT係数）に対応します。共鳴という物理的なイメージを持っておくと、後で出てくる「単位円上の極」という話も腑に落ちやすくなります。
ここで自然な疑問が生まれます。DFTの和の式が、どうして「共鳴する振り子（IIRフィルタ）」に化けるのでしょうか。その変形こそがゲーツェルアルゴリズムの本体です。次のセクションで、DFTの単一ビンの式をz変換し、フィルタの形へと一行ずつ書き換えていきます。
DFTの単一ビンからフィルタへ
ゴールの宣言
これから示すことは次の一点です。DFTの $k$ 番目のビン
$$ X[k] = \sum_{n=0}^{N-1} x[n] \, e^{-j 2\pi k n / N} $$
を、係数 $\cos(2\pi k/N)$ をもつ実数の2次IIRフィルタを $N$ サンプル走らせるだけで計算できる、ということです。フィルタ部分は乗算が実数1個分しか要らず、複素数の演算は最後の1回だけに押し込めます。これがFFTより軽くなる源です。
畳み込みの形に書き換える
まず、DFTの定義式を少し書き換えます。指数の符号を反転させたい（フィルタの因果的な形に合わせたい）ので、$X[k]$ に $e^{-j 2\pi k N / N} = e^{-j2\pi k} = 1$ を掛けても値は変わらないことを使います。$k$ は整数なので $e^{-j2\pi k}=1$ です。これを使うと
$$ X[k] = e^{-j 2\pi k N / N} \sum_{n=0}^{N-1} x[n] \, e^{-j 2\pi k n / N} = \sum_{n=0}^{N-1} x[n] \, e^{j 2\pi k (N – n) / N} $$
と書けます。指数の符号がプラスに変わり、肩には $N-n$ が現れました。ここで右辺は「ある列 $x[n]$ と $e^{j2\pi k m/N}$ の畳み込みを $m=N$ で評価したもの」の形をしています。これを正確にするため、次のような数列 $y_k[m]$ を導入します。
$$ y_k[m] = \sum_{n=-\infty}^{\infty} x[n] \, h_k[m-n], \qquad h_k[m] = e^{j 2\pi k m / N} \, u[m] $$
ここで $u[m]$ は単位ステップ関数（$m \ge 0$ で1、それ以外で0）です。$x[n]$ が $n=0,\dots,N-1$ の範囲だけ非ゼロだとすると、$m=N$ における畳み込みは
$$ y_k[N] = \sum_{n=0}^{N-1} x[n] \, e^{j 2\pi k (N – n) / N} = X[k] $$
となり、ちょうど $X[k]$ に一致します。つまり$X[k]$ は、インパルス応答 $h_k[m]=e^{j2\pi km/N}u[m]$ を持つフィルタに $x[n]$ を入れたときの、時刻 $m=N$ における出力だと分かりました。あとはこのフィルタを「再帰計算」で安く実装できればよいわけです。
ここまでで「DFTの1本のビンはフィルタの1サンプル出力である」という見方が手に入りました。次に、このフィルタの伝達関数 $H_k(z)$ を求め、それを実数係数の都合のよい形に整えます。
複素1次フィルタの伝達関数
インパルス応答 $h_k[m] = e^{j 2\pi k m / N} u[m]$ をz変換します。$W_N = e^{j 2\pi k / N}$ と略記すると、これは公比 $W_N$ の等比数列のz変換です。
$$ H_k(z) = \sum_{m=0}^{\infty} W_N^{\,m} z^{-m} = \sum_{m=0}^{\infty} \left( W_N z^{-1} \right)^{m} = \frac{1}{1 – W_N z^{-1}} $$
最後の等号は $|W_N z^{-1}| < 1$ で収束する等比級数の和の公式 $\sum r^m = 1/(1-r)$ を使いました。これは極が $z = W_N = e^{j 2\pi k/N}$ にある複素係数の1次フィルタです。極が単位円のちょうど上（$|W_N|=1$）にあることに注目してください。これが「ぴったり共鳴する振り子」の正体で、周波数 $\omega_k = 2\pi k/N$ の成分を際限なく積み上げる仕組みになっています。
しかし、係数 $W_N$ が複素数のままだと、毎サンプルで複素乗算が必要になり、せっかくの軽さが台無しです。そこで次のひと工夫として、極の複素共役 $W_N^{*}=e^{-j2\pi k/N}$ をわざと分子・分母に掛けて、分母を実数化します。
共役極を導入して分母を実数化する
$H_k(z)$ の分子・分母に $(1 – W_N^{*} z^{-1})$ を掛けます。
$$ H_k(z) = \frac{1}{1 – W_N z^{-1}} \cdot \frac{1 – W_N^{*} z^{-1}}{1 – W_N^{*} z^{-1}} = \frac{1 – W_N^{*} z^{-1}}{(1 – W_N z^{-1})(1 – W_N^{*} z^{-1})} $$
分母を展開します。$W_N + W_N^{*} = 2\cos(2\pi k/N)$ と $W_N W_N^{*} = |W_N|^2 = 1$ を使うと
$$ (1 – W_N z^{-1})(1 – W_N^{*} z^{-1}) = 1 – (W_N + W_N^{*}) z^{-1} + W_N W_N^{*} z^{-2} = 1 – 2\cos\!\Big(\tfrac{2\pi k}{N}\Big) z^{-1} + z^{-2} $$
となり、分母が完全に実数の2次多項式になりました。これがゲーツェルの肝です。分子はまだ複素数 $1 – W_N^{*}z^{-1}$ ですが、これは最後にまとめて処理します。$\cos\theta_k$ を $c_k = \cos(2\pi k/N)$ と置いて、伝達関数を
$$ H_k(z) = \frac{1 – W_N^{*} z^{-1}}{1 – 2 c_k \, z^{-1} + z^{-2}} $$
と書きます。ここで重要なのは、分母（フィードバック部）が実数係数 $2c_k$ だけで決まっているという点です。フィルタを「再帰部（分母）」と「出力部（分子）」に分けて2段で実装すれば、サンプルごとに走るのは実数の再帰部だけにできます。
ここまでで「実数係数の2次の再帰」という形が見えてきました。次に、この伝達関数を実際の差分方程式に落とし、毎サンプルの更新式 $s[n]=x[n]+2c_k s[n-1]-s[n-2]$ を導きます。
中間変数 $s[n]$ の漸化式
伝達関数 $H_k(z) = Y_k(z)/X(z)$ を、中間信号 $S(z)$ を介して2段に分けます。
$$ \frac{Y_k(z)}{X(z)} = \underbrace{\frac{1}{1 – 2 c_k z^{-1} + z^{-2}}}_{\text{再帰部 } S(z)/X(z)} \cdot \underbrace{(1 – W_N^{*} z^{-1})}_{\text{出力部 } Y_k(z)/S(z)} $$
まず再帰部 $S(z)/X(z) = 1/(1 – 2c_k z^{-1} + z^{-2})$ を取り出します。これは
$$ S(z)\,(1 – 2 c_k z^{-1} + z^{-2}) = X(z) $$
すなわち $S(z) = X(z) + 2 c_k z^{-1} S(z) – z^{-2} S(z)$ を意味します。z領域の $z^{-1}$ は時間領域の1サンプル遅延なので、逆z変換すると次の漸化式になります。
$$ \boxed{\; s[n] = x[n] + 2 \cos\!\Big(\tfrac{2\pi k}{N}\Big)\, s[n-1] – s[n-2] \;} $$
これがゲーツェルアルゴリズムの中心となる更新式です。初期条件は $s[-1]=s[-2]=0$ とします。注目すべきは、この1行に複素数が一切現れないことです。1サンプルあたり、乗算は $2c_k \cdot s[n-1]$ の1回（係数 $2c_k$ は前計算しておく）、加減算は2回だけです。$n=0,1,\dots,N-1$ の $N$ サンプル分この漸化式を回し、$s[N-1]$ と $s[N-2]$ を保持しておけば、あとは1回の仕上げ計算で $X[k]$ が得られます。
ここまでで毎サンプルの再帰が手に入りました。残るは出力部 $Y_k(z)/S(z) = 1 – W_N^{*}z^{-1}$ を使って、最後に $X[k]$ を取り出す式です。次のセクションでそれを導きます。
最終ステップで $X[k]$ を取り出す
出力部は $Y_k(z) = (1 – W_N^{*} z^{-1}) S(z)$ なので、逆z変換すると
$$ y_k[n] = s[n] – W_N^{*} \, s[n-1] = s[n] – e^{-j 2\pi k/N} s[n-1] $$
求めたいのは $X[k] = y_k[N]$ です。漸化式を $n=0$ から $n=N$ まで回したいところですが、実装上は $n=0,\dots,N-1$ の $N$ サンプルだけ再帰部を回し、最後に $s[N]$ を「入力ゼロでもう1ステップ進めた値」として扱うのが定石です。入力が尽きた後（$x[N]=0$）の漸化式は $s[N] = 2c_k\,s[N-1] – s[N-2]$ なので、
$$ X[k] = y_k[N] = s[N] – e^{-j2\pi k/N} s[N-1] = \big(2 c_k s[N-1] – s[N-2]\big) – e^{-j 2\pi k/N} s[N-1] $$
ここで $2c_k = e^{j2\pi k/N} + e^{-j2\pi k/N}$ を代入して整理すると、$e^{-j2\pi k/N}s[N-1]$ の項が一部相殺され、
$$ X[k] = e^{j 2\pi k/N} s[N-1] – s[N-2] $$
という実装に便利な複素数の式が得られます。実部・虚部に分けて書くと
$$ \operatorname{Re} X[k] = \cos\!\Big(\tfrac{2\pi k}{N}\Big) s[N-1] – s[N-2], \qquad \operatorname{Im} X[k] = \sin\!\Big(\tfrac{2\pi k}{N}\Big) s[N-1] $$
となります。複素乗算が登場するのはこの最後の1回だけで、$N$ 回の再帰中はずっと実数だけで済んでいることが、改めて確認できます。
ところで、DTMF検出のように「位相は要らず、その周波数のエネルギー（強さ）だけ知りたい」場合は、複素数の式すら省けます。次にそのパワーの式を導きます。
パワー（マグニチュード二乗）の直接計算
検出だけが目的なら $|X[k]|^2$ が分かれば十分です。$X[k] = e^{j2\pi k/N}s[N-1] – s[N-2]$ の絶対値の二乗を、$|a-b|^2 = |a|^2 – 2\operatorname{Re}(a b^{*}) + |b|^2$ の公式で展開します。$|e^{j2\pi k/N}s[N-1]|^2 = s[N-1]^2$、$|s[N-2]|^2 = s[N-2]^2$、そしてクロス項は
$$ 2\operatorname{Re}\!\big(e^{j2\pi k/N} s[N-1] \cdot s[N-2]\big) = 2\cos\!\Big(\tfrac{2\pi k}{N}\Big) s[N-1] s[N-2] $$
です（$s$ は実数なので共役は不要）。これらをまとめると
$$ \boxed{\; |X[k]|^2 = s[N-1]^2 + s[N-2]^2 – 2\cos\!\Big(\tfrac{2\pi k}{N}\Big) s[N-1]\, s[N-2] \;} $$
という美しい式が得られます。$N$ 回の再帰を回し終えた後、最後に保持していた2つの値 $s[N-1], s[N-2]$ と前計算済みの係数 $2c_k$ から、乗算3回・加減算2回だけでパワーが求まります。三角関数も複素数もループの中には一切現れません。これがゲーツェルが組込みで好まれる理由です。
ここまでで、毎サンプルの軽い再帰と最後の仕上げで $X[k]$ やそのパワーが得られることが分かりました。では、これは本当にFFTより少ない計算量なのでしょうか。次のセクションで演算量を数えて比較します。
演算量の比較とFFTより有利になる条件
ゲーツェルの演算量
ゲーツェルアルゴリズムで1つの周波数ビンを求めるコストを数えます。再帰部 $s[n]=x[n]+2c_k s[n-1]-s[n-2]$ は、係数 $2c_k$ を前計算しておけば、1サンプルあたり実数乗算1回・実数加減算2回です。これを $N$ サンプル繰り返すので、再帰全体で実数乗算 $N$ 回、実数加減算 $2N$ 回。最後の仕上げ（パワー計算）は乗算数回・加算数回の定数コストです。したがって1ビンあたりおよそ
$$ \text{ゲーツェル（1ビン）} \approx N \text{ 実数乗算} + 2N \text{ 実数加減算} $$
です。$M$ 個の周波数を検出したいなら、これを $M$ 回繰り返して $\approx MN$ 乗算となります。重要なのは、欲しい周波数の数 $M$ に比例し、$N$ には線形でしか効かない点です。
DFT直接計算とFFTの演算量
比較対象を確認します。DFTを定義どおり計算すると、1ビンあたり $N$ 回の複素乗算が必要で、全 $N$ ビンで $O(N^2)$ です。FFTはこれを劇的に改善し、全 $N$ ビンを $\approx \tfrac{N}{2}\log_2 N$ 回の複素乗算（実数乗算ではおよそ $2N\log_2 N$ 程度）で求めます。FFTは「全ビンを一度に出す」点で圧倒的に効率的ですが、裏を返すと、欲しいビンが少数でも常に全ビン分を計算してしまいます。
交差点：いくつまでの周波数ならゲーツェルが勝つか
$M$ 個のビンだけ欲しいとき、ゲーツェルのコストは $\approx MN$（実数乗算）、FFTは $\approx 2N\log_2 N$（実数乗算、全ビン）です。ゲーツェルが有利になる条件は
$$ M N < 2 N \log_2 N \quad \Longleftrightarrow \quad M < 2 \log_2 N $$
です。両辺を $N$ で割れる（$N>0$）ので $N$ が消え、判定は欲しいビン数 $M$ と $2\log_2 N$ の大小だけで決まります。たとえば $N=205$（後で出てくるDTMFの典型値）なら $2\log_2 205 \approx 15.4$ なので、欲しい周波数が15個くらいまでならゲーツェルのほうが乗算回数が少なくなります。DTMFは検出すべき周波数がわずか8個ですから、まさにゲーツェルの独壇場です。
加えて、ゲーツェルはFFTのような大きな配列バッファ（$N$ 点の複素配列やビット反転テーブル）を必要とせず、ビンごとに状態変数 $s[n-1], s[n-2]$ の2つを持つだけで済みます。メモリが乏しいマイコンでは、この省メモリ性が演算量以上に効いてきます。これらが、ゲーツェルが組込み用途で選ばれる実務的な理由です。
ここまでで理論と演算量の根拠がそろいました。次は実際にPythonで実装し、単一トーンの検出、複数トーンへの応答、そしてDTMF復号までを動かして確かめます。
Pythonでの実装
ゲーツェルフィルタの実装
まず、導出した漸化式とパワー式をそのままコードに落とします。検出したい物理周波数 $f_\text{target}$ とサンプリング周波数 $f_s$、サンプル数 $N$ を与え、最寄りのビン番号 $k$ を選びます。
import numpy as np
def goertzel_power(x, fs, f_target):
"""ゲーツェルアルゴリズムで f_target のパワー |X[k]|^2 を返す"""
N = len(x)
# 目標周波数に最も近いビン番号 k（実数でよい）
k = N * f_target / fs
omega = 2.0 * np.pi * k / N
coeff = 2.0 * np.cos(omega) # 再帰係数 2cos(2πk/N) を前計算
s_prev = 0.0 # s[n-1]
s_prev2 = 0.0 # s[n-2]
# 再帰部: 毎サンプル 乗算1回 + 加減算2回（複素数なし）
for n in range(N):
s = x[n] + coeff * s_prev - s_prev2
s_prev2 = s_prev
s_prev = s
# 仕上げ: パワー = s[N-1]^2 + s[N-2]^2 - 2cos(ω) s[N-1] s[N-2]
power = s_prev**2 + s_prev2**2 - coeff * s_prev * s_prev2
return power
def goertzel_complex(x, fs, f_target):
"""ゲーツェルで複素 DFT 係数 X[k] = e^{jω} s[N-1] - s[N-2] を返す"""
N = len(x)
k = N * f_target / fs
omega = 2.0 * np.pi * k / N
coeff = 2.0 * np.cos(omega)
s_prev = 0.0
s_prev2 = 0.0
for n in range(N):
s = x[n] + coeff * s_prev - s_prev2
s_prev2 = s_prev
s_prev = s
# X[k] = e^{jω} s[N-1] - s[N-2]
return np.exp(1j * omega) * s_prev - s_prev2
このコードでは、coeff に $2\cos(2\pi k/N)$ を一度だけ計算して入れ、ループ内では x[n] + coeff*s_prev - s_prev2 という乗算1回・加減算2回だけを回しています。導出どおり、ループ内に複素数も三角関数も登場しないことが確認できます。パワー版は最後に1行で $|X[k]|^2$ を返し、複素版は仕上げの複素式で $X[k]$ そのものを返します。
NumPyのFFTとの一致確認
導出が正しいかどうかは、同じ周波数ビンに対してFFTの結果とゲーツェルの結果が一致するかで確かめられます。整数ビン $k$ を選んで両者を比べてみます。
import numpy as np
# 検証用信号: ちょうど整数ビンに乗る周波数を合成
fs = 8000 # サンプリング周波数 [Hz]
N = 200 # サンプル数
n = np.arange(N)
k_test = 25 # 検証するビン番号
f_test = k_test * fs / N # このビンに対応する物理周波数 [Hz]
x = 1.3 * np.cos(2*np.pi*f_test*n/fs + 0.4)
# FFT による X[k]
X_fft = np.fft.fft(x)[k_test]
# ゲーツェルによる X[k]
X_goertzel = goertzel_complex(x, fs, f_test)
print(f"f_test = {f_test:.1f} Hz (bin k={k_test})")
print(f"FFT X[k] = {X_fft.real:+.4f} {X_fft.imag:+.4f}j")
print(f"Goertzel X[k] = {X_goertzel.real:+.4f} {X_goertzel.imag:+.4f}j")
print(f"差の絶対値 = {abs(X_fft - X_goertzel):.3e}")
実行すると、FFTとゲーツェルの $X[k]$ は実部・虚部とも一致し、差の絶対値は $10^{-10}$ 程度（浮動小数点の丸め誤差レベル）になります。これは、ゲーツェルが近似ではなくDFTの厳密な単一ビンを別の順序で計算しているだけであることを意味します。位相情報まで含めて完全に一致する点は、振幅だけ合っていればよいと思っていた読者には嬉しい確認でしょう。
整数ビンでぴったり一致することが分かったので、次に「目標周波数の周りでパワーがどう立ち上がるか」を周波数掃引で可視化し、ゲーツェルが共鳴フィルタとして振る舞う様子を見てみます。
単一トーン検出と周波数応答の可視化
ゲーツェルフィルタが特定周波数にどれだけ鋭く反応するかを見るため、検出周波数を固定し、入力トーンの周波数を掃引してパワーを測ります。
import numpy as np
import matplotlib.pyplot as plt
fs = 8000
N = 205
n = np.arange(N)
f_detect = 1209.0 # 検出したい周波数 [Hz]（DTMF列の1つ）
# 入力トーンの周波数を掃引してゲーツェルパワーを測る
freqs = np.linspace(900, 1500, 400)
powers = []
for f in freqs:
tone = np.cos(2*np.pi*f*n/fs)
powers.append(goertzel_power(tone, fs, f_detect))
powers = np.array(powers)
powers_db = 10*np.log10(powers / powers.max() + 1e-12)
plt.figure(figsize=(9, 5))
plt.plot(freqs, powers_db, 'b-', linewidth=1.8)
plt.axvline(f_detect, color='red', linestyle='--', alpha=0.7,
label=f'検出周波数 {f_detect:.0f} Hz')
plt.xlabel('入力トーン周波数 [Hz]')
plt.ylabel('正規化ゲーツェルパワー [dB]')
plt.title(f'Goertzel filter response (N={N}, fs={fs} Hz)')
plt.ylim([-40, 2])
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig('goertzel_response.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから2つのことが読み取れます。第一に、ゲーツェルパワーは検出周波数 1209Hz でちょうどピークを持ち、そこから離れるにつれて急速に減衰します。共鳴する振り子という直感どおり、狙った周波数だけに強く反応するバンドパスフィルタとして働いていることが分かります。第二に、ピークの両側にサイドローブ（小さな盛り上がり）が見えます。これは有限長 $N$ の信号を切り出したことによる窓（矩形窓）の効果で、DFTのスペクトル漏れと同じ現象です。$N$ を大きくすればピークは鋭くなり、周波数分解能 $f_s/N$ が細かくなります。
単一トーンへの鋭い応答が確認できました。実際のDTMFでは2つの周波数が同時に鳴ります。次に、複数のトーンが混ざった信号に対し、8個のゲーツェル検出器を並べてそれぞれのパワーを測ってみます。
複数トーン信号へのパワー検出
DTMFで使われる8つの周波数すべてに対してゲーツェル検出器を用意し、ある2周波が混ざった信号を入れて、どの検出器が強く反応するかを棒グラフで見ます。
import numpy as np
import matplotlib.pyplot as plt
# DTMF の 8 周波数 [Hz]（低群4 + 高群4）
DTMF_LOW = [697, 770, 852, 941]
DTMF_HIGH = [1209, 1336, 1477, 1633]
ALL_FREQS = DTMF_LOW + DTMF_HIGH
fs = 8000
N = 205
n = np.arange(N)
# 合成信号: 770 Hz と 1336 Hz の2トーン（= キー '5'）+ 雑音
np.random.seed(0)
signal2 = (np.cos(2*np.pi*770*n/fs) + np.cos(2*np.pi*1336*n/fs)
+ 0.2*np.random.randn(N))
# 8 周波それぞれのゲーツェルパワー
powers = [goertzel_power(signal2, fs, f) for f in ALL_FREQS]
powers = np.array(powers)
plt.figure(figsize=(9, 5))
colors = ['steelblue']*4 + ['darkorange']*4
bars = plt.bar([str(f) for f in ALL_FREQS], powers, color=colors)
plt.xlabel('検出周波数 [Hz] (青=低群, 橙=高群)')
plt.ylabel('ゲーツェルパワー')
plt.title('8-frequency Goertzel detection (tones: 770 Hz + 1336 Hz)')
plt.grid(True, axis='y', alpha=0.3)
plt.tight_layout()
plt.savefig('goertzel_dtmf_bars.png', dpi=150, bbox_inches='tight')
plt.show()
棒グラフを見ると、低群では770Hzの棒だけが、高群では1336Hzの棒だけが突出して高く、残り6本は雑音レベルに沈んでいます。雑音を加えたにもかかわらず、目標の2周波が明確に分離されています。これは、ゲーツェルが各周波数を独立した狭帯域フィルタで測っているため、他の周波数や広帯域雑音の影響を受けにくいからです。低群・高群それぞれで「最も強い1本」を選べば、押されたボタンが特定できそうだと直感的に分かります。次に、この判定ロジックを組み込んだDTMF復号器を実装します。
DTMF復号器の実装
DTMFは低群4周波・高群4周波の組み合わせ（$4\times4=16$ 通り）で、電話のキーパッド（0-9、*、#、A-D）を表します。低群と高群からそれぞれ最強の1本を選び、その組み合わせから鍵を引きます。誤検出を避けるため、「最強の周波数が次に強い周波数より十分大きいか」という簡単な妥当性チェックも入れます。
import numpy as np
DTMF_LOW = [697, 770, 852, 941]
DTMF_HIGH = [1209, 1336, 1477, 1633]
# (低群インデックス, 高群インデックス) -> キー
DTMF_KEYS = [
['1', '2', '3', 'A'],
['4', '5', '6', 'B'],
['7', '8', '9', 'C'],
['*', '0', '#', 'D'],
]
def decode_dtmf(x, fs, ratio_thresh=4.0):
"""1ブロックの信号からDTMFキーを1文字復号する。判定不能なら None"""
low_p = np.array([goertzel_power(x, fs, f) for f in DTMF_LOW])
high_p = np.array([goertzel_power(x, fs, f) for f in DTMF_HIGH])
i_low = int(np.argmax(low_p))
i_high = int(np.argmax(high_p))
# 妥当性チェック: 最強が2番目より ratio_thresh 倍以上強いか
def dominant(p, idx):
second = np.max(np.delete(p, idx))
return p[idx] > ratio_thresh * (second + 1e-12)
if dominant(low_p, i_low) and dominant(high_p, i_high):
return DTMF_KEYS[i_low][i_high]
return None
この関数は、まず8個のゲーツェルパワーを求め、低群・高群それぞれの最強ビンのインデックスを取り、DTMF_KEYS の表からキーを引きます。dominant は「最強が2番目の ratio_thresh 倍以上か」を確かめる単純なガードで、無音や雑音だけのブロックで誤って鍵を出さないようにする役割です。導出したパワー式が、この実用的な復号ロジックの土台になっていることに注目してください。
復号器ができたので、次に「123A」というキー列を実際にDTMF信号として合成し、ブロックごとに切り出して正しく復号できるかを確かめます。
キー列の合成と連続復号
複数キーを順に鳴らした信号を作り、各キーのブロックを切り出して decode_dtmf にかけ、元のキー列を復元できるかを検証します。
import numpy as np
import matplotlib.pyplot as plt
DTMF_LOW = [697, 770, 852, 941]
DTMF_HIGH = [1209, 1336, 1477, 1633]
KEY_FREQS = { # キー -> (低群, 高群)
'1': (697, 1209), '2': (697, 1336), '3': (697, 1477),
'4': (770, 1209), '5': (770, 1336), '6': (770, 1477),
'7': (852, 1209), '8': (852, 1336), '9': (852, 1477),
'0': (941, 1336), 'A': (697, 1633), 'B': (770, 1633),
}
fs = 8000
tone_dur = 0.05 # 各トーン長 [s]
N_tone = int(fs*tone_dur)
def make_tone(key, fs, N):
f_lo, f_hi = KEY_FREQS[key]
n = np.arange(N)
return np.cos(2*np.pi*f_lo*n/fs) + np.cos(2*np.pi*f_hi*n/fs)
sequence = "123A"
# 各キーをトーン → 無音(ギャップ) の順に連結
signal_parts = []
for key in sequence:
signal_parts.append(make_tone(key, fs, N_tone))
signal_parts.append(np.zeros(int(fs*0.02))) # 20ms の無音
full_signal = np.concatenate(signal_parts) + 0.1*np.random.randn(
sum(len(p) for p in signal_parts))
# ブロックに区切って各トーン区間を復号
block = 205
decoded = []
pos = 0
while pos + block <= len(full_signal):
seg = full_signal[pos:pos+block]
if np.mean(seg**2) > 0.2: # エネルギーがある区間だけ判定
key = decode_dtmf(seg, fs)
if key and (not decoded or decoded[-1][0] != key
or pos - decoded[-1][1] > block):
decoded.append((key, pos))
pos += block
print("送信キー列:", sequence)
print("復号キー列:", "".join(k for k, _ in decoded))
実行すると、復号キー列が送信した “123A” と一致します。ブロック長205サンプル（約25.6ms）で切り出し、エネルギーのある区間だけを判定し、同じキーの連続検出を1回にまとめる簡単な後処理を入れることで、雑音入りの連続信号からも正しく鍵列を取り出せています。これは、各キーが「低群1本・高群1本」というスパースな周波数構造を持ち、ゲーツェルがその少数の周波数だけをピンポイントで測れるからこそ実現できる処理です。次に、この「少数周波数だけ測る」ことの計算量的なメリットを、FFTと数値で比較してみます。
FFTとの演算量比較
理論で導いた「$M < 2\log_2 N$ ならゲーツェルが有利」という条件を、実際の乗算回数をモデル化して可視化します。横軸を欲しいビン数 $M$、縦軸を概算実数乗算回数として、ゲーツェルとFFTを比べます。
import numpy as np
import matplotlib.pyplot as plt
N = 205 # DTMF の典型ブロック長
M = np.arange(1, 33) # 欲しいビン数
goertzel_mults = M * N # ゲーツェル: 約 M*N 実数乗算
fft_mults = 2 * N * np.log2(N) # FFT: 全ビンで約 2N log2 N（M に依存しない）
fft_line = np.full_like(M, fft_mults, dtype=float)
cross = 2*np.log2(N) # 交差点 M = 2 log2 N
plt.figure(figsize=(9, 5))
plt.plot(M, goertzel_mults, 'o-', color='steelblue', label='Goertzel (~M·N)')
plt.plot(M, fft_line, 's-', color='darkorange',
label=f'FFT (~2N·log2 N, all bins)')
plt.axvline(cross, color='red', linestyle='--', alpha=0.7,
label=f'交差点 M=2log2 N ≈ {cross:.1f}')
plt.scatter([8], [8*N], color='green', zorder=5, s=80,
label='DTMF (M=8)')
plt.xlabel('検出したい周波数ビン数 M')
plt.ylabel('概算 実数乗算回数')
plt.title(f'Goertzel vs FFT computation (N={N})')
plt.grid(True, alpha=0.3)
plt.legend(fontsize=9)
plt.tight_layout()
plt.savefig('goertzel_vs_fft.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"交差点 M = 2 log2({N}) = {cross:.2f}")
print(f"DTMF (M=8): Goertzel {8*N} 乗算 vs FFT {fft_mults:.0f} 乗算")
このグラフから、ゲーツェルの乗算回数は $M$ に比例して直線的に増え、FFTは $M$ によらず一定（全ビン計算）であることが見て取れます。両者は $M = 2\log_2 N \approx 15.4$ 付近で交差し、それより少ないビン数ならゲーツェルが、多いビン数ならFFTが有利になります。緑の点で示したDTMF（$M=8$）は交差点の左側にあり、ゲーツェルのほうが乗算回数が少ない領域に明確に入っています。出力された数値でも、DTMFではゲーツェルがFFTのおよそ半分程度の乗算で済むことが確認できます。さらに前述のとおり、ゲーツェルは大きな配列バッファを要さない省メモリ性も併せ持つため、組込み機器での優位性は乗算回数の差以上に大きくなります。
これらの実装と比較により、ゲーツェルアルゴリズムが「少数の周波数を安く・省メモリで検出する」道具として実用的であることが、理論と数値の両面から確かめられました。
まとめ
本記事では、ゲーツェルアルゴリズムによる単一周波数検出とDTMF復号について解説しました。
- DFTの単一ビンはフィルタの1サンプル出力: $X[k]=\sum_n x[n]e^{-j2\pi kn/N}$ は、インパルス応答 $e^{j2\pi km/N}u[m]$ のフィルタに $x[n]$ を通したときの時刻 $N$ の出力 $y_k[N]$ に等しい
- 共役極で実数化: 複素1次フィルタ $1/(1-W_N z^{-1})$ の分子・分母に共役 $1-W_N^{*}z^{-1}$ を掛けると、分母が実数の2次多項式 $1-2\cos(2\pi k/N)z^{-1}+z^{-2}$ になる
- 再帰式: $s[n]=x[n]+2\cos(2\pi k/N)s[n-1]-s[n-2]$ を $N$ 回回す。ループ内は実数乗算1回・加減算2回で、複素数も三角関数も現れない
- パワーの式: 検出だけなら $|X[k]|^2=s[N-1]^2+s[N-2]^2-2\cos(2\pi k/N)s[N-1]s[N-2]$ を最後に1回計算するだけでよい
- 演算量: 1ビンあたり $\approx N$ 実数乗算。$M$ ビン欲しいときFFT（$\approx 2N\log_2 N$）より有利になる条件は $M<2\log_2 N$ で、$N$ に依存しない
- DTMF復号: 検出周波数はわずか8個なので交差点の十分内側に入り、ゲーツェルが省メモリ・低演算量で鍵を判定できる
ゲーツェルアルゴリズムは「DFTの一部だけを賢く計算する」という発想の好例であり、フィルタ理論・z変換・DFTの3つを橋渡しする教材としても優れています。共鳴フィルタとしての見方を押さえておくと、トーン検出やFSK復調、チューナーなど多くの応用にそのまま応用できます。
次のステップとして、以下の記事も参考にしてください。