2011-10-26

最初に、この翻訳は前半部分しかありません。後半の電子工作の話は、「今年こそArduinoとかやろぞ〜」といって数年経った人間ですので、訳語とかで迷惑をかけること必至で訳しませんでした。さて内容はプッシュ式電話機などのボタンを押したときに発せられるピ、ポ、パ音の周波数を検出するというもの。毎度ですが訳の精度は高くないので怪しい時は原文でご確認ください。おまけ、実際にコードを書いてみたい人はここを参考にされてはどうでしょう。
>>>それでは、はりきってどうぞ。

Goertzel algorithm DTMF Detection(日本語訳)

原文 http://bit.ly/vcjM23

GoertzelアルゴリズムはDSP ドメイン(訳注：時間ドメインや周波数ドメインなど使って解析する分野)でよくしられたものだ。それは周波数の存在やそれらがいくつあるのかを検出するのに使われる。少ない周波数を検出するつもりなら、DFT（FFT）のような他のテクニックに比べて実装するのに必要な計算が少なくてすむので、このアルゴリズムは人気がある。

わたしが今まで見たことがあるGoertzelアルゴリズムを実装しているプロジェクトをあげてみる:-
1.DTMF検出
2.ギターの弦のチューニング
3.オーディオ信号の特定の周波数を検出する

私のプロジェクトは、DTMFの周波数を検出するためにPIC18F4520とGoertzelアルゴリズムを使う。コードはすべてC言語で書かれている。たとえ私のプロジェクトの実施が良好だとしても、それが商用アプリケーションで使用できることを保証するものではない。

Goertzelアルゴリズムを理解する

Goertzelアルゴリズムの理解にはz平面上で単位円の線図(下図参照)を使って行うことができる。

DTMF検出には普通、サンプルリング周波数は8000Hzにするように選択される。しかし、私はサンプリング周波数を5908Hzとして使用するよう選択した。この特定の値を選択した理由については後述する。

我々が単位円上の2つの共役極（効果的に複雑な項を削除し、現実世界のマイクロコントローラを使用して実装を可能にするための共役極）を持つシステムを想定した場合、そのようなシステムの出力は、その周波数で鋭いピークを持つ。

5908Hzのサンプリング周波数は、z平面の単位円上の2πラジアンに当てる。また、πを超える周波数はサンプリング周波数の2分の1より多くなり、したがって、その後のナイキストレートを超えて当てた、そこを考慮しなくていい。

今、我々はDTMF上に次の周波数を持っている：-

697Hz
770Hz
852Hz
941Hz

1209Hz
1336Hz
1477Hz

z平面上に697Hzを配置するために、我々は、サンプリング周波数でこれを正規化しよう:-

5908Hz ---> 2π
だから 697Hz = (2π/5908) * 697 = 0.74126ラジアン

同様の変換に続いて、各周波数の値は次のとおりです:-

697Hz = 0.74126ラジアン
770Hz = 0.81889ラジアン
852Hz = 0.9061ラジアン
941Hz = 1.00713ラジアン

1209Hz = 1.28577ラジアン
1336Hz = 1.4208ラジアン
1477Hz = 1.57079ラジアン

我々は697Hzの周波数でのシステムのために2つの共役極、そして原点に2のゼロを配置する場合、z平面の図はこのようになる：-

だから、システムの伝達関数は次のようになる:-

$\frac{(z - 0)(z - 0)}{(z - e^{j0.74126})(z - e^{-j0.74126})}$

私たちの周波数を一般化するために、項0.74126を"θ"に置き換えよう：-

$\frac{(z - 0)(z - 0)}{(z - e^{j\theta})(z - e^{-\j\theta})}$

$\frac{z^2}{z^2 - ze^{-j\theta} - ze^{j\theta} + 1}$

次にオイラーの公式を使って代入

${e^{j\theta} = Cos\theta + jSin\theta}$

そして

${e^{-j\theta} = Cos\theta - jSin\theta}$

等式は以下のようになる:-

$\frac{z^2}{z^2 - z(Cos\theta - jSin\theta) - z(Cos\theta + jSin\theta) + 1}$

$\frac{z^2}{z^2 - zCos\theta - zjSin\theta - zCos\theta + zjSin\theta + 1}$

${ \frac{z^2}{z^2 - 2zCos\theta + 1} = \frac{Output}{Input}}$

"Z^-2"で分子と分母を乗算：-

${ \frac{z^2 * z^{-2}}{(z^2 - 2zCos\theta + 1) * z^{-2}} = \frac{Output}{Input}}$

${ \frac{1}{1 - 2z^{-1}Cos\theta + z^{-2}} = \frac{Output}{Input}}$

"Z^-1"を1つ目の遅延として、そして"Z-2"を2つ目の遅延に代入：-

Input = Output  - 2 * Cosθ * Prev_Output + Prev_Prev_Output
Output = Input  + 2 * Cosθ * Prev_Output - Prev_Prev_Output

上記の式は、Goertzelアルゴリズムの実装です。項（2 *cosθ）は"係数"と呼ばれ、各周波数に対して事前に計算されます： -

697Hz = 0.74126ラジアン = 1.4752
770Hz = 0.81889ラジアン = 1.3660
852Hz = 0.9061ラジアン = 1.2336
941Hz = 1.00713ラジアン = 1.0685

1209Hz = 1.28577ラジアン = 0.5623
1336Hz = 1.4208ラジアン = 0.2988
1477Hz = 1.57079ラジアン = 0

"1477Hz"の係数がゼロであることから確認できるように、基本的には1477Hzの検出時の計算を大幅に節約する。このような所以で私は5908Hzのサンプリング周波数を選択した。

これらの係数は、私のコードの"DTMF.h"ファイルに"#define"で定数として配置されている。 PICの計算をスピードアップするために、私は浮動小数点計算を避ける。したがって、それを整数に変換するために係数は256で乗算されている。後で、ほんとうの値を知るために結果を256で割る。この手で、私の計算は浮動小数点の実装よりも早い手法だといえる。

Goertzelの計算は、多くの取得したADC（デジタル変換回路）のサンプルで実行される。私は計算を実行するために100のサンプル数を選択している。これは周波数の検出分解能と処理時間のバランス。

計算は100のサンプルすべてに対して実行され、検出ピークの大きさは次のように計算される:-

Magnitude = Prev_Prev_Output 2 + Prev_Output2 - (Coefficient * Prev_Prev_Output * Prev_Output)

この大きさが閾値より>なら、その時、周波数が存在する。そうでなければ、周波数は存在しない。

if(magnitude > THRESHOLD)
{
// 周波数が存在する
}

閉村観光

GoertzelアルゴリズムでDTMF検出