← ブログ一覧

Goertzel アルゴリズム：1つの周波数だけを見たいとき

📡 信号・RF

FFT を知ってから、周波数解析といえばとりあえず FFT、という感じになっていた。でも少し前に「特定の1周波数成分だけを N サンプルごとに更新したい」という状況になって、FFT を使うのが明らかに過剰だと気づいた。そこで見つけたのが Goertzel アルゴリズムだった。

ロックインアンプを知っている人は、ひとまず「あれの1周波数・有限ブロック版」だと思ってもらっていい。関係はあとで整理する。

問題設定：電話のプッシュ音（DTMF）

例として、電話のプッシュボタン音、いわゆる DTMF（Dual-Tone Multi-Frequency）を考える。名前のとおり、1つのボタンで「2つの音（デュアルトーン）」を同時に鳴らす方式だ。

キーパッドを 4 行 × 4 列の格子だと思ってほしい。行に1つずつ「低い側」の周波数（697, 770, 852, 941 Hz）、列に1つずつ「高い側」の周波数（1209, 1336, 1477, 1633 Hz）を割り当てる。ボタンを押すと、その行の低音と、その列の高音が、1つずつペアで鳴る。

1209 Hz1336 Hz1477 Hz1633 Hz
697 Hz123A
770 Hz456B
852 Hz789C
941 Hz*0#D

たとえば「1」は左上だから 697 Hz ＋ 1209 Hz、「5」なら 770 Hz ＋ 1336 Hz、という具合だ。行が4本・列が4本なので、出てくる周波数は全部で 8 種類しかない。この 8 つだけを見張って「いま鳴っている低音1つ＋高音1つ」の組を読めば、4 × 4 = 16 通りのボタンが一意に決まる。

つまり検出側の仕事は、8 つの決まった周波数それぞれについて「今この成分のエネルギーはあるか？」を判定することになる。

N = 205 サンプル（8 kHz サンプリング）を集めて判定するとしよう。ここで FFT を使うと、205 点ぶんの複素スペクトルを全部計算することになる。欲しいのは 8 本だけなのに、だ。これが「過剰」の正体だ。

DFT の k 番目だけを取り出す

DFT の定義はこうだ：

X[k]=∑n=0N−1x[n]⋅e−j2πkn/NX[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j2\pi kn/N}

数式だけ見ると硬いけど、やっていることは素朴だ。各サンプル x[n]x[n] を角度 −2πkn/N-2\pi k n/N だけ回してから、全部足し合わせる。狙った周波数 kk と入力の周波数・位相がそろっているほど、回したベクトルが同じ向きに重なって、合計が長く伸びる。絵にするとこうなる：

DFTのk番目は、各サンプルを少しずつ回して継ぎ足した合計ベクトル

WN=ej2π/NW_N = e^{j2\pi/N} と略記すると X[k]=∑n=0N−1x[n] WN−knX[k] = \sum_{n=0}^{N-1} x[n] \, W_N^{-kn} と書ける。

ここで小さな事実を1つ使う：WN−kN=e−j2πk=1W_N^{-kN} = e^{-j2\pi k} = 1（kk は整数）。値が変わらないので全体に掛けておくと、指数の中身を N−nN-n の形に書き換えられる：

X[k]=WN−kN∑n=0N−1x[n] WN−kn=∑n=0N−1x[n] WN−k(N−n)X[k] = W_N^{-kN} \sum_{n=0}^{N-1} x[n] \, W_N^{-kn} = \sum_{n=0}^{N-1} x[n] \, W_N^{-k(N-n)}

この形が効いてくる。WN−k(N−n)W_N^{-k(N-n)} は「残り N−nN-n ステップぶん回す」という意味だから、結局こういう操作と同じになる——サンプルを1つ受け取るたびに、それまでの合計に WN−kW_N^{-k} を1回だけ掛けて、新しいサンプルを足す。漸化式で書くと：

y[n]=x[n]+WN−k y[n−1],y[−1]=0y[n] = x[n] + W_N^{-k} \, y[n-1], \quad y[-1] = 0

これを N 回まわして、最後の出力だけ読めば X[k]=y[N−1]X[k] = y[N-1] になる。

（きちんと言えば、上の和は「x[n]x[n] と WN−knW_N^{-kn} の畳み込みを n=Nn=N で評価したもの」と一致する、という話なのだけど、そこは難しく考えなくていい。要は「走っている合計に、毎回ちょっと回転をかけながら足していくだけ」だ。）

ただし、この WN−kW_N^{-k} は複素数だ。だから毎サンプル複素数の掛け算が要る＝ちょっと重い。ここを軽くするのが次の一手だ。

実数演算に落とす

さっきの漸化式は、極が単位円上にある1次フィルターになっている。WN−k=e−j2πk/NW_N^{-k} = e^{-j2\pi k/N} は大きさ1・角度 −2πk/N-2\pi k/N の複素数だから、その極（共鳴点）はちょうど狙った周波数 kk の真上に乗る。言い換えると、これは周波数 kk にぴったり同調した共振器で、その成分が入ってくると振動が育っていく。

問題は係数が複素数なこと。毎サンプル複素乗算（実数の掛け算4回ぶん）はマイコンには重い。そこで Goertzel のトリックが効く——この極を、その鏡像（複素共役 e+j2πk/Ne^{+j2\pi k/N}）とペアにする。単位円上の角度 ±2πk/N\pm 2\pi k/N にある2つの極を組み合わせると、係数が実数の2次フィルターになるのだ：

単位円上の共役な2つの極と、それに対応する実係数の2次フィルター

s[n]=x[n]+2cos⁡ ⁣(2πkN)s[n−1]−s[n−2]s[n] = x[n] + 2\cos\!\left(\frac{2\pi k}{N}\right) s[n-1] - s[n-2]

係数は 2cos⁡(2πk/N)2\cos(2\pi k/N) という実数1つと、−1-1 だけ。複素数がループから消えた。初期値 s[−1]=s[−2]=0s[-1] = s[-2] = 0 で N サンプル回したあと、最後の2つの状態 s[N−1]s[N-1]（最新）と s[N−2]s[N-2] から、複素の出力を1回だけ組み立てる：

X[k]=s[N−1]−e−j2πk/N⋅s[N−2]X[k] = s[N-1] - e^{-j2\pi k/N} \cdot s[N-2]

実部・虚部に展開すると：

Re[X[k]]=s[N−1]−cos⁡ ⁣(2πkN)s[N−2]\mathrm{Re}[X[k]] = s[N-1] - \cos\!\left(\tfrac{2\pi k}{N}\right) s[N-2] Im[X[k]]=sin⁡ ⁣(2πkN)s[N−2]\mathrm{Im}[X[k]] = \sin\!\left(\tfrac{2\pi k}{N}\right) s[N-2]

ループ内は実数の加減乗算だけ。複素演算はループが終わった後に1回かかるだけだ。位相が要らず、エネルギー（パワー）だけ欲しいなら、それすら要らない：

∣X[k]∣2=s[N−1]2+s[N−2]2−2cos⁡ ⁣(2πkN)s[N−1]⋅s[N−2]|X[k]|^2 = s[N-1]^2 + s[N-2]^2 - 2\cos\!\left(\tfrac{2\pi k}{N}\right) s[N-1] \cdot s[N-2]

DTMF 検出のように「その周波数が鳴っているか」を見るだけなら、このパワーで十分だ。

計算量の比較——ただし万能ではない

手法演算量（乗算数オーダー）
直接 DFT（全 N 点）O(N2)O(N^2)
FFT（基数2）O(Nlog⁡2N)O(N \log_2 N)
Goertzel（K 周波数）O(NK)O(NK)

K 個の周波数だけ必要なとき、Goertzel が FFT より速くなる条件は：

NK<Nlog⁡2N  ⟹  K<log⁡2NNK < N \log_2 N \implies K < \log_2 N

N = 1024 なら K < 10、N = 512 なら K < 9。DTMF の 8 周波数は、ちょうどこの境界あたりにいる。

ここで「思ったより速くないな、それなら FFT でよくない？」と感じたなら、その勘は正しい。Goertzel は「いつでも FFT より速い魔法」ではない。狙う周波数が増えて KK が log⁡2N\log_2 N に近づくと優位はすぐ消えるし、スペクトル全体が欲しいなら素直に FFT のほうが速い。

Goertzel が本当に効くのは、演算量の表に出てこない別の事情があるときだ：

だから「ごく少数の決まった周波数を、小さなハードで、ストリームのまま見張りたい」という DTMF のような場面が Goertzel の主戦場で、汎用のスペクトル解析はやっぱり FFT、という素直な棲み分けになる。

実装してみると

Python で書くとこんな感じだ：

import math

def goertzel(samples, k, N):
    coeff = 2 * math.cos(2 * math.pi * k / N)
    s1, s2 = 0.0, 0.0           # s1=s[n-1], s2=s[n-2]
    for x in samples:
        s0 = x + coeff * s1 - s2
        s2, s1 = s1, s0
    # パワー |X[k]|² だけ返す（複素演算なし）
    return s1**2 + s2**2 - coeff * s1 * s2

ループの中身は加算2回と乗算2回。マイコンの DSP 命令にそのまま乗せやすい形をしている。

最初の DTMF に戻ると、8 周波数ぶんこれを回して、低音グループと高音グループからいちばん強い1本ずつを選べば、押されたボタンが分かる：

LOW  = [697, 770, 852, 941]
HIGH = [1209, 1336, 1477, 1633]
KEYS = [['1', '2', '3', 'A'],
        ['4', '5', '6', 'B'],
        ['7', '8', '9', 'C'],
        ['*', '0', '#', 'D']]

def decode_dtmf(samples, fs=8000):
    N = len(samples)
    def power_at(f):
        k = round(f / fs * N)        # その周波数に一番近い DFT ビン
        return goertzel(samples, k, N)
    row = max(range(4), key=lambda i: power_at(LOW[i]))
    col = max(range(4), key=lambda j: power_at(HIGH[j]))
    return KEYS[row][col]

8 × N サンプルを舐めるだけ。FFT のように全ビンを計算しないぶん、こういう用途では素直に軽い。

ちょっとした歴史：FFT より古い

このアルゴリズムを発表したのは Gerald Goertzel で、1958 年の論文だ。面白いのは、これが Cooley–Tukey の FFT（1965 年）より7年も前だということ。当時はまだ「DSP（ディジタル信号処理）」という分野が立ち上がる前で、Goertzel の動機も信号処理というより、三角級数（フーリエ係数）を計算機で効率よく評価するという数値計算寄りの工夫だった。

それが後年、電話交換機の DTMF 検出という「少数の周波数を安いハードで見張る」用途にぴたりとはまって、DSP の定番として生き残った。生まれたときの目的と、有名になった理由がずれているアルゴリズム——というのが、個人的には少し好きだ。

ロックインと何が違う？

冒頭で「ロックインの1周波数版」と書いた。その答え合わせをしておく。

ロックインアンプ（別記事で書いた）は、信号にリファレンスの sin⁡/cos⁡\sin/\cos を掛けて積分する＝狙った周波数だけを直流に“折り畳んで”取り出す相関検出器だった。Goertzel がやっている ∑x[n] e−j2πkn/N\sum x[n]\, e^{-j2\pi kn/N} も、中身はまったく同じ「リファレンスとの相関」だ。狙った周波数の波と内積を取っている。

ではどこが“お得”かというと、リファレンスの正弦波を毎サンプル作らなくていいところ。素直なロックイン（や直接 DFT）は、サンプルごとに sin⁡\sin と cos⁡\cos の値を生成して掛ける。Goertzel は2次の漸化式がその回転を肩代わりしてくれるので、内側のループに要るのは定数 2cos⁡(2πk/N)2\cos(2\pi k/N) ひとつだけ。三角関数テーブルも複素乗算も、ループから消える。

ざっくり言うと——ロックインは連続（アナログでも回る）の相関検出器、Goertzel はそれを有限 N サンプルのブロックで、いちばん安い算術に落とした離散版、という関係だ。狙っているものは同じで、Goertzel は「実装が一番軽い形」を選んでいる、ということになるね。

何が面白いか

FFT はすべての周波数を一括で求める「グローバルな」変換だ。それに対して Goertzel は「特定の周波数を狙い撃ちにする局所的なフィルター」と捉えられる。同じ DFT を、全部まとめて解くか・狙った1本だけ引くか、という見方の違いだ。

もう一つ：kk を整数に限る必要はない。kk を実数にすれば、DFT グリッドの外の任意の周波数も計算できる。これは「ズーム FFT（Zoom FFT）」の発想にも繋がる話で、スペクトルの特定帯域だけを細かく見たいときに使える。

さらに、中間状態 s[n]s[n] を初期化せずに走らせ続ければ、ブロックを待たずスライディング窓でリアルタイム更新する拡張（Sliding DFT）にもなる。

小さなアルゴリズムなのに、意外と奥がある。

出典

— ランキン