デジタルフィルタ入門 #10｜検出のフィルタの設計：Goertzel・整合フィルタ・ヒルベルト変換で周波数・波形・包絡線を取り出す

はじめに

#2 から #9 までのフィルタは、どれも信号を「通す・消す・均す・推定する」道具でした。出てくるのは、入力と同じ形の波形か、状態の推定値です。ところが、手元で知りたいことが波形そのものではない場面があります。

• 電話のプッシュ音（DTMF）を聞いて、どのキーが押されたかを知りたい
• 超音波を送って、返ってきたエコーの到着時刻から距離を知りたい
• 振幅が変化する搬送波から、包絡線（振幅の変化の外形）や瞬時周波数（その瞬間に位相が進む速さ）を取り出したい

どれも答えは数です。キーの番号、時刻、振幅や周波数の値です。この記事で扱う検出のフィルタは、入力から「その周波数があるか」「その波形がいつ来たか」「振幅と位相はいくつか」を、判定に使える数にして返します。

次の図は、3 つの検出のフィルタが同じ「ノイズに埋もれた信号」から何を取り出すかを並べたものです。

(a) は、振幅 0.5 の 1 kHz の正弦波に、標準偏差 1 のノイズを足した 256 点の列です。FFT（高速フーリエ変換）は 128 本のビンをすべて計算しますが、知りたいのが 1 kHz の強さだけなら、Goertzel で 1 本だけを計算すれば足ります。値は FFT の同じビンと 3 \times 10^{-13} の差で一致しました。(b) は、40 kHz のバースト波（決まった周期だけ続く正弦波）を、1 サンプルあたりの SNR（信号とノイズの電力の比）−5 dB の雑音に埋めたものです。受信波形を見ても場所は分かりませんが、整合フィルタの出力は到着時刻にちょうどピークを立てました。(c) は 1 kHz の搬送波の振幅を 5 Hz でゆっくり変えた AM 波です。ヒルベルト変換で作った解析信号（元の信号と、その位相を 90° 回した信号を、実部と虚部に組んだ複素数の列）の絶対値が、外形をそのままなぞります。

この記事のゴールは、手元の検出の問題について、知っているものが周波数か波形か帯域かを見分けて方式を選び、ブロック長やテンプレートやしきい値を式と数値から決め、誤警報と見逃しの確率で設計を判定できるようになることです。Goertzel、整合フィルタ、ヒルベルト変換の順に、式を 1 段ずつ導いて手計算で確かめ、周波数と制御工学の道具で読みます。そのうえで、DTMF の判別、超音波の到着時刻、包絡線と FM の復調を設計します。実機は使わず、すべて計算と数値シミュレーションで示します。

この記事で分かること： Goertzel がなぜ実数の 2 次 IIR 1 本で DFT の 1 ビンを計算できるのか、極が単位円の上にあるのになぜ使えるのか、ビンの周波数応答と漏れ、整合フィルタがなぜテンプレートを時間反転したインパルス応答になるのか（コーシー・シュワルツの不等式）、そのとき出力 SNR がなぜ波形の形によらずエネルギーだけで決まるのか、誤警報と検出の確率の関係（ROC）、チャープで距離の分解能が上がる理由、ヒルベルト変換の理想特性と FIR で近似したときの誤差、瞬時周波数が意味を持つ条件が分かります。

読み終えたときにできること： 周波数・波形・帯域のどれを知っているかから方式を選び、Goertzel のブロック長 N、整合フィルタのテンプレートとしきい値、ヒルベルト FIR の長さを決め、Python と C で動かせます。FFT との一致、出力 SNR、ROC、包絡線の誤差で設計を判定できます。

🔎 検出のフィルタとは

検出という言葉は、日常ではただ「見つける」ことですが、信号処理ではもう少し狭い意味で使います。定義を 1 文で書きます。

これまでのフィルタは「どの周波数を通すか」を設計しました。検出のフィルタは「何を探すか」を設計します。探すものを周波数で与えるのが Goertzel、波形で与えるのが整合フィルタ、狭い帯域にあることだけを与えて振幅と位相を取り出すのがヒルベルト変換です。

語をほどきます。

• DFT（離散フーリエ変換）：N 個のサンプルを、N 本の周波数の成分（ビン）の和に分解する計算。#1 の PSD を計算するときの土台。FFT はこれを速く計算する方法
• ビン：DFT が値を返す周波数の 1 本 1 本。k 番目のビンの周波数は k f_s / N
• 仮説：「信号がある」「ない」のように、データを見る前に立てておく 2 つ以上の説明。検出はどちらが正しいかを選ぶ問題
• 判定量（検定統計量）：仮説を選ぶために、データから計算する 1 つの数
• 誤警報：信号がないのに「ある」と判定すること。その確率を P_{FA} と書く
• 見逃し：信号があるのに「ない」と判定すること。あると判定できる確率を検出確率 P_D と書く

3 つのフィルタを、それぞれ 1 文で定義します。

• Goertzel アルゴリズム：DFT のうち 1 本のビンだけを、実数の係数を持つ 2 次の IIR を N サンプル回して計算する方法
• 整合フィルタ（matched filter）：探したい波形 s[n] を時間反転したものをインパルス応答に持つ FIR。白色ノイズのもとで、出力の SNR を最大にする
• ヒルベルト変換：すべての周波数の位相を −90° 回し、振幅は変えない線形フィルタ。実数の信号から解析信号（複素数の信号）を作り、包絡線と瞬時位相を取り出す

何を渡すと何が返るのか

これまでのフィルタとの差を並べます。

出力が波形でなく判定量になること、周波数や波形という「仮説」を前もって与えることが、これまでとの違いです。ヒルベルト変換は判定そのものはしませんが、包絡線と位相という、検出に使う量を作る道具です。

名前の由来

Goertzel アルゴリズムは、Gerald Goertzel が 1958 年に American Mathematical Monthly 誌に載せた 2 ページの短い論文「An Algorithm for the Evaluation of Finite Trigonometric Series」に由来します。題名どおり、三角関数の有限和（DFT の 1 項）を少ない計算で求める方法でした。電話のプッシュ音の判別に広く使われ、Texas Instruments のアプリケーションノートでも DTMF の検出の方法として説明されています。

整合フィルタの考え方は、レーダーの研究から生まれました。D. O. North は 1943 年に RCA 研究所の報告として、パルスを雑音の中から見分ける要因を解析し、検波（振幅を取り出す処理）の前に置く最適な選択度（どの周波数をどれだけ通すかというフィルタの特性）を定式化しました。この報告は 1963 年に Proceedings of the IEEE に再録されています。Van Vleck と Middleton は 1946 年の論文で、最適な中間周波数（受信機で搬送波を一度変換した周波数）のフィルタが「パルスのフーリエ変換の共役」になることを、目視、耳、メータのどの受信方法でも成り立つと示しました。Turin の 1960 年の解説「An introduction to matched filters」は、整合フィルタの定義と性質と合成法をまとめたものです。

ヒルベルト変換は数学者 Hilbert の名にちなみます。信号処理で使う解析信号の考え方は、Gabor の 1946 年の論文「Theory of communication」にさかのぼります。瞬時周波数をどう定義し、どこまで直感と合うかは、Boashash の 1992 年の解説論文が歴史とともに整理しています。

前提

この記事の方法がうまく働くのは、次の条件がそろうときです。

1. 探すものが前もって分かる：Goertzel なら周波数、整合フィルタなら波形、ヒルベルトなら搬送波のおおよその帯域
2. ノイズが白色に近い：白色ノイズは、どの周波数にも同じ強さで広がり、サンプルごとに独立なノイズ。整合フィルタの最適性はこれが前提
3. 観測の区間が有限：Goertzel は N サンプル、整合フィルタは M サンプルの区間で判定する
4. ヒルベルトは狭帯域の信号：帯域が搬送波の周波数より十分狭い 1 成分の信号でないと、瞬時周波数が意味を持たない

1 から 4 が破れたときに何が起きるかは、それぞれの節と設計例で「破る例」として見ます。

🎹 Goertzel アルゴリズム

1 本の周波数の強さだけが知りたいとき、FFT ですべてのビンを計算するのは遠回りです。この節では、DFT の 1 項を書き出し、それが 2 次の IIR の再帰で計算できることを導き、終端処理、極の位置、周波数応答、計算量の順に見ます。

DFT の 1 ビンを書き出す

N 個のサンプル x[0], …, x[N−1] の DFT の k 番目のビンは、次の和です。

この和は、x[n] に「周波数 \theta_k で回る単位円上の点」の共役を掛けて足したものです。x[n] にその周波数の成分が含まれていれば、回転がそろって和が大きくなり、含まれていなければ打ち消し合います。そのまま計算すると、1 ビンあたり N 回の複素数の掛け算と、N 個の e^{-j\theta_k n} の値が要ります。

和を再帰に書き直す

和を、1 サンプルずつ回る再帰の形に直します。まず、次の複素数の数列 y[n] を考えます。

初期値は y[−1] = 0 です。これは極が e^{j\theta_k} の 1 次の IIR です。1 サンプル進むごとに、それまでの値を角度 \theta_k だけ回して、新しい入力を足します。展開すると次の形になります。

指数を e^{j\theta_k(N-1)} \cdot e^{-j\theta_k n} と分けると、和の部分はちょうど X[k] です。k が整数なら e^{j\theta_k N} = e^{j2\pi k} = 1 なので、 e^{j\theta_k(N-1)} = e^{-j\theta_k} です。よって、

で、N サンプル回した y の最後の値から X[k] が得られます。ただし、これでは毎ステップ複素数の掛け算が要ります。

そこで、伝達関数で分母を実数にします。y の伝達関数の分母と分子に 1 - e^{-j\theta_k} z^{-1} を掛けると、分母の共役複素数の積が実数になります。

右辺の分母は実数の係数だけでできた 2 次の IIR、分子は 2 タップの FIR です。分母の部分を v[n] と書けば、Goertzel の再帰が出ます。

1 サンプルあたりの計算は、実数の掛け算 1 回と足し算 2 回です。入力が実数なら状態も実数のままで、三角関数の表もいりません。

終端処理

分子の FIR 1 - e^{-j\theta_k} z^{-1} は、毎ステップ計算する必要がありません。欲しいのは最後の y[N−1] だけだからです。

これを X[k] = e^{j\theta_k} y[N-1] に入れると、

です。電力（絶対値の 2 乗）だけが欲しいなら、複素数を作らずに実数だけで書けます。 |a e^{j\theta} - b|^2 = a^2 + b^2 - 2ab\cos\theta を使います。

これが終端処理です。N サンプルの間は 2 次の IIR を回し、最後に 3 回の掛け算で電力を出し、状態を 0 に戻して次のブロックへ進みます。位相が要らない検出（あるかないか）では、この電力だけで判定できます。

手計算：N = 8 で 1 周

式が本当に DFT と一致するかを、手で 1 周回して確かめます。fs = 8 kHz、N = 8、k = 1（1 kHz のビン）です。入力は 1 kHz の正弦波を粗く 3 値で表した 8 個です。

係数は \theta_1 = 2\pi/8 = \pi/4 より、 2\cos(\pi/4) = \sqrt{2} \approx 1.41421 です。再帰を回します。

終端処理で電力を出します。

複素数の値は X[1] = e^{j\pi/4}(-6.82843) + 4.82843 = -4.82843\,j です。FFT で計算した X[1] も −4.82843j、電力は 23.314 で一致します。残りのビンの電力は、k = 7 は k = 1 の鏡像で同じ 23.314、k = 3 と 5 が 0.686、k = 0、2、4、6 が 0 で、粗い 3 値の波形に含まれる高調波の分だけが残っています。同じ計算を Python で確かめるコードは、実装の節にあります。

極は単位円の上：条件付き安定

Goertzel の 2 次の IIR の伝達関数を、#1 の極と零点で見ます。

分母を 0 にする z は z = e^{\pm j\theta_k} で、極はちょうど単位円の上にあります。#1 で見たとおり、極が単位円の内側にあれば安定、外側なら発散、単位円の上は安定の境目です。この IIR を止めずに回し続けると、ビンの周波数の入力に対して状態が限りなく大きくなります。

それでも Goertzel が使えるのは、N サンプルで止めて読み、状態を 0 に戻すからです。この性質を条件付き安定（決まった有限の時間だけ回すなら値が有限にとどまる）と呼びます。次の図で、状態の伸び方を見てください。

図の (c) で、ビンの周波数の入力（振幅 1）は |y[n]| が 1 サンプルに 1/2 ずつ直線的に増え、n = N − 1 = 63 で 32 = N/2 になりました。半ビンずれた入力は、回転がそろわないので約 20 で頭打ちになり、そのあと減っていきます。3 ビン離れた入力はほとんど積み上がりません。つまり Goertzel は、その周波数の成分だけを N サンプル積分する「有限時間の積分器」で、積分の時間 N が周波数の分解能を決めます。

ビンの周波数応答と漏れ

ビン k に周波数 f の正弦波が入ったときの |X[k]| を、f の関数として描くと、ビンの周波数応答になります。矩形の窓（N サンプルを同じ重みで足す）の DFT なので、形はディリクレ核（sinc 関数を周期的にした形）です。

負の周波数の成分の重なりは、ビンが 0 Hz とナイキスト周波数（fs/2）から十分離れていれば無視できるので、式では省いています。次の図の (a) は、fs = 8 kHz、N = 205 のビン 18（702.4 Hz）の応答です。点線は DTMF の 697、770、852 Hz です。

(a) で読み取れることを並べます。

• 主ローブの幅：ビン間隔 f_s/N = 39.0 Hz の 2 倍。隣のビンの中心で応答がちょうど 0 になる
• 第 1 サイドローブ：サイドローブは主ローブの外側に並ぶ小さな山で、その最初の山が矩形窓で −13.3 dB。隣の DTMF のトーン 770 Hz は、ビン 18 に −17.2 dB で漏れ込む
• ハン窓：入力にハン窓（両端を 0 に絞る重み）を掛けてから Goertzel を回すと、第 1 サイドローブは −31.5 dB に下がり、770 Hz の漏れは −23.6 dB になる。代わりに主ローブの幅が 2 倍になる

(b) は「破る例」です。N = 64 で、k = 10 の整数ビンに乗る正弦波は 1 本のビンにだけ立ちます。k = 10.5 のように半ビンずれると、最寄りの 2 本のビンがどちらも −4.0 dB に下がり（理論値は 20\log_{10}(2/\pi) = -3.9 dB、スカロッピング損失と呼ぶ）、遠く離れたビンにも −27〜−30 dB の漏れが広がります。ハン窓を掛けると、最寄りのビンの損失は −1.5 dB（理論値 −1.4 dB）に減り、遠いビンへの漏れは −70 dB まで下がります。周波数を狙い撃ちしたいなら、一般化 Goertzel で k = 10.5 をそのまま指定すれば、損失は 0 dB です。

窓と漏れの関係は、#1 の Welch 法で PSD を計算するときに窓を掛けたのと同じ話です。Goertzel は DFT の 1 項なので、DFT の性質がそのまま当てはまります。

計算量：何本までなら FFT より軽いか

N_b 本のビンが要るとき、Goertzel を N_b 本並べるのと、FFT を 1 回回すのとでは、どちらが軽いかを数えます。掛け算は実数の回数で、複素数の掛け算は実数 4 回と数えます。

この数え方では、 N_b \lt 2\log_2 N 、N = 256 なら 16 本より少ないとき Goertzel が軽くなります。実数入力の FFT は工夫で掛け算を約半分にできるので、境目は 8 本前後まで下がります。数のほかにも違いがあります。Goertzel は入力が来るたびに 1 サンプル分を計算できるので、N 点を貯めるバッファがいりません。N を 2 のべき乗に限る必要もなく、ビンの周波数を自由に選べます。

🎯 整合フィルタ

Goertzel は「周波数」を知っているときの検出でした。整合フィルタは「波形」を知っているときの検出です。この節では、検出を仮説の選択として書き、出力 SNR を最大にするフィルタをコーシー・シュワルツの不等式から導き、誤警報と検出の確率、周波数応答、有色ノイズ、パルス圧縮の順に見ます。

検出を仮説の選択として書く

受信した M 個のサンプル x[0], …, x[M−1] について、2 つの仮説を立てます。

ガウスノイズは、値の分布が正規分布（釣り鐘形）に従うノイズです。この 2 つのどちらが正しいかを、x から計算した 1 つの数で決めたい、というのが検出の問題です。

最初の候補は、x とテンプレートの内積（対応するサンプルを掛けて足したもの）です。

x にテンプレートと同じ形が含まれていれば、掛け算の符号がそろって Λ は大きくなり、ノイズだけなら正負が打ち消し合います。既知の波形と白色ガウスノイズのときは、この内積をしきい値と比べる方法が、誤警報の確率を決めたうえで検出の確率を最大にする意味で最適であることが知られています（Kay の『Detection Theory』の教科書が導出を与えています）。ここでは、もう少し直接的に「出力 SNR を最大にする FIR」として同じ答えを導きます。

コーシー・シュワルツで出力 SNR を最大にする

長さ M の FIR h[m] を x に掛けたとき、テンプレートの最後のサンプルが入った時刻の出力は、

です。信号の分とノイズの分に分けます。逆順のテンプレートを g[m] = s[M-1-m] と置くと、

出力 SNR を、信号の分の 2 乗とノイズの分散の比で定義します。

ここでコーシー・シュワルツの不等式を使います。

2 つのベクトルの内積は、長さの積に向きのそろい具合（角度の余弦）を掛けたもので、余弦は 1 を超えない、という幾何の事実を式にしたものです。a = h、b = g とすると、分子は \sum h^2 \sum g^2 以下なので、

です。等号は h が g の定数倍のとき、つまり次の形のときだけ成り立ちます。

これが整合フィルタです。インパルス応答は、テンプレートを時間反転したものです。定数 c は SNR に影響しないので、c = 1 とします。

結論を 2 つ書き出します。

• 出力 SNR の上限は E_s/\sigma^2 ：テンプレートのエネルギーとノイズの分散の比だけで決まり、波形の形には依らない。矩形のパルスでも正弦波のバーストでも、エネルギーが同じなら同じ SNR になる
• 1 サンプルあたりの SNR が低くても、長さで稼げる：振幅 A の正弦波なら E_s = MA^2/2 なので、M を 2 倍にすると出力 SNR は 3 dB 上がる

相関器と同じもの

整合フィルタの出力を、時刻 n の式で書き直します。

右の形は、x の n−M+1 から始まる M 個とテンプレートとの内積です。つまり、テンプレートを x の上で 1 サンプルずつずらしながら内積を取る相互相関（2 つの列の、ずれごとの内積を並べたもの）と同じ計算で、ピークは「テンプレートの始まり＋M−1」の時刻に立ちます。整合フィルタは「FIR として畳み込む」見方、相関器は「ずらして内積を取る」見方で、出力の列は同じです。SciPy の scipy.signal.correlate との一致は、実装の節で確かめます。

次の図の (a) は、長さ 64 の正弦バーストを、1 サンプルあたりの SNR −4.1 dB（ E_s/\sigma^2 = 14 dB）でノイズに埋めたものです。

(b) は、エネルギーをそろえた 3 つの形（矩形パルス、正弦バースト、13 ビットの Barker 符号（自分とのずれた内積が小さくなるように選んだ ±1 の列）を 5 サンプルずつ伸ばしたもの）で、出力 SNR を 2 万回の計算から求めたものです。すべての点が E_s/\sigma^2 の直線に ±0.1 dB で乗りました。(c) は、正弦バーストに別の h を当てたときの損失です。周波数が 2 % ずれると −1.5 dB、5 % ずれると −12.8 dB です。8 周期のバーストでは、5 % のずれが最後までに 0.4 周期の位相のずれになり、前半と後半で掛け算の符号が逆になるからです。位相が 45° ずれると \cos^2 45^\circ で −3 dB、テンプレートの前半しか使わないとエネルギーが半分で −3 dB です。位相が 90° ずれると、内積はちょうど 0 になります。

手計算：長さ 4 のテンプレート

小さな例で、ピークの位置を手で求めます。テンプレートは長さ 4 の Barker 符号です。

受信は 10 個で、n = 3 から s にノイズが少し乗ったものが入っています。

始まりの位置 n₀ をずらしながら、内積 c[n_0] = \sum_{l} s[l]\, x[n_0+l] を計算します。

最大は n₀ = 3 の 4.0 で、 E_s = 4 とほぼ同じ値です。整合フィルタとして因果に回すと、このピークは n = n₀ + M − 1 = 6 の出力に出ます。Python の matched_filter の出力は (0, 0.2, −0.3, 1.4, −0.2, −0.7, 4.0, −1.2, 0, 1.5) で、n = 3 以降が上の表の c[n−3] と一致します。

しきい値と ROC

出力のピーク Λ をしきい値 η と比べて「ある」「ない」を決めます。白色ガウスノイズのもとでは、Λ は正規分布に従います。

平均の差は E_s 、標準偏差は \sigma\sqrt{E_s} なので、誤警報と検出の確率は次のとおりです。

η を動かすと P_{FA} と P_D が一緒に動きます。その組を描いた曲線を ROC（受信者動作特性）と呼びます。η を下げれば見逃しは減りますが誤警報が増え、その交換の割合を d が決めます。

(a) の丸は、長さ 32 のテンプレートで 40 万回ずつ判定した結果です。 P_{FA} が 1e−4 から 1e−1 の 4 点すべてで、 P_D の計算値は理論と 0.001 以内で一致しました。(b) から、 P_D = 0.9 に要る E_s/\sigma^2 は、 P_{FA} = 10^{-2} で 11.1 dB、 10^{-4} で 14.0 dB、 10^{-6} で 15.6 dB です。

しきい値は、許せる誤警報の確率から決めるのが手順です。1 秒に何回判定するかを掛けると、1 時間あたりの誤警報の回数になります。たとえば 1 秒に 1000 回判定して、誤警報を 1 時間に 1 回までにしたいなら、 P_{FA} \approx 2.8 \times 10^{-7} です。σ は前もって分からないことが多いので、信号のない区間でノイズの分散を推定してから η を決めます。

周波数応答と、有色ノイズで破れる例

整合フィルタの周波数応答は、テンプレートのスペクトル S(e^{j\theta}) で書けます。s が実数なら、時間反転は周波数では共役になります。

振幅は |H| = |S| で、テンプレートが強く持っている周波数ほどよく通します。つまり整合フィルタは「テンプレートの帯域を通すバンドパス」でもあります。#4 のバンドパスとの違いは、振幅の形までテンプレートに合わせていることと、位相を打ち消していることです。

(a) は、40 kHz の 8 周期のバースト（fs = 400 kHz）について、整合フィルタの振幅と、出力 SNR が最大になるように Q 値を選んだ #4 のバンドパスを重ねたものです。バンドパスの出力のピークで判定したときの SNR は、最良の Q 値 19.75 でも整合フィルタより 1.5 dB 低くなりました。バンドパスは主ローブの形しか似せられず、位相もそろえないからです。

(b) は、前提 2（白色ノイズ）を破る例です。ノイズを AR(1)（1 つ前の値の 0.95 倍に新しい乱数を足す、低い周波数が強いノイズ）にしました。白色を仮定した整合フィルタの出力は、ノイズの低い周波数の揺れがそのまま出て、正しい時刻のピークが埋もれます。理論の出力 SNR は 1.4 dB で、400 回のうちピークが正しい時刻に立ったのは 7 % でした。

直し方は、#8 の前白色化と同じ発想です。ノイズのスペクトルを平らにする固定のフィルタ 1 - 0.95z^{-1} を、受信とテンプレートの両方に掛けてから整合を取ります。白色化した後のノイズは白色なので、整合フィルタの最適性が戻ります。理論の出力 SNR は 15.2 dB（有色ノイズのもとでの上限）まで上がり、正しい時刻のピークは 98 % になりました。ノイズのスペクトルの形は、信号のない区間を #1 の Welch 法で測って決めます。

チャープでパルスを圧縮する

出力 SNR はエネルギーで決まりますが、ピークの鋭さ、つまり時刻の分解能はテンプレートの帯域で決まります。同じ長さとエネルギーで帯域を広げる波形がチャープ（周波数が時間とともに変わる波形）です。整合フィルタに通すと、長いチャープが短いピークに縮むので、これをパルス圧縮と呼びます。レーダーや超音波で、送信の電力を上げずに距離の分解能を上げる原理です。

(b) を #1 の群遅延で読みます。チャープは低い周波数が先、高い周波数が後に来るので、その群遅延（各周波数の成分がどれだけ遅れて来るか）は周波数とともに増えます。整合フィルタの位相は -\angle S - \theta(M-1) なので、群遅延は逆に周波数とともに減ります。2 つを足すと、どの周波数でも M - 1 = 79 サンプル（197.5 μs）で一定になりました。

この式が同じことを言っています。整合フィルタを通した後のスペクトルは、位相が直線（群遅延が一定）になり、すべての成分が同じ時刻 M−1 に頂上をそろえます。整合フィルタは、群遅延を一定にして「ピークを立てる」操作です。

(c) は、100 μs 離れた同じ強さの 2 つのエコーです。単独のエコーの整合出力の包絡線の幅（頂上の半分以上の区間）は、バーストで 202.5 μs、チャープで 62.5 μs でした。バーストの山は 2 つのエコーをつないで 1 つにし、チャープは 2 つに分けます。チャープの幅は帯域 B = 20 kHz の逆数 50 μs の程度で、長さ 200 μs には依りません。音速 343 m/s の往復なら、62.5 μs は約 1.1 cm です。

🌀 ヒルベルト変換と解析信号

整合フィルタの callout で、位相が分からないときは sin と cos の 2 つの出力の 2 乗和を取ると書きました。どんな実数の信号にも、その「相方の sin」を作るのがヒルベルト変換です。この節では、理想特性からインパルス応答を導き、FIR で近似し、解析信号から包絡線と瞬時周波数を取り出し、FFT での実装と比べます。

理想特性：すべての周波数を −90° 回す

ヒルベルト変換の周波数応答は、次の形で定義します。

振幅は全周波数で 1、位相は正の周波数で −90°、負の周波数で +90° です。cos を入れると sin が、sin を入れると −cos が出ます。どの周波数の正弦波も、振幅を変えずに 4 分の 1 周期だけ遅らせる操作です。

インパルス応答は無限長で、因果でない

周波数応答を逆変換して、インパルス応答を求めます。

正の側と負の側に分けて積分します。n ≠ 0 なら \int_0^{\pi} e^{j\theta n} d\theta = \bigl(e^{j\pi n} - 1\bigr)/(jn) で、 e^{j\pi n} = (-1)^n です。−j を掛けた正の側と、+j を掛けた負の側は同じ値 (1 - (-1)^n)/n になり、足して 2π で割ると次の形になります。

n = 0 では、被積分関数が奇関数なので 0 です。この応答には 2 つの困った性質があります。

• 無限に続く： 1/n でしか減らないので、どこで切っても切り捨てた分が残る
• 因果でない：n が負の側にも同じ大きさの値があり、未来の入力を使う

窓で切って遅らせた FIR

実際に使うときは、−D から D までの 2D+1 個で切り、端の切り捨ての影響を窓で和らげ、D サンプル遅らせて因果にします。

m は 0 から L−1 までです。

係数は中央に対して奇対称（ h_L[D+m] = -h_L[D-m] ）で、中央から偶数だけ離れた係数は 0 です。この形の FIR は、遅れ D を除いた位相がどの周波数でもちょうど −90° になります。誤差はすべて振幅に出て、θ = 0 とナイキスト周波数では振幅が必ず 0 に落ちます。どの帯域で振幅が 1 に近いかは、次の「周波数と制御工学の道具で見る」の節で、L を変えたボード線図で見ます。

掛け算の数は、対称性と 0 の係数を使うと (D+1)/2 回まで減ります。L = 63 なら 16 回です。

解析信号、包絡線、瞬時周波数

元の信号とヒルベルト変換を、実部と虚部に組みにした複素数の列を解析信号と呼びます。

周波数で見ると、 Z = X(1 + \mathrm{sgn}\,\theta) なので、正の周波数は 2 倍、負の周波数は 0 になります。実数の信号のスペクトルは正と負で対称なので、負の側を捨てても情報は失われません。このおかげで、z[n] の絶対値と偏角から、振幅と位相を 1 サンプルごとに読めます。

瞬時周波数は位相の差 \psi[n] - \psi[n-1] ですが、ψ は ±π で飛ぶので、そのまま引くと #6 の角度の巻き戻しと同じ問題が起きます。 z[n]\overline{z[n-1]} の偏角を取れば、差が最初から −π から π に収まるので、巻き戻しの処理がいりません。差分はノイズを持ち上げるので、ノイズが多いときは #9 の SG 微分で位相を微分する手もあります。

手計算：4 点の DFT で解析信号を作る

FFT で解析信号を作る手順を、4 点で手計算します。手順は「DFT → 負の周波数のビンを 0、正のビンを 2 倍（直流とナイキストはそのまま）→ 逆 DFT」です。N = 4 では、ビン 0 が直流、ビン 1 が正、ビン 2 がナイキスト、ビン 3 が負の周波数です。

まず fs/4 の cos、 x = (1, 0, -1, 0) です。DFT は X = (0, 2, 0, 2) で、ビン 1 を 2 倍、ビン 3 を 0 にすると Z = (0, 4, 0, 0) です。逆 DFT は z[n] = \tfrac{1}{4} \cdot 4\, e^{j\pi n/2} なので、

です。虚部 (0, 1, 0, −1) は sin で、包絡線は振幅 1 で一定です。

次に、直流 1 を足した x = (3, 1, -1, 1) = 1 + 2\cos(\pi n/2) です。DFT は X = (4, 4, 0, 4) 、手順を当てると Z = (4, 8, 0, 0) 、逆 DFT は z[n] = 1 + 2e^{j\pi n/2} です。

振幅 2 の cos なのに、包絡線は 1 から 3 まで揺れました。直流は周波数 0 で正でも負でもないので、ヒルベルト変換では相方を作れず、実部にだけ残るからです。包絡線を取る前に、#3 のハイパスで直流を落としておく必要があります。どちらの値も、Python の analytic_fft と scipy.signal.hilbert の出力と一致しました。

FFT での実装と FIR での実装

SciPy の scipy.signal.hilbert は、上の手計算と同じ「FFT で負の周波数を 0、正を 2 倍」で解析信号を返します。ブロック全体を一度に処理するので、ブロックの中では FIR の窓による誤差がありません。代わりに、DFT はブロックを周期的につないだものとして扱う（円状畳み込み）ので、ブロックの端と端が不連続だと、そのつなぎ目で包絡線が乱れます。

(a) の 1030 Hz は 256 点に整数周期が収まらないので、FFT の包絡線は最初の 16 サンプルで最大 4.3 % ずれました。中央では 0.13 % です。FIR（L = 63）は、最初の 2D = 62 サンプル（窓が入力で満たされるまで）が過渡で、それ以降は 0.15 % 以内でした。(b) の 1000 Hz は 256 点にちょうど 32 周期収まるので、FFT の包絡線は端まで 10^{-14} の誤差です。

2 つの実装の使い分けは、ブロックで後から処理するなら FFT（端の数 % を捨てるか、前後に余白を取る）、1 サンプルずつ流れてくるなら FIR（遅れ D と最初の 2D サンプルの過渡を受け入れる）です。

瞬時周波数が意味を持つ条件

瞬時周波数は、1 つの成分だけでできた狭帯域の信号でないと、直感どおりの値になりません。Boashash は、1 成分の信号（monocomponent）と多成分の信号（multicomponent）を区別し、瞬時周波数の解釈がこの区別に依ることを論じています。次の図は、満たす例と破る例と、包絡線の誤差の条件です。

(a) の FM（周波数が 1 kHz を中心に ±100 Hz、5 Hz で揺れる）では、計算した瞬時周波数は真の値と 0.2 Hz 以内で重なりました。

(b) は前提 4 を破る例で、1000 Hz と 1150 Hz の 2 つの正弦波の和です。2 つ目の振幅を a とすると、解析信号は e^{j\theta_1 n}(1 + a\,e^{j\Delta\theta\, n}) で、その偏角の進みは次の値の間を往復します。

Δf = 150 Hz です。a = 0.5 では 850〜1050 Hz で、計算値も同じでした。a = 0.9 では理論の最小が −350 Hz で、2 つの成分が打ち消し合う瞬間に瞬時周波数が負になります。サンプルした点の最小は −196 Hz でした（打ち消しの谷が 1 サンプルより鋭いため）。どちらの成分の周波数とも違う値で、「その瞬間の周波数」という直感には合いません。

(c) は包絡線の条件です。1 kHz の搬送波を変調周波数 f_m で AM 変調すると、 1000 \pm f_m Hz に側帯波（変調で搬送波の両側に立つ成分）が立ちます。下側帯波が正の周波数にあるうち（ f_m \lt 1000 Hz）は、FFT の包絡線は 10^{-12} の誤差です。 f_m が 1000 Hz を超えて下側帯波が負の周波数に回ると、ヒルベルト変換はそれを捨ててしまい、誤差は 0.5 に跳ねます。FIR（L = 63）は振幅が 1 に近い帯域が約 207 Hz より上なので、下側帯波がそれを下回る f_m = 800 Hz から誤差が増え、900 Hz で 0.07 になりました。包絡線が正しく取れる条件は、信号の帯域全体が 0 Hz より上、FIR ならその通過域の中にあることです。

📊 周波数と制御工学の道具で見る

3 つの検出のフィルタを、これまでの回の道具で読み直します。極の位置、群遅延、ボード線図、そして制御ループの中での使われ方です。

Goertzel は Q 値が無限大のバンドパスを N で止めたもの

Goertzel の極は単位円の上にありました。#4 の 2 次のバンドパスの極は、単位円の少し内側にあり、内側に入るほど帯域が広く Q 値が低くなります。極を単位円の上まで押し出したバンドパスは Q 値が無限大で、止めなければ発散します。Goertzel は、それを N サンプルだけ回して読むことで、帯域を「N で決まる幅」に区切ったものと読めます。

ビンの −3 dB の全幅は約 0.886\, f_s/N です。N = 205 のビン 18（702.4 Hz）では 34.6 Hz で、同じ幅のバンドパスなら Q 値は約 20 です。Q 値 20 のバンドパスの振幅は、時定数（Q 値を \pi f_0 で割った値）9.1 ms で落ち着くので、3 倍の 27 ms でほぼ定常になります。Goertzel の 1 ブロック 25.6 ms とほぼ同じ時間です。違いは、バンドパスは 1 サンプルごとに出力を出し続けるのに対し、Goertzel は N ごとに区切って 1 個の数を出すこと、そしてビンの応答の形（sinc 形のサイドローブ）です。隣のトーンの漏れを窓で調整できるのは Goertzel の側です。

整合フィルタは群遅延を一定にしてピークを立てる

パルス圧縮の節で見たとおり、整合フィルタを通した後のスペクトルは |S|^2 e^{-j\theta(M-1)} で、位相が直線、群遅延が M−1 で一定です。#1 で、群遅延が周波数で違うと波形が崩れると見ました。整合フィルタはその逆で、テンプレートの群遅延の凸凹を打ち消してすべての成分を同じ時刻に集め、波形をわざと 1 点に潰します。振幅も |S|^2 と強い成分ほど重く足すので、ノイズに対して得になる成分から順に集めていることになります。

ヒルベルト FIR のボード線図

ヒルベルト FIR の近似の質を、ボード線図（振幅と位相を周波数に対して描いた図）で見ます。fs = 8 kHz、ハミング窓で、長さ L = 31、63、127 を比べました。位相は遅れ D の分を除いて描いています。

振幅が 1 から 1 % 以内に収まる帯域は、L = 31 で 432〜3568 Hz、L = 63 で 207〜3793 Hz、L = 127 で 102〜3898 Hz でした。L を 2 倍にすると、使えない低域の幅がほぼ半分になります。位相は、奇対称の FIR なので 3 つとも −90° からのずれが 10^{-12} 度以下でした。設計の手順は、使いたい最も低い周波数から L を決めることです。fs = 8 kHz で 250 Hz から上を使いたいなら、1 % 帯域の下端が 207 Hz の L = 63（遅れ 31 サンプル、3.9 ms）で足ります。

振幅を全帯域で 1 にして位相の側に誤差を逃がす作り方もあり、オールパス（振幅が 1 で位相だけを変える IIR）を 2 本並べて、位相の差が 90° になるように設計します。これは #11 のオールパスの節で扱います。

制御ループの中の検出器：位相検出

検出のフィルタは、測定の最後に置くだけでなく、制御ループの中にも入ります。代表が PLL（位相同期ループ。発振器の位相を入力の位相に追従させるループ）の位相検出器です。入力 x の解析信号 z を作り、ループ側の発振器の位相 \psi_{\mathrm{ref}} との差を \arg\bigl(z\, e^{-j\psi_{\mathrm{ref}}}\bigr) で直接読めば、それが位相誤差です。掛け算型の位相検出器と違って 2 倍の周波数の成分が出ないので、後段のループフィルタ（位相誤差から発振器の制御量を作るローパス）で落とす必要がありません。ただし、ヒルベルト FIR の遅れ D はループの中の遅れになるので、#2 で見たとおり、ゲイン交差周波数 f_{gc} （開ループのゲインがちょうど 1 になる周波数）で 360^\circ f_{gc} D T だけ位相余裕（位相が −180° までどれだけ残っているか）を減らします。

🛠️ 設計例

ここまでの式を使って、3 つの検出を要件から設計します。どれも「要件 → 方式 → パラメータ → 検証」の順に進め、破る例を 1 つずつ添えます。

設計例 1：DTMF（プッシュ音）の判別

要件。 電話のプッシュ音は、低い群（行）の 4 つの周波数から 1 つ、高い群（列）の 4 つから 1 つを同時に鳴らす 2 音で、16 個のキーを表します。

fs = 8 kHz で、どのキーかを判定します。周波数の許容差は、Texas Instruments のアプリケーションノート（SPRA096A）が Bellcore の仕様として挙げている値を使います。公称から ±1.5 % 以内のトーンは必ず受け入れ、±3.5 % 以上ずれたトーンは拒否します。SNR はノイズに対する 2 音の合計の電力で数えます。

方式。 8 本の周波数だけが要るので、FFT ではなく Goertzel を 8 本並べます。

N の選び方。 最も近い 2 音（697 と 770 Hz）の間隔は 73 Hz なので、ビン間隔 f_s/N をそれより十分狭くします。そのうえで、8 つのトーンができるだけビンの中心に乗る N を探しました。N を 90 から 259 まで動かし、8 つの k = f N / f_s と最寄りの整数の差の最大を比べると、最小は N = 205 の 0.27 ビンでした。次は 113（0.29）、185（0.30）で、2 のべき乗の 256 は 0.36 です。N = 205 はビン間隔 39.0 Hz、1 ブロック 25.6 ms で、k は 18、20、22、24、31、34、38、42 です。公称の周波数での自分のビンの損失は、最大で 770 Hz の −1.0 dB でした。

判定のルール。 Goertzel の電力を (N/2)^2 で割ると、振幅 a のトーンでほぼ a^2 になります。この正規化した電力を使い、次の 3 つをすべて満たしたときだけキーを返します。

1. エネルギー：行と列の最大のビンの和が、ブロックの平均電力の 2 倍（2 音の a^2 の和に当たる）の半分以上
2. 相対ピーク：各群で、最大のビンが 2 番目より 6 dB 以上大きい
3. ツイスト：行と列の電力の比が ±8 dB 以内（ツイストは、行の音と列の音の大きさの差。TI のノートは向きで 8 dB と 4 dB を分けているが、ここでは両方 8 dB にした）

次の図の (a) は、キー 5（770＋1336 Hz、それぞれ +1 % と −1 % ずらした）を SNR 10 dB で入れたときの 8 本のビンです。

(a) では、770 Hz と 1336 Hz のビンが 0 dB 近くに立ち、残りのビンは −19 dB 以下です。相対ピークの 6 dB は十分に満たします。

破る例：固定の N と % の許容差。 (b) は、各トーンを公称からずらしたとき、自分のビンの応答がどう落ちるかを、ずれを % で横軸に取って描いたものです。灰色の 8 本が N = 205 です。ビン間隔は Hz で一定なのに、許容差は % で決まるので、高い周波数ほど同じ % のずれが多くのビンに当たります。

右の列の「／」の左が ±1.5 % の最悪、右が ±3.5 % の最良です。N = 205 では、1633 Hz の ±1.5 %（受け入れたい）が −11.9 dB まで落ちる一方、770 Hz の ±3.5 %（拒否したい）は −2.8 dB しか落ちません。770 Hz のビン 20 の中心 780.5 Hz は公称より上にあり、+3.5 % の 797 Hz は中心を挟んだ反対側の 16.5 Hz（0.42 ビン）の位置です。公称の 770 Hz も中心から 10.5 Hz（0.27 ビン）離れているので距離に大差がなく、0.42 ビンのずれでは −2.75 dB しか落ちません。1 つのしきい値では両立しません。

直し方：トーンごとに N を変える。 一般化 Goertzel は k が整数でなくてよいので、トーンごとに N を変えられます。窓の中に入る周期の数をそろえて N_i = \mathrm{round}(24 f_s / f_i) とすると、どのトーンでも「1 % のずれ＝0.24 ビン」になり、応答が % のずれだけで決まります。N は 697 Hz で 275（34.4 ms）、1633 Hz で 118 です。表の右の列のとおり、8 つすべてで ±1.5 % が −1.9 dB、±3.5 % が −15 dB にそろいました。

(c) は、±1.5 % 以内でランダムにずらしたトーンで、正しいキーを返した割合です（各点 600 回）。N = 205 は SNR を上げても 0.90 で頭打ちになり、ずれが大きい組み合わせでエネルギーの条件を満たせません。24 周期の方式は、SNR 6 dB で 0.998、9 dB 以上で 1.0 でした。±3.5 % ずらしたトーンを誤って受け入れた割合は、どちらの方式もすべての SNR で 0 でした。SNR 20 dB で ±1.5 % ちょうどずらしたときの受け入れは、N = 205 で 43 %、N = 105 で 58 %、24 周期で 100 % です。N = 105 は ±3.5 % を 28 % 受け入れてしまいました。TI の SPRA096A が、列の周波数ごとに ±9 Hz ずらした 2 本のビンを置いて主ローブを広げているのも、同じ問題への対処です。

隣のキーへの漏れ。 行の 4 つのトーンを公称の周波数で入れたとき、隣の行のビンに漏れ込む最大の値です。

どちらも 17 dB 以上の差があり、相対ピークの 6 dB に十分な余裕があります。N = 205 では隣のトーンがサイドローブの谷の近くに来るものがあり、漏れが −17〜−25 dB とばらつきます。24 周期の方式では −18〜−19 dB にそろいました。

計算量。 8 本の Goertzel は、N = 205 で 1640 回、24 周期で \sum N_i = 1504 回の掛け算（それぞれ終端処理の 24 回を足す）です。256 点の FFT は複素の掛け算 1024 回、実数で 4096 回です。トーンごとに N が違うときは、8 本の Goertzel を最新のサンプルで同時に終わるように始点をずらし、ブロックの終わりをそろえて判定します。

実用に足すもの。 TI のノートは、このほかに 2 つの確認を入れています。音声や音楽を誤ってキーと判定しないように、検出した 2 つのトーンの 2 倍の周波数（2 次高調波）のビンを計算し、それが小さいことを確かめる確認と、トーンが決まった時間以上続いたことを 2 ブロック続けて確かめる長さの確認です。この記事の数値シミュレーションには入れていません。

設計例 2：超音波の到着時刻（TOF）

要件。 40 kHz の超音波を 8 周期（200 μs）送り、返ってきたエコーの到着時刻をノイズの中から求めます。fs = 400 kHz（1 サンプル 2.5 μs）、エコーの振幅は 1、ノイズは白色ガウスです。1 サンプルあたりの SNR を (A^2/2)/\sigma^2 で数えます。テンプレートの長さは M = 80 なので、 E_s/\sigma^2 は 1 サンプルあたりの SNR より 19.0 dB 大きくなります。誤差が搬送波の半周期（±12.5 μs）を超えたものを外れと数えます。

比べる 3 つの方式。

• バンドパス＋包絡線のしきい値：#4 のバンドパス（Q 値 8）を通し、ヒルベルト FIR（L = 31）で包絡線を取り、ノイズのないときの頂上の半分を最初に超えた時刻。そのときの遅れ（18 サンプル）は前もって測って引いておく
• 整合フィルタ（バースト）：送った 8 周期をテンプレートにして、出力の最大の時刻から M−1 を引く
• 整合フィルタ（チャープ）：30→50 kHz のチャープを同じ長さとエネルギーで送り、同じように求める

(b) の外れ率を表にします（各 600 回）。

外れなかったものの RMS 誤差は、バンドパスで 5〜8 μs、整合フィルタは SNR 0 dB で 0.15〜0.18 μs、2.5 dB 以上ではサンプル単位で 0 でした。1 サンプル 2.5 μs は、音速 343 m/s の往復で約 0.43 mm の距離に当たります。

結果の読み方を 3 つ書きます。

• バンドパスのしきい値は、雑音で先に超えるか、立ち上がりの傾きの分だけずれる：包絡線がしきい値を横切る時刻は、ノイズが乗ると前後に動く。立ち上がりがゆるやかなほど、同じノイズで時刻が大きく動く
• 整合フィルタの外れは、1 周期ずれ：整合フィルタのノイズのない出力は、テンプレートの自己相関（自分自身とのずれごとの内積）になる。8 周期のバーストの自己相関は、1 周期ずれた位置に頂上の 7/8 の山がある。SNR が低いと、ノイズでその隣の山が勝ち、25 μs（約 4.3 mm）ずれる
• チャープは隣の山が低い：帯域が広いぶん自己相関の頂上が鋭く、隣の山が低いので、1 周期ずれが起きにくい。外れ率が 5 % になる SNR は、バーストの約 −1.3 dB に対してチャープは約 −4.9 dB で、3.6 dB 低かった

設計例 3：AM の包絡線と FM の復調

要件。 1 kHz の搬送波（fs = 8 kHz）について、(1) 振幅が 20 Hz で揺れる AM（変調度 0.5、振幅は 0.5〜1.5）の包絡線を取る、(2) 周波数が 10 Hz で ±200 Hz 揺れる FM の周波数を取り出す。ノイズは標準偏差 0.02 です。

包絡線の 2 つの方式。 1 つは、全波整流（絶対値を取る）して #2 の 2 次バターワース 100 Hz のローパスを通し、整流の平均の比 \pi/2 を掛ける方法です。正弦波の絶対値の平均は振幅の 2/\pi 倍なので、これで振幅に戻ります。もう 1 つは、ヒルベルト FIR（L = 63）で解析信号を作り、その絶対値を取る方法です。

遅れをそろえて比べた包絡線の RMS 誤差は、整流＋ローパスが 0.049（遅れ 18 サンプル、2.25 ms）、ヒルベルト FIR が 0.020（遅れ 31 サンプル、3.9 ms）でした。整流＋ローパスの誤差のほとんどは、振幅の目盛りのずれです。1 周期が 8 サンプルしかないと、サンプル点での |cos| の平均は 2/\pi = 0.637 ではなく 0.604 になり、包絡線が約 5 % 低く出ます（図 (a) のオレンジの頂上が 1.43 にとどまるのはこのため）。この値は搬送波とサンプルの位相の関係でも変わります。ノイズのない包絡線そのものを同じローパスに通した場合との差は RMS 0.0013 で、ローパスの遅れによる形のゆがみはわずかでした。整流で 2 kHz 以上に生じる成分も、100 Hz のローパスで −56 dB 以下に落ちています。fs が搬送波の数十倍あれば、目盛りのずれは小さくなります。

包絡線だけが目的で、変調の周波数が搬送波より十分低く、fs が搬送波より十分高いなら、整流＋ローパスで足りることが多く、計算も軽く済みます。ヒルベルトが要るのは、包絡線と一緒に位相が要るとき、変調が速くて整流＋ローパスのリプル（ローパスで落としきれない搬送波の 2 倍の成分）と遅れの釣り合いが取れないとき、そして次の FM のように瞬時周波数が要るときです。

FM の復調。 (b) では、ヒルベルト FIR の解析信号から \arg(z[n]\overline{z[n-1]}) で瞬時周波数を出しました。そのままの値は真の周波数との RMS 誤差が 26 Hz で、差分がノイズを持ち上げています。16 点の移動平均（#2）を通すと 2.3 Hz になりました。

破る例。 狭帯域の条件を破る 2 成分の信号では、瞬時周波数が負にもなることを、ヒルベルト変換の節の図 (b) で見ました。包絡線も、下側帯波が FIR の通過域（L = 63 で約 207 Hz 以上）を下回ると崩れます（同じ図の (c)）。搬送波が低いときや変調が速いときは、L を伸ばすか、ブロックで処理できるなら FFT の実装に切り替えます。

✅ 検証

本文の数値が正しいかを、別の計算と突き合わせます。

Goertzel と FFT

ランダムな 205 点の列について、205 本すべてのビンを Goertzel で計算し、NumPy の FFT と比べました。差の最大は 6.9 \times 10^{-12} （最大のビンの値に対して 2 \times 10^{-13} ）で、浮動小数点の丸めの範囲です。整数でない k = 17.86 の一般化 Goertzel も、和を直接計算した値と 1.3 \times 10^{-13} で一致しました。

整合フィルタの出力 SNR と ROC

出力 SNR は、3 つの形すべてで 0〜20 dB の全点が理論の E_s/\sigma^2 と ±0.1 dB で一致しました（2 万回ずつ）。ROC は、d = 1、2、3 と P_{FA} = 1e−4〜1e−1 の 12 点で、 P_D の計算値と理論の差がすべて 0.001 以内でした（40 万回ずつ）。

ヒルベルト FIR の長さと誤差

fs = 8 kHz、300〜3700 Hz の範囲で、振幅の誤差の最大を長さごとに求めました。

L = 63 からは、300 Hz で 1 % を下回ります。L = 127 から先の誤差は、ハミング窓の通過域のリプルで決まり、長さを延ばしてもゆっくりとしか減りません。解析信号の FFT の実装（analytic_fft）は、scipy.signal.hilbert と偶数長でも奇数長でも 2 \times 10^{-15} 以内で一致しました。

シードを変えても

乱数のシードを 3 通りに変えて、主な数値を計算し直しました。

P_D は d = 2、 P_{FA} = 10^{-2} の値です。どのシードでも結論は変わりません。

C と Python

実装の節の C（float）を MSVC でビルドし、同じ入力で Python と比べました。Goertzel の電力は 20 ブロックで相対誤差 1.6 \times 10^{-6} 、整合フィルタの出力は最大 6.4 \times 10^{-6} （ピーク 13.7 に対して）、ヒルベルトの虚部は 6.1 \times 10^{-7} で、単精度の丸めの範囲で一致しました。

チェックリスト

自分の設計を確かめるときの 6 項目です。

1. Goertzel の値が FFT と一致するか：同じ N、同じ k で、差が丸めの範囲か
2. トーンの周波数がビンに乗っているか： k = fN/f_s と最寄りの整数の差。0.5 に近いなら一般化 Goertzel か窓
3. 整合フィルタのピークが M−1 の位置に立つか：ノイズなしのテンプレートを通して確かめる
4. しきい値が誤警報の確率から決まっているか：信号のない区間で P_{FA} を数えて、狙いと合うか
5. ノイズが白色か：信号のない区間の PSD が平らか。平らでなければ白色化してから整合
6. ヒルベルトの通過域に信号の帯域が入っているか：下側帯波が FIR の 1 % 帯域の下端より上か

🕳️ 落とし穴

Goertzel の N と周波数が合わない

周波数がビンの中心から半ビンずれると、矩形窓では −3.9 dB 落ち、ほかのビンにも漏れます。k を周波数に合わせて非整数にする（一般化 Goertzel）か、ハン窓を掛けて −1.4 dB に抑えます。ハン窓は主ローブが 2 倍に太るので、隣のトーンとの間隔と相談して決めます。

固定小数点で係数と状態がはみ出す

係数 2\cos\theta_k は −2 から 2 の値をとるので、Q15（−1 から 1 未満を表す 16 ビットの固定小数点）には入りません。Q14（−2 から 2 未満を表す形式）で持つか、1 回の掛け算を「係数の半分を 2 回足す」形にします。係数の丸めは周波数のずれになり、Q14 の刻み 2^{-14} は、N = 205 の 697 Hz のビンで 0.07 Hz、100 Hz 付近のビンでは 0.4 Hz です（sin θ が小さいほど大きくずれる）。状態 v[n] は N に比例して伸び、振幅 1 の 697 Hz では最大 189、100 Hz では 857 になりました。入力の振幅に対して 8〜10 ビットの余裕を、32 ビットの累算で確保します。

整合フィルタの遅れ

整合フィルタのピークは、テンプレートの始まりではなく終わり（始まり＋M−1）に立ちます。到着時刻はピークの時刻から M−1 を引いて求めます。長いテンプレートほど判定が遅れるので、制御ループの中で使うなら M−1 サンプルの遅れを位相余裕で見積もります。

ノイズが白色でないと偽のピークが立つ

有色ノイズでは、白色を仮定した整合フィルタの出力に、ノイズの強い周波数の揺れがそのまま出て、正しくない時刻にピークが立ちます。設計例の AR(1) では、正しい時刻のピークは 7 % でした。信号のない区間で PSD を測り、白色化してから整合を取ります。

ヒルベルトの端の過渡

FFT の実装はブロックを周期的につないだものとして扱うので、ブロックの端で包絡線が数 % 乱れます。FIR の実装は、最初の 2D サンプルが過渡です。どちらも、端の区間は判定に使わないか、前後に余白を取ってから処理します。

瞬時周波数の負の値と巻き戻し

位相を引き算するときは、 \arg(z[n]\overline{z[n-1]}) で差を直接取り、±π の飛びを避けます。それでも負の値や帯域外の値が出たら、2 つ以上の成分が混じっている合図です。バンドパスで 1 成分に絞ってから瞬時周波数を取ります。

しきい値をノイズの推定なしに決める

P_{FA} の式は σ を知っている前提です。σ を決め打ちすると、ノイズが増えた日に誤警報が増えます。信号のない区間でノイズの分散を推定し続け、η を \sigma\sqrt{E_s}\,Q^{-1}(P_{FA}) に合わせて更新します。推定に使うサンプルが少ないと推定値自体がばらつくので、そのぶん η を少し上げます。

DC オフセットが Goertzel とヒルベルトに漏れる

ADC の入力には、たいてい直流のオフセットが乗っています。整数の k なら、直流は DFT の別のビン（k = 0）で、ちょうど 0 しか漏れません。一般化 Goertzel の非整数の k では漏れが出て、同じ振幅のトーンに対して 697 Hz（N = 275）で −49 dB、1633 Hz（N = 118）で −42 dB でした。オフセットがトーンの何十倍もあると、これが無視できなくなります。ヒルベルトでは、手計算 3 のとおり直流が包絡線を揺らします。どちらも、前段で #3 のハイパスを通して直流を落としておきます。

💻 実装

C：Goertzel・整合フィルタ・ヒルベルト FIR

#8 の LMS と同じく、状態を構造体に持って 1 サンプルずつ処理します。整合フィルタとヒルベルトの入力の履歴は、循環バッファ（書き込み位置を回して使い回す配列）に持ちます。

#ifndef DETECT_H
#define DETECT_H

/* ---- Goertzel：1 ビンを N サンプルごとに ---- */
typedef struct {
    float coeff;          /* 2cos(2πk/N) */
    float v1, v2;         /* v[n-1], v[n-2] */
    int   n, N;
} goertzel_t;

void goertzel_init(goertzel_t *g, float k, int N);
int  goertzel_push(goertzel_t *g, float x, float *power);  /* N 個目で 1 を返し |X[k]|^2 を書く */

/* ---- 整合フィルタ：テンプレートの長さ M ---- */
#define MF_M 80
typedef struct {
    float s[MF_M];        /* テンプレート s[0..M-1]（そのままの向きで持つ） */
    float x[MF_M];        /* 入力の循環バッファ */
    int   head;           /* 最新の入力の位置 */
} mf_t;

void  mf_init(mf_t *m, const float *s);
float mf_step(mf_t *m, float x);   /* y[n] = Σ s[l] x[n-M+1+l] */

/* ---- ヒルベルト FIR：長さ L = 2D+1（奇数）---- */
#define HB_D 31
#define HB_L (2 * HB_D + 1)
typedef struct {
    float g[HB_D + 1];    /* g[m]（m = 1, 3, 5, ... だけ使う。g[-m] = -g[m]） */
    float x[HB_L];        /* 入力の循環バッファ */
    int   head;
} hilbert_t;

void hilbert_init(hilbert_t *h);
void hilbert_step(hilbert_t *h, float x, float *re, float *im);  /* re = x[n-D]、im = x̂[n] */

#endif

#include <math.h>
#include <string.h>
#include "detect.h"

#define PI_F 3.14159265358979f

void goertzel_init(goertzel_t *g, float k, int N)
{
    g->coeff = 2.0f * cosf(2.0f * PI_F * k / (float)N);
    g->v1 = g->v2 = 0.0f;
    g->n = 0;
    g->N = N;
}

int goertzel_push(goertzel_t *g, float x, float *power)
{
    float v0 = x + g->coeff * g->v1 - g->v2;   /* 1 回の掛け算と 2 回の足し算 */
    g->v2 = g->v1;
    g->v1 = v0;
    if (++g->n < g->N) return 0;
    /* 終端処理：|X[k]|^2 = v1^2 + v2^2 - coeff v1 v2 */
    *power = g->v1 * g->v1 + g->v2 * g->v2 - g->coeff * g->v1 * g->v2;
    g->v1 = g->v2 = 0.0f;                      /* 次のブロックへ */
    g->n = 0;
    return 1;
}

void mf_init(mf_t *m, const float *s)
{
    memcpy(m->s, s, sizeof(m->s));
    memset(m->x, 0, sizeof(m->x));
    m->head = 0;
}

float mf_step(mf_t *m, float x)
{
    m->head = (m->head == 0) ? MF_M - 1 : m->head - 1;
    m->x[m->head] = x;                         /* x[head] が最新の x[n] */
    /* 最新から数えて M-1-l 個前の入力に s[l] を掛ける＝ h[i] = s[M-1-i] の FIR */
    float y = 0.0f;
    int i = m->head;
    for (int l = MF_M - 1; l >= 0; l--) {
        y += m->s[l] * m->x[i];
        if (++i == MF_M) i = 0;
    }
    return y;
}

void hilbert_init(hilbert_t *h)
{
    memset(h, 0, sizeof(*h));
    for (int m = 1; m <= HB_D; m += 2) {
        /* 理想応答 2/(πm) にハミング窓 0.54 + 0.46 cos(πm/D) */
        float w = 0.54f + 0.46f * cosf(PI_F * (float)m / (float)HB_D);
        h->g[m] = 2.0f / (PI_F * (float)m) * w;
    }
}

void hilbert_step(hilbert_t *h, float x, float *re, float *im)
{
    h->head = (h->head == 0) ? HB_L - 1 : h->head - 1;
    h->x[h->head] = x;
    /* x[n-i] は x[(head + i) % L]。中央は i = D */
    int c = h->head + HB_D;
    if (c >= HB_L) c -= HB_L;
    float acc = 0.0f;
    for (int m = 1; m <= HB_D; m += 2) {      /* 奇数の m だけ、対称性で掛け算を半分に */
        int a = c + m;  if (a >= HB_L) a -= HB_L;   /* x[n-D-m] */
        int b = c - m;  if (b < 0) b += HB_L;       /* x[n-D+m] */
        acc += h->g[m] * (h->x[a] - h->x[b]);
    }
    *re = h->x[c];
    *im = acc;
}

ヒルベルトの hilbert_step は、中央から ±m の 2 つの入力の差に 1 つの係数を掛けるので、L = 63 でも掛け算は 16 回です。ハミング窓の式 0.54 - 0.46\cos(2\pi i/(L-1)) は、中央からの距離 m で書くと 0.54 + 0.46\cos(\pi m/D) になります。整合フィルタは M 回の積和で、M が数百を超えるなら、入力を区切って FFT で畳み込み、結果をつなぎ合わせる重畳加算法のほうが軽くなります。

固定小数点で書くときは、落とし穴の節のとおり、Goertzel の係数を Q14 にし、状態を 32 ビットで持ちます。Arm の CMSIS-DSP には Goertzel の関数はありませんが、相互相関の arm_correlate_f32 と複素 FFT の arm_cfft_f32 があり、整合フィルタと FFT 版の解析信号に使えます。Espressif の ESP-DSP には、テンプレートとの相関の dsps_corr_f32、相互相関の dsps_ccorr_f32、基数 2 の FFT の dsps_fft2r_fc32 があります。

Python：7 つの関数

記事の図と数値は、次の関数で計算しました。SciPy の scipy.signal.hilbert と scipy.signal.correlate、NumPy の numpy.fft と同じ結果になることを、検証の節で確かめています。

import numpy as np


def goertzel(x, k, N=None):
    """x の先頭 N 点で、ビン k の DFT  X[k] = Σ x[n] e^{-j2πkn/N}  を返す。k は実数でもよい。"""
    x = np.asarray(x, float)
    N = len(x) if N is None else N
    theta = 2 * np.pi * k / N
    c = 2 * np.cos(theta)
    v1 = v2 = 0.0                       # v[n-1] と v[n-2]
    for n in range(N):
        v0 = x[n] + c * v1 - v2         # 実数の 2 次 IIR
        v2, v1 = v1, v0
    # 終端処理：y[N-1] = v[N-1] - e^{-jθ} v[N-2]、X = e^{-jθ(N-1)} y[N-1]
    return np.exp(-1j * theta * (N - 1)) * (v1 - np.exp(-1j * theta) * v2)


def goertzel_power(x, k, N=None):
    """|X[k]|^2 を実数だけで求める（終端処理の 3 回の掛け算）。"""
    x = np.asarray(x, float)
    N = len(x) if N is None else N
    c = 2 * np.cos(2 * np.pi * k / N)
    v1 = v2 = 0.0
    for n in range(N):
        v0 = x[n] + c * v1 - v2
        v2, v1 = v1, v0
    return v1 * v1 + v2 * v2 - c * v1 * v2


def matched_filter(x, s):
    """テンプレート s の整合フィルタ（因果の FIR、h[n] = s[M-1-n]）。
    出力 y[n] のピークは、テンプレートが x の n-(M-1) から始まる位置に立つ。"""
    h = np.asarray(s, float)[::-1]
    return np.convolve(np.asarray(x, float), h)[:len(x)]


def hilbert_fir(L):
    """長さ L（奇数）のヒルベルト変換器。理想応答 2/(πm)（m 奇数）にハミング窓を掛ける。遅れ D=(L-1)/2。"""
    assert L % 2 == 1
    D = (L - 1) // 2
    m = np.arange(-D, D + 1)
    h = np.zeros(L)
    odd = (m % 2) != 0
    h[odd] = 2.0 / (np.pi * m[odd])
    return h * np.hamming(L)


def analytic_fir(x, L):
    """FIR で解析信号 z = x + j x̂ を作る。x を D サンプル遅らせて x̂ と時刻をそろえる。"""
    h = hilbert_fir(L)
    D = (L - 1) // 2
    x = np.asarray(x, float)
    xh = np.convolve(x, h)[:len(x)]
    xd = np.concatenate([np.zeros(D), x[:len(x) - D]])
    return xd + 1j * xh


def analytic_fft(x):
    """FFT で解析信号を作る：負の周波数を 0、正の周波数を 2 倍（直流とナイキストはそのまま）。"""
    x = np.asarray(x, float)
    N = len(x)
    X = np.fft.fft(x)
    H = np.zeros(N)
    H[0] = 1.0
    if N % 2 == 0:
        H[N // 2] = 1.0
        H[1:N // 2] = 2.0
    else:
        H[1:(N + 1) // 2] = 2.0
    return np.fft.ifft(X * H)


def inst_freq(z, fs):
    """解析信号から瞬時周波数 [Hz]。位相の差を z[n] z*[n-1] の偏角で取るので、巻き戻しが要らない。"""
    dpsi = np.angle(z[1:] * np.conj(z[:-1]))
    return dpsi / (2 * np.pi) * fs

goertzel は、整数でない k でも位相まで正しい値を返すように、終端で e^{-j\theta(N-1)} を掛けています。整数の k なら、これは本文の e^{j\theta_k} と同じです。

手計算を確かめる

本文の 3 つの手計算を、上の関数で確かめるコードです。

import numpy as np
from scipy import signal

np.set_printoptions(precision=4, suppress=True)

# 手計算 1：Goertzel、N = 8、k = 1
x = np.array([0, 1, 1, 1, 0, -1, -1, -1], float)
print("Goertzel X[1] :", np.round(goertzel(x, 1), 5))
print("FFT      X[1] :", np.round(np.fft.fft(x)[1], 5))
print("|X[1]|^2      :", round(goertzel_power(x, 1), 4))

# 手計算 2：長さ 4 のテンプレート
s = [1, 1, -1, 1]
xr = [0, 0.2, -0.1, 1.1, 0.8, -0.9, 1.2, 0.1, -0.2, 0.0]
y = matched_filter(xr, s)
print("matched y     :", y, " peak n =", int(np.argmax(y)))
print("correlate     :", signal.correlate(xr, s, mode="full")[:len(xr)])

# 手計算 3：4 点の解析信号
for xx in ([1, 0, -1, 0], [3, 1, -1, 1]):
    z = analytic_fft(xx)
    print(xx, "z =", np.round(z, 4), " |z| =", np.round(np.abs(z), 4))
print("scipy hilbert :", np.round(signal.hilbert([3, 1, -1, 1]), 4))

出力は次のとおりです。

Goertzel X[1] : -4.82843j
FFT      X[1] : -4.82843j
|X[1]|^2      : 23.3137
matched y     : [ 0.   0.2 -0.3  1.4 -0.2 -0.7  4.  -1.2 -0.   1.5]  peak n = 6
correlate     : [ 0.   0.2 -0.3  1.4 -0.2 -0.7  4.  -1.2 -0.   1.5]
[1, 0, -1, 0] z = [ 1.+0.j  0.+1.j -1.+0.j  0.-1.j]  |z| = [1. 1. 1. 1.]
[3, 1, -1, 1] z = [ 3.+0.j  1.+2.j -1.+0.j  1.-2.j]  |z| = [3.     2.2361 1.     2.2361]
scipy hilbert : [ 3.+0.j  1.+2.j -1.-0.j  1.-2.j]

Goertzel の値と電力は表の手計算と、整合フィルタの出力は相関の表と、解析信号は 4 点の手計算と一致しています。signal.correlate の mode="full" の先頭の 10 個が matched_filter と同じになるのは、相関器と整合フィルタが同じ計算だからです。

💡 活用事例

• DTMF とトーンの検出：プッシュ音の判別は、Texas Instruments のアプリケーションノートが Goertzel の題材にしている代表的な用途。決まった周波数のトーンがあるかどうかを見る検出は、どれも同じ考え方で書ける
• レーダー、ソナー、超音波の距離の測定：送った波形をテンプレートにして整合フィルタで到着時刻を求める。チャープによるパルス圧縮で、送信の電力を上げずに分解能を上げる
• GNSS の受信機：衛星ごとに決まった符号（拡散符号）の複製を作り、受信信号との相関のピークから到着の時刻のずれを測る。これも既知の波形との相関による検出
• AM・FM・SSB の変復調：解析信号から包絡線（AM）と瞬時周波数（FM）を取り出す。片側の側帯波だけを作る SSB（単側波帯）にも、ヒルベルト変換を使う。Gabor の 1946 年の論文も、伝送の問題への応用の例として単側波帯の直接の生成を挙げている
• 心電図の R 波の検出：Pan と Tompkins の 1985 年のアルゴリズムは、専用のバンドパスで雑音を落とし、波形の傾き、振幅、幅で QRS 波（心電図の鋭い山）を見分け、しきい値を自動で調整する
• 機械の振動の包絡線：回転機械の振動を、共振の帯域でバンドパスしてから包絡線を取ると、ベアリングの傷のような周期的な衝撃の繰り返し周波数が見えやすくなる。包絡線解析と呼ばれる手法

🧭 使い分けの勘所：Goertzel・整合フィルタ・ヒルベルト

ここまでの結果を、判断の材料としてまとめます。方式を選ぶ前に、何を知っていて何が欲しいのかを確かめます。

判断の軸

方式を決める前に答えておく 7 つの問いです。上から順に、手元の検出の問題について確かめます。

⑥が抜けがちです。設計例 1 のとおり、固定の N のビンは Hz で等幅なので、% で決まった許容差は 1 つのしきい値で満たせないことがあります。

5 つの方式の比較

この記事の 3 つに、FFT と #4 のバンドパス＋検波を加えて並べます。まず周波数を測る側の 3 つです。

次に波形と包絡線の側の 2 つです。

決定木

同じ判断を 1 本の流れにします。上から順に答えていくと、方式と最初に決める値に着きます。

周波数も波形も分からず、帯域だけを知っていて位相が要らないなら、ヒルベルトより先に整流＋ローパスを試します。設計例 3 のとおり、包絡線だけなら誤差 0.05 程度で足りることが多く、計算も軽く済みます。

用途別の早見表

よくある場面ごとに、最初に試す方式と値の目安を並べます。

勘所

• ビンが数本なら Goertzel、十数本を超えたら FFT：数え方にもよるが、境目は 2\log_2 N 本前後。Goertzel はバッファがいらず、N を自由に選べる
• 周波数がビンに乗らないなら一般化 Goertzel：k を非整数にすれば、半ビンずれの −3.9 dB の損失が消える。許容差が % で決まるなら、窓に入る周期の数をそろえて N を変える
• テンプレートが分かるなら、整合フィルタに勝る線形フィルタはない：白色ガウスノイズのもとで出力 SNR の上限 E_s/\sigma^2 に届く。#4 のバンドパスは最良の Q 値でも 1.5 dB 劣った
• SNR が足りなければ、振幅ではなく長さで稼ぐ：出力 SNR はエネルギーで決まるので、M を 2 倍にすれば 3 dB 上がる。分解能を落としたくなければ、長さはそのままでチャープにする
• しきい値は誤警報の確率から決める： \eta = \sigma\sqrt{E_s}\,Q^{-1}(P_{FA}) 。σ は信号のない区間で推定し続ける
• ノイズが白色か、先に PSD を見る：有色なら白色化してから整合。そのままでは偽のピークが立つ
• 包絡線は整流＋ローパスで足りることが多い：位相が要るとき、変調が速いとき、瞬時周波数が要るときにヒルベルト
• 瞬時周波数は 1 成分の狭帯域のときだけ：2 つの成分が混じると、帯域の外や負の値になる。先にバンドパスで 1 成分に絞る

✅ まとめ

• 検出のフィルタは、周波数や波形という仮説を前もって与え、入力から「あるか・いつか・どれだけか」を判定量にして返す。周波数を知っているなら Goertzel、波形を知っているなら整合フィルタ、帯域だけを知っていて振幅と位相が欲しいならヒルベルト変換
• Goertzel は、DFT の 1 ビンを実数の 2 次 IIR v[n] = x[n] + 2\cos\theta_k v[n-1] - v[n-2] で計算し、最後に 3 回の掛け算で電力を出す。極は単位円の上で、N で止めるから使える。k を非整数にした一般化 Goertzel で、ビンに乗らない周波数をそのまま測れる
• ビンの応答は sinc 形で、第 1 サイドローブは矩形窓で −13 dB、ハン窓で −31 dB。半ビンずれると −3.9 dB 落ち、全ビンに漏れる
• 整合フィルタは、テンプレートを時間反転したインパルス応答を持つ FIR で、コーシー・シュワルツの不等式から出力 SNR の上限 E_s/\sigma^2 に届く。SNR は形によらずエネルギーで決まり、ROC は P_D = Q(Q^{-1}(P_{FA}) - d) 。有色ノイズでは白色化してから整合を取る
• 整合フィルタは群遅延を一定にしてすべての成分を 1 時刻に集める。チャープを使うと、同じ長さとエネルギーで頂上の幅が 202.5 μs から 62.5 μs に縮み、100 μs 離れたエコーを分けられた
• ヒルベルト変換は全周波数で −90° 回す。インパルス応答は無限長で因果でないので、窓で切って D 遅らせた FIR で近似する。位相は −90° ちょうどで、誤差は振幅の両端に出る。L = 63 なら fs = 8 kHz で 207 Hz から 1 % 以内
• 解析信号から包絡線と瞬時周波数を取れるが、瞬時周波数は 1 成分の狭帯域のときだけ意味を持ち、2 成分では負にもなる。直流も包絡線を揺らす
• 設計例では、DTMF で固定の N = 205 が % の許容差を満たせず（±1.5 % の受け入れ 43 %）、窓に 24 周期の一般化 Goertzel で 100 % と誤受け入れ 0 を両立した。超音波の到着時刻では、SNR 0 dB で整合フィルタの外れ率が 2.2 %（バンドパス＋しきい値は 49 %）、チャープはさらに約 3.6 dB 低い SNR まで外れなかった。AM の包絡線は整流＋ローパスで誤差 0.049、ヒルベルトで 0.020

ここまでで、知っているもの（周波数・波形・帯域）から方式を選び、N、テンプレート、しきい値、L を式と数値から決め、FFT との一致、出力 SNR、ROC、包絡線の誤差で判定できるようになりました。次の一手は、手元で鳴らせるトーン（Arduino の tone 関数など）を ADC で取り込み、この記事の goertzel_power を 1 本ずつ当てて、周波数を少しずつずらしたときの電力が (b) の曲線のように落ちるかを確かめることです。

次回の #11 は「レート変換」です。サンプリング周波数を下げたり上げたりするとき、エイリアシングを防ぎながら軽く計算する CIC フィルタとコムフィルタ、位相だけを変えるオールパスフィルタを、同じ流儀で式から設計します。この記事のヒルベルト変換をオールパスの組で作る方法も、そこで扱います。

よくある質問（FAQ）

Q. Goertzel と FFT は、どちらを使えばいいですか？

A. 必要なビンが数本なら Goertzel、十数本を超えるなら FFT です。掛け算の数を実数で数えると、Goertzel は N_b 本で約 N_b N 回、基数 2 の FFT は約 2N\log_2 N 回なので、N = 256 なら 16 本前後が境目です。Goertzel は 1 サンプルずつ処理できて N を自由に選べ、FFT は全体を一度に見られます。DTMF のように 8 本だけなら Goertzel が軽く済みます。

Q. Goertzel の N はどう決めますか？

A. 3 つの条件から決めます。ビンの間隔 f_s/N を、見分けたい 2 つの周波数の間隔より十分狭くすること。1 ブロックの時間 N/f_s が、信号が続く時間より短いこと。そして、測りたい周波数ができるだけビンの中心に乗ることです。乗らないなら、一般化 Goertzel で k を非整数にします。許容差が % で決まっているなら、窓に入る周期の数をそろえて、周波数ごとに N を変えます。

Q. 整合フィルタのテンプレートが少し違っていたら、どれくらい損をしますか？

A. 周波数のずれは、テンプレートの長さの間に何周期分の位相がずれるかで決まります。8 周期のバーストでは、2 % のずれで −1.5 dB、5 % で −12.8 dB でした。位相が 45° ずれると −3 dB、90° ずれると内積が 0 になります。位相が分からない受信では、sin と cos の 2 つのテンプレートで整合を取り、出力の 2 乗和の平方根（包絡線）で判定します。

Q. 包絡線を取るのに、ヒルベルト変換は必要ですか？

A. 包絡線だけなら、多くの場合は全波整流＋ローパスで足ります。1 kHz の搬送波を 20 Hz で変調した例では、整流＋ローパスの誤差は 0.049、ヒルベルトは 0.020 でした。整流側の誤差のほとんどは、1 周期 8 サンプルと粗いことによる約 5 % の目盛りのずれで、fs を上げれば小さくなります。ヒルベルトが要るのは、位相や瞬時周波数も欲しいとき、変調が速くてローパスのリプルと遅れが釣り合わないときです。

Q. 瞬時周波数が負になりました。計算の誤りですか？

A. 計算の誤りではなく、信号に 2 つ以上の成分が混じっている合図であることが多いです。1000 Hz と、その 0.9 倍の振幅の 1150 Hz の和では、理論で −350 Hz まで落ちます。瞬時周波数は 1 成分の狭帯域の信号でだけ直感どおりの値になるので、先にバンドパスで 1 成分に絞ります。直流が残っていても位相がおかしくなるので、ハイパスで落としておきます。

関連記事

この連載：

• #1 基礎編 — 極と零点、群遅延、白色ノイズの PSD と Welch 法
• #2 ローパスフィルタの設計 — 整流＋ローパスの包絡線と、制御ループの位相余裕
• #3 ハイパスフィルタの設計 — Goertzel とヒルベルトの前段で直流を落とす
• #4 バンドパスとノッチの設計 — 整合フィルタと比べたバンドパスと Q 値
• #5 カルマンフィルタの設計 — モデルと分散で切り分ける
• #6 EKF・UKF の設計 — 角度の差の巻き戻し
• #7 パーティクルフィルタの設計 — 分布をサンプルで持つ
• #8 適応フィルタの設計 — 前白色化と、線スペクトルを取り出す ALE
• #9 平滑化 — 外れ値を落とし、形を保つ。位相の微分に使える SG 微分
• #11 レート変換のフィルタの設計 — CIC・コム・オールパス。ヒルベルト変換をオールパスの組で作る

信号の入口：

• Arduino の AD 変換（ADC） — 検出の前にサンプルを取り込む
• Arduino の tone 関数と圧電ブザー — Goertzel で測るトーンを鳴らす
• 加速度センサ KXR94-2050 で作るデジタル水平器 — センサの値を読んで処理する題材

参考

取得日はいずれも 2026-11-10 です。図と数値はすべて筆者が Python（NumPy、SciPy、Matplotlib）で計算したものです。乱数のシードは固定です。

Goertzel アルゴリズムと DTMF

• G. Goertzel「An Algorithm for the Evaluation of Finite Trigonometric Series」（The American Mathematical Monthly, 65(1), 34–35, 1958）
• G. Schmer「DTMF Tone Generation and Detection: An Implementation Using the TMS320C54x」（Texas Instruments, SPRA096A, 2000）
• C. J. Chen「Modified Goertzel Algorithm in DTMF Detection Using the TMS320C80 DSP」（Texas Instruments, SPRA066, 1996）
• ITU-T Recommendation Q.23「Technical features of push-button telephone sets」、Q.24「Multifrequency push-button signal reception」

整合フィルタと検出

• D. O. North「An Analysis of the factors which determine signal/noise discrimination in pulsed-carrier systems」（Proceedings of the IEEE, 51(7), 1016–1027, 1963。1943 年の RCA 研究所の報告の再録）
• J. H. Van Vleck, D. Middleton「A Theoretical Comparison of the Visual, Aural, and Meter Reception of Pulsed Signals in the Presence of Noise」（Journal of Applied Physics, 17(11), 940–971, 1946）
• G. L. Turin「An introduction to matched filters」（IRE Transactions on Information Theory, 6(3), 311–329, 1960）

ヒルベルト変換と瞬時周波数

• D. Gabor「Theory of communication. Part 1: The analysis of information」（Journal of the IEE, Part III, 93(26), 429–441, 1946）
• B. Boashash「Estimating and interpreting the instantaneous frequency of a signal. I. Fundamentals」（Proceedings of the IEEE, 80(4), 520–538, 1992）

活用事例

• J. Pan, W. J. Tompkins「A Real-Time QRS Detection Algorithm」（IEEE Transactions on Biomedical Engineering, BME-32(3), 230–236, 1985）

教科書

• A. V. Oppenheim, R. W. Schafer「Discrete-Time Signal Processing」（Pearson）
• S. M. Kay「Fundamentals of Statistical Signal Processing, Volume II: Detection Theory」（Prentice Hall）
• S. L. Hahn「Hilbert Transforms in Signal Processing」（Artech House）

ライブラリ

• SciPy「scipy.signal.hilbert」
• SciPy「scipy.signal.correlate」
• NumPy「numpy.fft」
• Arm CMSIS-DSP「arm_correlate_f32.c」
• Arm CMSIS-DSP「arm_cfft_f32.c」
• Espressif ESP-DSP「dsps_corr.h」
• Espressif ESP-DSP「dsps_fft2r.h」