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

問題設定

電話のプッシュボタン音（DTMF）を思い浮かべてほしい。「1」を押すと 697 Hz と 1209 Hz の正弦波が同時に鳴る。検出側がやるべきことは、8 つの周波数（697, 770, 852, 941, 1209, 1336, 1477, 1633 Hz）それぞれに「今この周波数のエネルギーはあるか？」を判定すること。

N = 205 サンプル（8 kHz サンプリング）を集めて判定するとして、FFT を使うと 205 点すべての複素スペクトルを計算することになる。8 周波数しか必要ないのに。

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

DFT の定義は：

$X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j2\pi kn/N}$

$W_N = e^{j2\pi/N}$ と書けば $X[k] = \sum_{n=0}^{N-1} x[n] \, W_N^{-kn}$。

ここで $W_N^{-kN} = e^{-j2\pi k} = 1$（k は整数）という事実に注目する。これを乗じても値は変わらない：

$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)}$

右辺は「$x[n]$ と $W_N^{-kn}$ の畳み込みを $n = N$ で評価したもの」に見える。だから次のフィルターを走らせて、最後の出力だけ読み出せばいい：

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

$X[k] = y[N-1]$ になる。ただしこれは複素係数なので、実装が少し重い。

実数演算に落とす

$W_N^{-k} = e^{-j2\pi k/N}$ の極は単位円上にある。この1次フィルターは安定しているが、複素乗算を含む。Goertzel のトリックは、これを共役極と組み合わせて実係数の2次フィルターにすること。

$s[n] = x[n] + 2\cos\!\left(\frac{2\pi k}{N}\right) s[n-1] - s[n-2]$

初期値 $s[-1] = s[-2] = 0$。N サンプル処理した後：

$X[k] = s[N] - e^{-j2\pi k/N} \cdot s[N-1]$

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

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

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

$|X[k]|^2 = s[N]^2 + s[N-1]^2 - 2\cos\!\left(\tfrac{2\pi k}{N}\right) s[N] \cdot s[N-1]$

位相が要らなければ複素演算ゼロで完結する。

計算量の比較

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

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

$NK < N \log_2 N \implies K < \log_2 N$

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

実装してみると

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

import math

def goertzel(samples, k, N):
    coeff = 2 * math.cos(2 * math.pi * k / N)
    s1, s2 = 0.0, 0.0
    for x in samples:
        s0 = x + coeff * s1 - s2
        s2, s1 = s1, s0
    # パワーだけ欲しい場合
    power = s2**2 + s1**2 - coeff * s1 * s2
    return power

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

何が面白いか

FFT はすべての周波数を一括して求める「グローバルな」変換だ。それに対して Goertzel は「特定の周波数を狙い撃ちにする局所的なフィルター」と捉えられる。

面白いのは、Goertzel フィルターを走らせながらサンプルを受け取るとき、中間状態 $s[n]$ が更新され続けるという点だ。N サンプル集まったら出力を取り出して初期化する、というブロック処理をするのが基本だけど、スライディング窓でリアルタイムに更新する拡張もある（Sliding DFT と呼ばれる）。

もう一つ：$k$ を整数に限る必要はない。$k$ を実数にすれば、DFT グリッド外の任意の周波数を計算できる。これは「ゾーム FFT（Zoom FFT）」の一形態にも繋がる話で、スペクトルの特定帯域を高分解能で見たいときに使える。

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

— ランキン