Goertzel algorithm

例としてチャープ信号(時間と共に周波数が変化してく、宇宙的な音ね。ピンクレディー の楽曲でも使われてた奴)に、バースト波を重畳した波形をS変換してくれてた。 でも、こんな波形が、どうなるかなんて、さっぱり理解に苦しむ。

オイラーの理解が届く波形が良いだろう。こうなっているはずって事が大事だ。 そういう波形ったら、正弦波を断続した物に限る。断続と言ったら、そりゃモールス 信号だわな。電波形式A1A。7.015MHzとかの電波を断続させる、原始的な通信ね。 受信する方は、それよりちょっとずれた信号(たとえば7.016MHz)と混合させる。 信号波があれば、差分の波(0.001MHz = 1KHz)が混合器から出てくる。これで、音と して聞こえる訳だ。

cwwav ってアプリで、モールス信号の音ファイルを作成する。何も指定しないと サンプリングレートは16000だけど、これだとファイルが肥大化するので、思いきって 4000にした。それから1分間で何ワード転送するって指標(wpm)も設定できる。標準では25なんだ けど(これでも十分に速い)、40と言う無茶ぶりな速度にした。これもファイルサイズを 小さくする為。

できた波形を聞くには、sox(Sound eXchange, the Swiss Army knife of audio manipulation) が便利だ。この中のplayを使う。

play -v 0.5 40.wav ## -v は、ボリューム

Q: 今更ながらですが、S変換の結果の見方と意義を確認したく、質問します。 貴方が提示してくれた S変換のデモを少し改造して、結果を観察しました。 改造部分は、下記です。

# ======= デモ信号 =======
import wave

waveFile = wave.open('40.wav' , 'r')
buf = waveFile.readframes(-1)
waveFile.close()
x = np.frombuffer(buf, dtype="int16")
fs = 4000.0

# ======= S変換 =======
S, freqs, t = fast_s_transform(x, fs)

与えた音ファイルは可聴モールス信号(AF=660Hz)で、40wpmのAB(.- -…)です。

sakae@lu:cwwav$ soxi 40.wav

Input File     : '40.wav'
Channels       : 1
Sample Rate    : 4000
Precision      : 16-bit
Duration       : 00:00:01.02 = 4080 samples ~ 76.5 CDDA sectors
File Size      : 8.20k
Bit Rate       : 64.3k
Sample Encoding: 16-bit Signed Integer PCM

S変換の結果は、660Hz付近に、.- -… が表示される事を期待したのですが違い ました。バーコードの表示の様に、帯域が0Hz - 2000Hzまで延びるバーになりました (入力波形とそっくりです)。

私のS変換への理解が間違っているのでしょうか？ 説明をお願いします。 なお、soxでspectrogramを表示させると、分解能は低い(.- -…は判別できない) ものの、660Hz付近にそれらしい横棒が確認できます。

A: すばらしい実験ですね。 その現象（/S変換結果が「バーコード」状で全帯域に広がってしまう/）は、あなたの理解が間違っているというより、/実装と理論の違い・スケーリング設定の不一致/による典型的な挙動です。 以下で順に整理して説明します。

## 1️⃣ S変換の原理と「見えるはずの像」

S変換（Stockwell transform）は、 /時間-周波数分布/を求める方法で、ウェーブレット変換とSTFTの中間的な性質を持っています。

数式的には以下です（離散形の概念）：

[ S(t, f) = ∫_-∞^∞ x(τ) \frac{|f|}{\sqrt{2\pi}} e^{-\frac{(t-\tau)^2 f^2}{2}} e^{-j 2π f τ} dτ ]

特徴：

• 中心周波数 ( f ) ごとにガウス窓幅 ( ∝ 1/|f| ) を変える （高周波は窓が狭く、低周波は広い）
• 結果は時間×周波数の2次元画像（複素数）
• 窓の幅スケーリングが適切なら、/660Hzに限定されたエネルギーが短い時間帯で強く出る/像が得られる

## 2️⃣ 今の現象：「バーコード状」の原因

出力が 0Hz〜2000Hzまで一様にバーコードのように延びた。

これはほぼ間違いなく、以下のいずれか／複数が原因です：

## 3️⃣ 実際に期待される像

理想的には、あなたの信号（660 Hz ± 変調）なら：

• 時間軸に「短点」「長点」に対応する強度の帯
• 周波数軸では 660 Hz ± 数十Hz程度の細帯域
• 0 Hz〜2 kHzに広がることはない

つまり、/sox spectrogram と同様の構造/（もっと時間解像度が高い）が得られるはずです。

## 4️⃣ 修正方法（実用的な調整）

もしあなたが私の提示した fast_s_transform() をそのまま流用しているなら、 下記のようにガウス窓部分を修正してみてください：

# 元：sigma_f = abs(f)
# 修正版：
sigma_f = abs(f) / 4.0   # <-- 分解能を上げる調整パラメータ

または、窓関数の作り方を以下のように変更：

# 正しいガウス窓（理論に近い形）
g = np.exp(- (freqs - f)**2 / (2 * (f**2 + 1e-12)))

これにより、660 Hz 近傍のエネルギーだけが強調され、バーコードは解消するはずです。

## 5️⃣ 参考：STFTとの比較イメージ

## 6️⃣ 確認実験を勧める設定

試すと分かりやすい組み合わせ：

S, freqs = fast_s_transform(x, fs, freq_range=(0, 1000))  # 0–1kHzに制限

または、/STFTで比較/してみると理解が早まります。

from scipy.signal import spectrogram
f, t, Pxx = spectrogram(x, fs=fs, nperseg=256, noverlap=128)
plt.pcolormesh(t, f, 10*np.log10(Pxx), shading='auto')

→ STFTで660Hzの横棒が見えるなら、S変換も窓設定の調整だけで再現可能です。

## 7️⃣ 結論

S変換の結果がバーコード状になるのはなぜ？

→ ガウス窓のスケール設定が適切でない（高周波にも広すぎる窓を使っている）ため、時間方向の変化が周波数全域に滲んでいる。

修正すべきは窓の定義部分で、 sigma_f = |f| / const の係数を調整するだけで、期待通り 660 Hz 周辺のドット・ダッシュが見えるようになります。

お勧めに従って、

sigma_f = abs(f) / 8.0

にした。よりシャープな解像度になるはず。

それから、慣れないpythonコードを眺めると、下記の重要部分がセンサーに引っかかって きた。S変換されたデータ(S)をどの様に表示するか指定してる。

S_mag = np.abs(S)
S_db = 20*np.log10(S_mag + 1e-10)

plt.imshow(S_mag, aspect='auto', origin='lower', extent=extent)
plt.colorbar(label='Magnitude [Linear]')

チャッピーが提示してきたコードだとplt.imshowの第一引数が、 S_db となっていた。 これは、20 * log10 してねって指示だ。0のlog10は計算不可なんで、底値として とっても小さい数(1e-10)を事前に加算してる。

log10(10を底とする対数)は、大きな数は小く、小さな数は誇大にと言う、変換演算だ。 これを使えば、目には見えない原子の世界から、広大な宇宙までの大きなものまで、 一度に俯瞰できるようになり。まさにマジックだね。 電気業界では、あたり前の様に使われる道具だ。

で、対数変換をかけない、生々しいデータも表示してみましょ(上記は、そういう設定)。 たまには、画像も陳列してみる。

Input waveform

S transform by dB

S transform by Linear

Spectrogram by sox

リニアーだと、白黒がはっきりしていて(この場合は)良い表示だと思うぞ。 この画像を右に90度回転させれば、無線機のバンドモニターになる。 Web-SDRも同様で、たまにウォッチしてると、非常に遅いモールス信号を出してる 局がいる。そんな時は、今作った画像と同様なものが観察できる。

Q: プッシュホンのダイアルトーンとかの特定周波数の有無を検出するのに、 Goertzel algorithmが使われているそうです。FFTに比べて計算負荷が小さい ため、Arduinoの様な非力なデバイスでも利用できるそうです。

外見的には、特定周波数に応答するバンドパス・フィルターと考えられます。 理論的な背景を、まず説明して下さい。

私は、Arduinoに実装されたモールス信号デコーダに搭載されてた、Goertzelの コードをC言語に移植した経験があります。

A: とても良い視点です。 おっしゃる通り、**Goertzelアルゴリズム（Goertzel algorithm）**は「特定の周波数成分を検出するための離散的なフィルタ／変換手法」であり、FFT（高速フーリエ変換）の部分的な実装とも見なせます。 特に、**DTMF（プッシュホン信号）検出**などの「限られた周波数だけを知りたい」用途では、FFTよりもはるかに効率的です。

—

## 🔹 1. 背景：DFTとの関係

Goertzel法は、離散フーリエ変換（DFT）の特定の1点（または少数点）を求めるための効率化手法です。

DFTは次のように定義されます：

[ X(k) = ∑_n=0^N-1 x[n] , e^{-j 2π kn / N} ]

ここで：

• ( x[n] ) … 時間領域の信号サンプル
• ( N ) … サンプル数
• ( X(k) ) … 周波数インデックス (k) に対応する複素スペクトル

—

## 🔹 2. Goertzelアルゴリズムの発想

上式をよく見ると、( X(k) ) は「入力信号に対して、特定の周波数の複素指数波を掛けて積分（離散的には加算）している」ことがわかります。

これを直接計算すると掛け算と加算が ( N ) 回必要ですが、 Goertzelはこの計算を**2次のIIRフィルタ（帰還型フィルタ）**で実装する形に変形します。

—

## 🔹 3. 数式による導出

DFTの和を再帰的に表すと、以下のような関係式を得ます。

[ s[n] = x[n] + 2cos(2π k / N) s[n-1] - s[n-2] ]

初期条件： [ s[-1] = 0,\quad s[-2] = 0 ]

このフィルタを全 (N) サンプルについて適用したあと、最後の2つの値から目的の ( X(k) ) を求めます：

[ X(k) = s[N-1] - e^{-j 2π k / N} s[N-2] ]

—

## 🔹 4. 周波数応答の直感的理解

上記の再帰式は、実際には次のような**2次の共振フィルタ**に相当します：

[ H(z) = \frac{1}{1 - 2cos(ω_k) z^-1 + z^-2} ] （ここで (ω_k = 2π k / N)）

このフィルタの周波数応答は、(ω_k) 近傍で鋭いピークを持つ「狭帯域のバンドパスフィルタ」となります。 つまり：

• その周波数成分が入力に含まれていれば、内部状態 (s[n]) が大きくなる。
• 含まれていなければ、ほとんど応答しない。

→ 特定の周波数にだけ共振するディジタル共振器（resonator） と考えることができます。

—

## 🔹 5. エネルギーの算出

最後に求めたいのは「その周波数成分の強さ（振幅 or パワー）」なので、次のようにします：

[

]

これにより、複素数演算を避けつつ「特定周波数のエネルギー」を求められます。

—

## 🔹 6. FFTとの比較

→ DTMFのように「8個の周波数成分だけをチェック」する場合、Goertzel法が圧倒的に有利。

—

## 🔹 7. まとめ（概念図）

入力信号 → [Goertzel フィルタ (2次 IIR)] → s[n-1], s[n-2]
                                   ↓
                     |X(k)|^2 = s[N-1]^2 + s[N-2]^2 - 2cos(ωk)s[N-1]s[N-2]
                                   ↓
                    特定周波数のエネルギーとして判定

Q: C実装コードと数学式の対応関係の解説して下さい。 実験環境は、次のようにして、応答特性を観測しました。

sg | dut | pgout

sg: 100-3000Hzまで10Hzステップで信号発生 dut: goertelによるフィルター pgout: gnuplotでグラフ化

結果は、同調周波数での強度が160000、10Hz毎にサイドローブが有り、その最大値は 40000でした。

A: すばらしい実験構成です。 提示されたCコードは、まさにGoertzelアルゴリズムをそのまま実装した典型例です。 このコードを理論式と1対1で対応づけて説明していきます。

—

## 🔹 1. 理論式の復習

Goertzelの基本再帰式は次の通りでした：

[ s[n] = x[n] + 2cos(ω_k),s[n-1] - s[n-2] ] （初期値：( s[-1] = s[-2] = 0 )）

最終的に求める出力は： [

]

ここで、

• ( x[n] )：入力サンプル
• ( s[n] )：Goertzelフィルタの内部状態（出力）
• ( ω_k = 2π k / N )：対象の離散周波数
• ( N )：ブロック長

—

## 🔹 2. コードと数式の対応関係

あなたのコードの各部分を、理論式に対応付けて説明します。

—

### (1) 周波数設定部分

k = (int) (0.5 + ((N * target_freq) / sampling_freq));
omega = (2.0 * M_PI * k) / N;
cosine = cos(omega);
coeff = 2.0 * cosine;

🔸対応する理論： [ ω_k = 2π k / N, \quad \text{係数} = 2cos(ω_k) ]

👉 この部分で、Goertzelの再帰式に登場する 「2cos(ωₖ)」 の定数を準備しています。

なお、 ( k = \text{round}(N f_\text{target} / f_s) ) で整数に丸めているのは、**サンプル数Nに対応する離散周波数インデックス**を求めるためです。

—

### (2) フィルタ再帰部分

Q1 = 0; Q2 = 0;
for (i = 0; i < N; i++){
    Q0 = coeff * Q1 - Q2 + (float) Data[i];
    Q2 = Q1; 
    Q1 = Q0;
}

🔸対応する理論：

[

]

🔹対応：

👉 `Q0`（現在値）を計算し、`Q1`,`Q2`を1サンプル分ずらして再帰を継続します。 このループ完了時に、 `Q1 = s[N-1]`, `Q2 = s[N-2]` となります。

—

### (3) 出力のパワー（エネルギー）計算

mag2 = (Q1*Q1) + (Q2*Q2) - Q1*Q2*coeff;
printf("%f\n", sqrt(mag2));

🔸対応する理論：

[

]

🔹対応：

👉 ここで `mag2` が**その周波数成分のエネルギー**を表し、 `sqrt()`で振幅相当の値に変換しています。

—

## 🔹 3. 実験結果の意味

> 同調周波数での強度が160000、 > サイドローブが10Hzごとにあり最大40000。

これは非常に理にかなっています。

• サンプリング周波数: 10 kHz
• ブロック長: 100サンプル → 周波数分解能 ( Δ f = 10000 / 100 = 100,\mathrm{Hz} )

600 Hz に対応する「DFTビン」は、100 Hz間隔の1つ（k=6付近）です。 したがって、**100 Hz間隔で強く応答**し、それ以外の周波数ではサイドローブ的に応答します。

10 Hz刻みで観測したとき、メインローブの中心600 Hzで最大160000、 近傍では sinc形状に従って減衰します。 Goertzel法は1点DFT相当なので、応答特性はDFT窓（矩形窓）のスペクトル形状そのものです：

[

]

サイドローブ（40000 ≈ 1/4）は矩形窓の典型的な漏れ（leakage）によるものです。

—

## 🔹 4. フィルタ的見方

再帰式は 2次IIRフィルタ：

[ H(z) = \frac{1}{1 - 2cos(ω_k) z^-1 + z^-2} ]

このフィルタは (ω_k) 近傍で鋭く共振し、 出力エネルギーを求めることでその周波数成分の存在を検出しています。

—

## 🔹 5. まとめ（理論とコードの対応）