Arduino Goetzel アルゴリズム 周波数検出器

　Raspberry Pi PicoでCW解読器（モールス符号解読器）を作ろうと思っています。CW解読器として有名なのは、Arduinoを使ったOZ1JHMのCWデコーダです。Arduino言語（C言語）で書かれているプログラムをPi用のmicropythonに移植する作業にとりかかりました。

　核になる部分を移植したのですが、Picoでは、うまく動きませんでした。やはりOZ1JHM CW decoderの動作原理を理解したうえで移植しないといけないと思い、Arduinoプログラムの解析にかかりました。

　OZ1JHM CW decoderの特徴は、アナログ入力にCWトーン信号を入れて、特定の周波数が入力されたときにマーク、それ以外（信号がない時も含め）の時にスペースと判断するプログラムが実装されていることです。

　これを行っているのが、Goetzel(ゲーツェル) アルゴリズムを使った周波数検出器です。FFTよりも簡易的に高速に特定の周波数を検出できるアルゴリズムとのことですので、その特性を調べてみることにしました。Goetzel アルゴリズムの説明はこちら。（英文のpdfです）

　回路図です。アナログ入力A1に信号を入れ、検出された大きさ（magnitude）をシリアルモニタに出力します。

　スケッチです。OZ1JHM CW decoderのスケッチから周波数検出部だけを抜き出しました。

-----------------------------------------------------------------------------------------

int audioInPin = A1;

float magnitude ;

float coeff;

float Q1 = 0;

float Q2 = 0;

float sine;

float cosine;  

float sampling_freq=8928.0;

float target_freq=558.0; /// adjust for your needs see above

float n=48.0;  //// if you change  her please change next line also 

int testData[48];

void setup() {

   int  k;

   float omega;

  k = (int) (0.5 + ((n * target_freq) / sampling_freq));

   omega = (2.0 * PI * k) / n;

   sine = sin(omega);

   cosine = cos(omega);

   coeff = 2.0 * cosine;

   Serial.begin(115200);

}

void loop() {

for (char index = 0; index < n; index++)

  {

     testData[index] = analogRead(audioInPin);

  }

  for (char index = 0; index < n; index++){

     float Q0;

     Q0 = coeff * Q1 - Q2 + (float) testData[index];

     Q2 = Q1;

     Q1 = Q0;  

  }

 float magnitudeSquared = (Q1*Q1)+(Q2*Q2)-Q1*Q2*coeff;  // we do only need the real part //

   magnitude = sqrt(magnitudeSquared);

  Q2 = 0;

  Q1 = 0;

   Serial.println(magnitude);

   delay(1000);

}

-----------------------------------------------------------------------------------------

　1secごとに各周波数（100Hz～1100Hz）におけるmagnitudeの値を読み取りグラフにしてみました。magnitudeは入力信号の大きさによって変化します。また、入力レベルを一定にしていてもある程度の範囲で大きくなったり小さくなったりします。各周波数におけるmagnitudeの値の最小値、最大値を読み取りました。

　OZ1JHM CW decoderのオリジナルスケッチで使用されている

　target周波数=558Hz、sample周波数=8928Hzの場合の周波数に対するmagnitudeの値のグラフです。確かにtarget周波数にピークがあり、帯域幅200Hz程度のフィルタになっています。8928/558=16で、入力信号１周期に対して16回のサンプリングが行われます。M(Low)が各周波数のmagnitudeの最小値。M(High)が最大値です。

target=625Hz、sample=10000Hzの場合のグラフです。上のグラフとほぼ同じ形になりました。

　

　target=800Hz、sample=10000Hzの場合です。

　maginitudeのピークが800Hzに移動しています。

　Goetzelアルゴリズムを使った周波数検出器をRaspberry Pi Pico　micropythonで同じことができるかどうか実験してみましょう。