JH7UBCブログ

アマチュア無線 電子工作 家庭菜園など趣味のブログです

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

2021-12-27 20:49:13 | Arduino
 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で、入力信号1周期に対して16回のサンプリングが行われます。M(Low)が各周波数のmagnitudeの最小値。M(High)が最大値です。

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

 target=800Hz、sample=10000Hzの場合です。
 maginitudeのピークが800Hzに移動しています。


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

Raspberry Pi Pico MicroPython RTC I2C LCD2004表示

2021-12-20 18:34:07 | Raspberry Pi Pico
 Raspberry Pi Pico MicroPythonプログラミングで、RTC(リアルタイムクロック)をI2C LCD2004に表示してみます。

 LCD2004の表示は、今回はライブラリを使います。I2C LCD1602/2004用のライブラリは、ネットやYouTubeで調べると複数あるようです。
 その中で、T-622 /RPI-PICO-I2C-LCDというライブラリを使ってみます。GitHubのこちらのページにアクセスします。 
 lcd_api.py とpico_i2c_lcd.py を開いて、Raspberry Pi Picoに同じファイル名でコピペします。
 回路図です。LCD1602/2004は5V動作ですが、直接I2Cバスに接続することもできるのですが、安全のため、3.3V/5VのレベルコンバータモジュールPCA9306を入れています。(秋月電子のモジュールで、3.3V側、5V側両方にプルアップ抵抗が入っています)I2Cは、id=1,sda=Pin(14),scl=Pin(15)を使います。

 スクリプトです。main.pyとして保存します。1secごとにTimer割込みを発生させて、RTCの表示を更新しています。Pico内蔵LEDは、1secごとに点滅します。Thonnyを経由して、パソコンの時計の値が書き込まれますので、Picoをパソコンに接続している場合は、現在の時刻が表示されます。
 Thonnyを使わないで、電源だけを入れると2021/01/01 00:00:00と表示され、カウントアップしていきます。
-----------------------------------------------------------------------------------------------------
from machine import I2C,Pin,Timer
import utime
from lcd_api import LcdApi
from pico_i2c_lcd import I2cLcd

I2C_ADDR      = 0x27
I2C_NUM_ROWS = 4
I2C_NUM_COLS = 20

i2c = I2C(1, sda=Pin(14), scl=Pin(15), freq=400000)
lcd = I2cLcd(i2c, I2C_ADDR, I2C_NUM_ROWS, I2C_NUM_COLS)
timer=Timer()
led = Pin(25,Pin.OUT)

lcd.move_to(3,0)
lcd.putstr("Real Time Clock")

def clock(t):
     global led
     led.toggle()
     lcd.move_to(0,2)
     time=utime.localtime()
     lcd.putstr("{year:>04d}/{month:>02d}/{day:>02d} {HH:>02d}:{MM:>02d}:{SS:>02d}".format(
             year=time[0], month=time[1], day=time[2],
             HH=time[3], MM=time[4], SS=time[5]))

timer.init(freq=1,mode=Timer.PERIODIC,callback=clock)

while True:
     pass
-----------------------------------------------------------------------------------------------------
ブレッドボードです。



Raspberry Pi Pico MicroPython レシプロカル(逆数)カウンタ+CWトーンインジケータ

2021-12-19 12:37:02 | Raspberry Pi Pico
 Raspberry Pi Pico MicroPythonの組み合わせで、レシプロカル(逆数)カウンタを作ってみます。
 信号の周期T(sec)を測定して、f(Hz)=1/T(sec)の式で計算して、周波数を求めます。
 回路図です。周波数カウンタと同じです。



 スクリプトです。以前にmicro:bitでレシプロカル・カウンタを作った時のものをPico用に移植しました。CWトーン・インジケータとして使えるように、700Hz±50HzでLEDが点灯するようにプログラミングしました。サンプリング間隔を短くする場合は、下から2行目のutime.sleep(0.5)の値0.5を0.2などに変更します。
-------------------------------------------------------------------------------------------------
from machine import Pin,I2C
import utime
import ssd1306

i2c=I2C(0,sda=Pin(16),scl=Pin(17),freq=400000)
oled=ssd1306.SSD1306_I2C(128,64,i2c)

sig=Pin(15,Pin.IN)
led=Pin(25,Pin.OUT)

oled.text("Reciorocal",16,0)
oled.text("Counter",56,8)
oled.text("Hz",88,32)
oled.show()

old_value = 0
flag = 0

def display():
     global freq
     oled.fill_rect(48,32,40,8,0)
     oled.show()
     oled.text(str(freq),48,32)
     oled.show()
     if freq>=650 and freq<=750:
         led.value(1)
     else:
         led.value(0)

#main loop
while True:
     now_value =sig.value()
     if old_value == 1 and now_value == 0:
         if flag == 0:
             start_time = utime.ticks_us()
             flag = 1
         else:
             end_time = utime.ticks_us()
             freq = int(1000000 / (end_time-start_time)+0.5)
             display()
             flag = 0
             utime.sleep(0.5)
     old_value = now_value

-------------------------------------------------------------------------------------------------
 ブレッドボードです。
 自作のSGから1000Hzを出力して、測定しています。



 ±50Hzくらいの誤差があり、測定値が細かく変化します。
 utime.ticks_us()の値のふらつきが原因でしょうか。

Raspberry Pi Pico MicroPython HC-SR04 距離の測定

2021-12-13 12:35:36 | Raspberry Pi Pico
 Raspberry Pi Picoで、超音波センサーHC-SR04を使って距離を測定するプログラムをMicroPythonで組んでみます。

 MicroPythonでのプログラムは、既にmicro:bitで実験済みです。記事はこちら。測定原理はその記事を参照してください。なお、今回は、Pico内蔵の温度センサーから温度を読み取り、音速を計算します。

回路図です。表示はOLEDを使い、I2Cは、id=0,sda=Pin(16),scl=Pin(17)とします。HC-SR04のTRIGはGPIO14に、ECHOはGPIO15に接続します。HC-SR04からのECHO出力は、5Vなので、抵抗で分圧してPicoに入力します。



 TRIGに加えるパルスは、10us以上となっていますので、
 trigger=Pin(14,Pin.OUT)  とし、
 trigger.value(1)
 utime.sleep_us(2)
 rtigger.value(0)
 で、triggerパルスをCH-SR04に出力します。この時のパルスの波形です。

 約20us幅のパルスが発生しました。(パルス幅は若干広くなったり狭くなったり、ゆらぎます。)utime.sleep_us()を省略するとパルス幅は、10us以下になってしまいます。

 HC-SR04から帰ってくるECHOパルスを測定してみました。


 echoパルスの立ち上がりの時刻t1と立下りの時刻t2をutime.ticks_us()で、us単位で測定し、t=t2-t1からechoパルスの幅(us)を計算することができます。

 tは、HC-SR04から発せられた超音波が戻ってくるまでの時間で、距離の2倍を音が進む時間です。従って、音速をVs(m/s)、HC-SR04と物体の距離をL(m)とするとL(m)=Vs(m/s)×t(s)/2で計算できます。
距離をL(cm)、t(us)とすると、L(cm)=(Vs(m/s)×100×t(us)/1000000)/2=Vs(m/s)×t(us)/20000で計算することができます。なお、音速は温度によって変わり、Vs(m/s)=331.5+0.6×temp(℃)で計算できます。

 スクリプトです。
-----------------------------------------------------------------------------------------------
from machine import ADC,Pin,I2C
import utime
import ssd1306

trigger=Pin(14,Pin.OUT)
echo=Pin(15,Pin.IN)

i2c=I2C(0,sda=Pin(16),scl=Pin(17),freq=400000)
oled=ssd1306.SSD1306_I2C(128,64,i2c)
oled.text("RP2 HC-SR04 TEST",0,8)
oled.text("TEMP:      C",18,24)
oled.text("L(cm):",18,40)
oled.show()

trigger.value(0)#Triggerの初期値0

coeff=3.3/65535
a=ADC(4)

def mesurement():
     global temp,L
     #Trigger pulse >10us
     trigger.value(1)
     utime.sleep_us(2)
     trigger.value(0)
    #反射時間の測定
     while echo.value()==0:
         t1=utime.ticks_us()
     while echo.value()==1:
         t2=utime.ticks_us()
     t=t2-t1

    #気温の測定(pico内部の温度計を利用)
     v=a.read_u16()*coeff
     temp=round((34-(v-0.706)/0.001721),1)
    #距離の計算
     vs=331.5+0.6*temp#音速の計算
     L=round((vs*t/20000),1)#L(cm)

def display():
     global temp,L
     oled.fill_rect(72,24,32,8,0)
     oled.show()
     oled.text(str(temp),72,24)
     oled.show()
     if L < 2000:
         oled.fill_rect(72,40,56,8,0)
         oled.show()
         oled.text(str(L),72,40,)
         oled.show()

while True:
     mesurement()
     display()
     utime.sleep(1)

-----------------------------------------------------------------------------------------------
 Picoの内蔵温度計の値はADC(4)の電圧として読み取れます。その値を温度に変換する式は
 27-(v-0.706)/0.001721 ですが、実際より低い値になってしまいますので、便宜的に27を34に変更しました。
 1secごとに距離を測定して、OLEDに表示します。

 ブレッドボードです。


 温度の値は1℃以内で細かく変動します。

 実際に物体との距離を測定してみました。距離の表示も1cm以内で細かく変動します。

測定値です。
距離(cm) 測定値(cm)
 10    10.3
 20    20.6
 30    30.1
 40    39.9
 50    49.9
 60    59.4
 70    69.5
 80    79.0
 90    88.8
100    98.9

 micro:bitの時より、かなり精度が良くなっています。1m以内なら±5mm程度の誤差で測定ができそうです。


Raspberry Pi Pico MicroPython 周波数カウンタ

2021-12-10 20:23:14 | Raspberry Pi Pico
 Raspberry Pico Pi をMicroPythonでプログラミングして、周波数カウンタを作ってみます。

 回路図です。入力信号を簡単なアンプで増幅して、GPIO15に加えます。
 入力波形の立ち上がり(0→1)をカウントします。Timer割込みで、1secごとに1sec間のカウント数をOLEDに表示します。


ブレッドボードです。入力ゲートが開いているとき、内蔵のLED(Pin25に接続されている)が点灯します。


 スクリプトです。1sec間ゲートを開け(gate_flag=1)、信号の立ち上がりをでカウントします。1sec経ったら、ゲートを閉じ(gate_flag=0)、count数をOLEDに表示します。ゲートのコントロールはTimer割込みで、1secごとに割込みを発生させています。
 OLEDを使いますので、ssd1306ライブラリをインストールしておき、このスクリプトは、main.pyという名前でPicoに保存します。
 数字を右寄せで表示したいのですが、MicroPythonでは、右揃え関数rjust()が使えないので、文字数を数え、文字数が足りないときは、左側にスペースを入れるプログラムにしました。
-----------------------------------------------------------------------------------------------------
from machine import Timer,Pin,I2C
import utime
import ssd1306

led = Pin(25, Pin.OUT)
sig = Pin(15,Pin.IN)

i2c=I2C(0,sda=Pin(16),scl=Pin(17),freq=400000)
oled=ssd1306.SSD1306_I2C(128,64,i2c)
oled.text("RP2 Fre Counter",0,8)
oled.text("Hz",96,24)
oled.show()
led.value(0)

#Timer生成
timer = Timer()
#counter関係設定
counter = 0
old_sig = 0
gate_flag=0

def gate_ctrl(timer):
     global led,counter,gate_flag
     gate_flag = 1 - gate_flag
     if gate_flag == 0:
         oled.fill_rect(48,24,40,8,0)
         oled.show()
         count=str(counter)
         if len(count)<5:
             for i in range(5-len(count)):
                count=" " + count
         oled.text(count,48,24)
         oled.show()
         counter=0
     led.value(gate_flag)

timer.init(freq=1, mode=Timer.PERIODIC, callback=gate_ctrl)

while True:
     if gate_flag == 1:
         sig_val=sig.value()
         if sig_val==1 and old_sig==0:
             counter += 1
         old_sig = sig_val
-----------------------------------------------------------------------------------------------------
  さて、どの程度の周波数まで測定できるでしょうか。自作のSGから信号を加え、少しずつ周波数を上げていきます。



10KHzまでは、ほぼ正確な値をカウントします。MicroPythonでのプログラミングでは、この辺が限界のようですが、オーディオ周波数なら実用になりそうです。