Goertzel algorithm
酒飲みの為の、酒呑みによる、 「日本酒カレンダー」 を教えてもらったので謹んでメモしておく。
probe goertzel
前回やったCW信号のデコード。この核心は、goertzelというアルゴリズムだった。 単一周波数が存在するか否かを効率よく検出できるとな。これでもう、FFTとはおさらばできる 。
内部でどんな計算がおこなわれているか、プローブもとえgdbをあててみたくなった。そこで困るのが、信号がon/offしてる事。(CWなんであたりまえ)
そんな訳で、連続波が欲くなったんだ。確かsoxでも発生できたはず。sox 正弦波 発生 とか 叫んでみたよ。
$ sox -r 10000 -b 16 -c 1 -n sweep.wav synth 10 sine 100-3000 gain -10
こんな風にすると、100から3000Hzまで10秒かけて掃引しながら正弦波を発生してくれるとな。なんだか、これってGNU RadioにVCOって名前で登録されてたな。
gain A power gain in dB. Zero gives no gain; less than zero gives an
attenuation.
gainの-10って何かなと調べてみた。ここには説明ないけど、絶対基準は、CDのフルスケールの音量みたい。大体こんな大音量だとクリッピングしちゃうので、-3dB (~70%)が実用的らしい。
ああ、検索結果に、こんな楽しいコンテンツも出てきた。
作者さんは、オイラーと同郷のsakura村の方。オイラーよりは富豪な方で、FreeBSDをバリバリ使っておられた。
こうなると、同じような事を、やってみたくなるのが人情というものです。
あれ？ Goertzelはどうした？
とりあえず、ワシントン大学の資料をみましょ。
work bench
どうも、上の掃引波形と言い、gnuplotによる表示機構といい、ネットワークアナライザーを作って実習しろと言う圧力があるな。
#include <stdio.h>
#include <math.h>
float sampling_freq=10000.0;
float target_freq= 600.0;
#define N 100 // BW = sampling_freq / N
short int Data[N];
float omega,coeff,mag2;
float Q0,Q1,Q2,sine,cosine;
int rc,i,k;
int main (){
k = (int) (0.5 + ((N * target_freq) / sampling_freq));
omega = (2.0 * M_PI * k) / N;
cosine = cos(omega);
coeff = 2.0 * cosine;
for (;;){
rc = fread(&Data, sizeof(short int), N, stdin);
if (rc < 1)
break;
Q1 = 0; Q2 = 0;
for (i = 0; i < N; i++){
Q0 = coeff * Q1 - Q2 + (float) Data[i];
Q2 = Q1; Q1 = Q0;
}
mag2 = (Q1*Q1) + (Q2*Q2) - Q1*Q2*coeff;
printf("%f\n", sqrt(mag2));
}
}
cwdecoder.cから核心部分を抜き出してきました。これ、中心周波数600Hz,帯域幅100Hzのバンドパスフィルターとして機能するのね。FFTみたいに、データ個数を2のべき乗にする必要なしってのが、便利でよさそう。
$ cat pg.out #! /usr/local/bin/gnuplot #set terminal png set nokey set grid #set logscale y plot '< cat -' with line pause 10
こちらは、gnuplot script です。pauseで10秒待つようにしました。こうしておかないと、 表示してすぐにグラフが消えてしまうからです。10の代りに -1 とすると改行入力がされるまで、停止するはずなんですが、機能しませんでした。まあ、スクリプトにチョッカイができない のは、あたり前の事だわな。
下記のようにパイプの終端で使っているので、セッションリーダーはcatが担当。catが終了すると、パイプの後段も全て死ぬので、悠長にキー入力なんて待っていられない？ そう思って、単独で動かしてみたけど、やっぱり無視されてましたよ。(-pなんてオプションも機能しなかった)
$ cat sweep.raw | ./dut | ./pg.out
dutてのは、上のプログラムね。信号源から被試験デバイスのフィルターに信号を入力。その応答を可視化するために、オシロスコープを使う。これってネットワークアナライザーの原型ですよ。
現役時代は、そのオシロにポラロイドカメラを取り付けて、写真にとる。オシロの輝度やら露光を間違えると、白けた写真になっちゃう。今なら、pg.outスクリプトのset terminal pngを有効にするだけで、綺麗なpngが得られる。( .. | ./pg.out > picture.png 等とする)
signal generator
soxに作らせた掃引音源が、どうも周波数をリニアに発生していないようだ。低い周波数の時は長時間、高い周波数になるにつれて、短時間になってるっぽい。
予想するに、発生波形の個数を一定になるようにしてないかな。そうすると、(ある)周波数の持続時間は、周波数に反比例するな。正確に調べるならsoxのソース嫁。
面倒なんで、信号発生器を自作しちゃえ。
ob$ cat sg.c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define SAMPLE_FREQ (10 * 1000)
#define LEN 100
void gen(int freq){
double out;
int t;
short int b16;
for (t=0; t<LEN; t++){
out = sin( 2 * M_PI * freq * t / SAMPLE_FREQ );
b16 = (short int) (out * 3276);
fwrite(&b16, sizeof(short int), 1, stdout);
}
}
int main( int argc, char *argv[] ){
int f;
//a = atof(argv[1]);
for (f=100; f <= 3000; f += 10){
gen(f);
}
}
arduinoのA/Dの変換スピードに敬意を表して、サンプリング周波数は、10KHzにした。なんてのは、表向きの理由で、単に暗算が楽だからです。
サンプリングの個数と言うか長さは100個。時間にすると、10ミリ秒だけの期間で各周波数を発生。
発生させる周波数は、100Hzから3000Hzまで、10Hzステップという、プリセット版のスィープジェネレーターです。信号レベルは、ほぼ -10dB です。( 3276 / 32768 )
DUTに直結できる様に、信号出力は、バイナリーなrow形式にしています。
goertzel の Q
ob$ ./sg | ./dut | ./pg.out > scan.png ob$ gpicview scan.png
ピークが500Hzにあって、あれって思わないでください。X軸の最左が100Hzになってますんで、あしからず。
このフィルターのQはどれぐらいになるんですかね？ センター周波数の信号レベルから3dB落ちた所の帯域幅との比って定義でしたっけ？
Freq BW=100 BW=50 BW=25
------ -------- -------- --------
41 0 0 0
42 1769 1639 1203
43 3667 2724 878
44 5536 3047 1041
45 7345 2339 1794
46 9151 0 0
47 10967 3585 2783
48 12668 7460 2336
49 14021 11162 3546
50 14823 14025 11257
51 14999 14999 14999 <== Center
52 14633 13956 11403
53 13880 11453 3512
54 12853 7752 2301
55 11550 3520 2929
56 9910 0 0
57 7930 2264 1940
58 5732 3338 1009
59 3540 3017 840
60 1585 1576 1350
61 0 0 0
BWを変えて(N=100,200,400)中心周波数+-100Hzの範囲(10Hz step)のレスポンスを表にしてみた。
サンプル数を変えると、値も変わるので、Nで除して正規化してある。また、入力振幅によって値が変わるのは常識なんで(線形応答)、見易い数値になるように、sgの振幅係数を、30000に設定しておいた。
sox synth
冒頭の所であげたsoxのsynth機能での掃引周波数発生。低い周波数はゆったりと、高い周波数はせわしなくってのが、どんな具合になってるか確認。どうやって？
無線家の常識、F = 1 / T という定義を思いだす。周期Tの逆数が、周波数Fだよってあれね。 周期なら、ちょっとの工夫で測定できるだろう。
#include <stdio.h>
#define SAMPLE_FREQ (10 * 1000)
#define N (SAMPLE_FREQ / 10)
short int Data[N];
int m2p(int offset){
int t;
for(t = offset; t<N; t++)
if (Data[t] < 0) break;
for(t; t<N; t++)
if (Data[t] > 0) break;
return (t);
}
void decide(){
int t,t0,t1;
t0 = m2p(0);
t1 = m2p(t0 + 1);
t = t1 - t0;
if (t)
printf("%d\n", SAMPLE_FREQ / t);
else
printf("0\n");
}
int main (){
int rc;
for (;;){
rc = fread(&Data, sizeof(short int), N, stdin);
if (rc < 1)
break;
decide();
}
}
0.1秒ごとにデータの塊を切り出す。その塊(配列)の中をスキャンして、マイナスからプラスに転じるポイントを2点みつける。その差が、周期になる。ポイント間は、1/SAMPLE_FREQという時間になってるからね。
すなわち、0.1秒間隔で、その時の周波数を求めているんだ。X軸は、発生開始からの経過時間(相当)、Y軸は周波数だ。Y軸を対数表記にしたら、直線になったよ。すなわち、経過時間に対して、周波数を対数的に発生させてるんですな。
音を扱う世界の対面には、人間がいる。物理的な刺激量と知覚量は、対数関係にあるっていう、ウェーバーヘィヒナの法則ってのが良く知られている。周波数検知も同様なの？
世の中には、絶対音感と言う特技の持ち主がいるそうだ。そういう人って、人間の格好をした周波数カウンターと思うんだけど、どうよ。
sin wave gen in golang
package main
import (
"encoding/binary"
"math"
"os"
)
const PI = 3.1415926
func main() {
var (
data float64
frequency float64 = 440 * math.Pow(2, 3/12) // 523 Hz
rate float64 = 44100
duration float64 = 3 // 3 secs
)
for n := 0; n < int(rate*duration); n++ {
// Calculate amplitude at the n / rate phase.
radius := frequency * 2 * PI * float64(n) / rate
data = math.Sin(radius)
// Render audio data to STDOUT.
binary.Write(os.Stdout, binary.LittleEndian, data)
}
return
}
go run main.go | play -r 44100 -c 1 -t f64 -