モールス信号decoderの製作
暫く避暑を兼ねて田舎に行っていた。もう田舎は秋。いい加減な格好でいると、寒いぐらい だったよ。
帰りに土産で、ララベリージャムと言うのを買ってきた。桑の実のジャムだ。子供の頃、 学校の行き帰りに道草をして、桑の実を食べたなあと、懐かしくなったからだ。桑の 葉っぱはお蚕さまの食べ物、実は子供達のよきオヤツ。
親からは、腹を壊すから、食べちゃいけないと、きつく言い渡されていたものだけど、 よく食べた。家に帰ってきて母に、食べただろうと言われ、食べて無いと言い張っても 歴然とした証拠が残っていて、すぐにバレる。口の中や舌が真っ赤に染色されてしまうのだ。
子供ごころに、バレないように、うがいをしたり、奥歯でそっと咬んでそのまま飲み干す等を したけど、だめだった。そのうち、母は諦めて何も言わなくなった。 あれ？ 染色されちゃうやつは、桑の実じゃなくて、ヘビ苺だったかなあ。
肝心のジャムであるが、やたらに甘くて砂糖みたいだったよ。うつろな記憶では、少し 酸っぱくて、ワイルドな味だったような。。。
モールス信号のdecoder
あれから、随分と日が経過しちゃったけど、取りあえず版を書いてみた。 プログラムは、FreeBSD用で、cwd.c と lookup.c の２本１組になる。
まずは、cwd.c の方
// Morse code decoder
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/uio.h>
#include <unistd.h>
//#define SOURCE "/dev/dspW" // Morse tone from audio device
#define SOURCE "sound.au" // Morse tone from raw file
#define DATAN 128 // FFT point #
#define MAXELE 10 // max dot dash element (par char)
#define DUMP(fmt,val) do{ if (argc == 2 ){ printf(fmt, val); } } while(0)
double WaveL[DATAN]; // tW input wave
double SpecL[DATAN]; // fA output spectoram abs
double cmpXL[DATAN]; // fR output spectoram real part
double cmpYL[DATAN]; // fI output spectoram image part
void fft(double *tW, double *fA, double *fR, double *fI) {
int N = DATAN;
int N_2 = N/2;
int i,j,k,kp,m,h,P;
double mdd;
double w1, w2;
double t1,t2,s1,s2;
double tri[DATAN];
double fReal[DATAN], fImag[DATAN];
// Get P, where N = 2 ** P
i = N; P = 0;
while (i != 1) {
i = i / 2;
P++;
}
// copy real wave into buffer
for (i = 0; i < N; i++) {
fReal[i] = tW[i]; fImag[i] = 0.0;
}
// make sin,cos table
for ( i = 0; i < N_2; i++ ) {
tri[i] = cos( 2 * i * M_PI / N );
tri[i + N_2] = (-1.0) * sin( 2 * i * M_PI / N );
}
// sort
j = 0;
for ( i = 0; i <= N - 2; i++ ) {
if (i < j) {
t1 = fReal[j]; fReal[j] = fReal[i]; fReal[i] = t1;
t2 = fImag[j]; fImag[j] = fImag[i]; fImag[i] = t2;
}
k = N_2;
while (k <= j) {
j = j - k; k = k/2;
}
j = j + k;
}
// do butterfly
for ( i = 1; i <= P; i++ ) {
mdd = pow(2.0, (double)i);
m = (int)mdd;
h = m/2;
for ( j = 0; j < h; j++ ) {
w1 = tri[j*(N/m)];
w2 = tri[j*(N/m) + N_2];
for( k = j; k < N; k+=m ) {
kp = k + h;
s1 = fReal[kp] * w1 - fImag[kp] * w2;
s2 = fReal[kp] * w2 + fImag[kp] * w1;
t1 = fReal[k] + s1; fReal[kp] = fReal[k] - s1; fReal[k] = t1;
t2 = fImag[k] + s2; fImag[kp] = fImag[k] - s2; fImag[k] = t2;
}
}
}
// fft result copy for output
for ( i = 0; i < N; i++ ) {
fA[i] = sqrt(fReal[i] * fReal[i] + fImag[i] * fImag[i])/DATAN;
fR[i] = fReal[i];
fI[i] = fImag[i];
}
}
load_wave(double *tW, int fd){
int i;
signed short buf[DATAN];
read(fd, &buf, sizeof(signed short) * DATAN);
for (i = 0; i < DATAN; i++){
tW[i] = (double)buf[i];
}
}
main(int argc, char* argv[]){
int count = 491520 / DATAN; // 491520 = samples @ 8000/sec * 60sec
int fd;
double th = 800.0; // thereshold for 1/0
int dot = 5; // dot length
int onoff; // Signal on(1) or off(0)
int ms = 0; // mark(1) or space(0)
int mcnt = 0; // mark counter
int scnt = 0; // space counter
char val[MAXELE]; // dot(1),dash(3)
int x = 0; // index for val
int ws= 1; // word-space print enable
fd = open(SOURCE, O_RDONLY);
while(count--){
load_wave(WaveL, fd);
fft(WaveL, SpecL, cmpXL, cmpYL);
onoff = SpecL[10] > th ? 1 : 0; // SpecL[10] is amplitude of 625Hz
DUMP("%d",onoff);
if (onoff != ms ){ // mark->space or space->mark
ms = onoff;
if (ms) { mcnt = scnt = 0; ws = 1; } else { scnt = 0; }
}
if (ms) {
mcnt++;
} else {
scnt++;
if ( mcnt && (scnt > dot) ){ // in char
val[x++] = mcnt > dot * 2 ? '3' : '1';
mcnt = 0;
if (x == MAXELE){ exit(1); } // fatal, broken morce code.
}
if ( x && (scnt > dot * 2)){ // detect char
val[x] = '\0';
printf("%s", lookup(atol(val)));
DUMP("%s", "\n");
x = 0;
}
if ( ws && (scnt > dot * 5) ){ // detect word
printf("%s", "_\n");
ws = 0;
}
}
}
close(fd);
printf("%s", "\n");
exit(0);
}
次は、lookup.cだ。
// lookup morse code to char
#include <stdio.h>
struct kv {
long key;
char val[8];
};
struct kv tbl[] = {
{ 13, "A" }, { 3111, "B" }, { 3131, "C" },
{ 311, "D" }, { 1, "E" }, { 1131, "F" },
{ 331, "G" }, { 1111, "H" }, { 11, "I" },
{ 1333, "J" }, { 313, "K" }, { 1311, "L" },
{ 33, "M" }, { 31, "N" }, { 333, "O" },
{ 1331, "P" }, { 3313, "Q" }, { 131, "R" },
{ 111, "S" }, { 3, "T" }, { 113, "U" },
{ 1113, "V" }, { 133, "W" }, { 3113, "X" },
{ 3133, "Y" }, { 3311, "Z" },
{ 33333, "0" }, { 13333, "1" }, { 11333, "2" },
{ 11133, "3" }, { 11113, "4" }, { 11111, "5" },
{ 31111, "6" }, { 33111, "7" }, { 33311, "8" },
{ 33331, "9" },
{ 131313, "." }, { 331133, "," }, { 113311, "\?"},
{ 31113, "=" }, { 311113, "-" }, { 333111, ":" },
{ 133331, "'" }, { 31131, "/" }, { 133131, "@" },
{11111111, "<HH>"}, { 13131, "<AR>" },
{ 0, "<\?\?>" } // case not found. (don't delete this line)
};
char* lookup(long k){
int i = 0;
while(1){
if (tbl[i].key == 0){ return tbl[i].val; }
if (tbl[i].key == k){ return tbl[i].val; }
i++;
}
}
/******* for test ******************************************************
main(){
printf("%s\n", lookup(atol("3111"))); // "B"
printf("%s\n", lookup(atol("33331113"))); // "<??>"
}
*************************************************************************/
おまけで、Makefile
a.out: cwd.o lookup.o cc -lm lookup.o cwd.o -o a.out cwd.o: cwd.c cc -g -c cwd.c lookup.o: lookup.c cc -g -c lookup.c
設計方針
前回書いたグラフで、モールス信号のある所（以後、マークと呼称）、無い所（スペースと 呼称）の時間を計り、ドット、ダッシュの検出を行う。それを元に、文字の検出、ワードの 抽出を行う事にする。
上で、時間を計ると表現したが、どのくらいの時間分解能が望ましいのだろう？ 細かければ 細かいほど、良いに決まっているが、FFTとの兼ね合いもあるので、好きに設定する事は 出来ない。（細かさと、周波数分解能は、反比例する）
FFTをするためには、ある程度の測定データが揃わなければならない。通常は、データの個数を ２のべき乗にする。
たとえば、１２８個のデータを揃えるには、8000Hzでサンプリングしてるとすると、 128 * (1/8000) となり、これを計算すると、16msとなる。すなわち、16ms待っていると、 １回FFTが出来るデータが揃うのだ。
一方、FFTした時の周波数分解能は、サンプリング周波数/FFTポイント数となる。上の例を 当てはめると、8000/128 = 62.5Hzだ。 今回は、この設定で手を打つ事にした。
なお、8000サンプリング/秒は、FreeBSDの/dev/dspWのデフォルトに引きずられているだけで 他意は無い。（単なる、手抜きです）
上のうらはらの関係は、周波数分解能を上げようとすると、長時間に渡ってデータを採取する 必要があると言う、「ハイゼンベルグの不確定性原理」によるものだ
cwdのデフォルト値について
本当は、プログラム起動時にでも（理想的には、DSCWのように、信号に自動追従）変更 出来るようにすべきだが、面倒になるのでやめた。変更が必要ならソース上で修正して 再コンパイルしてください。
まずは、抽出周波数を固定とした。以前は440Hzを使っていたけど、マイク感度が悪かった ので、気持ち周波数を上げる事にした。600Hzぐらいを中心に検出して欲しいので、600/62.5 = 10 。 これが、 SpecL[10]; の根拠だ。
次は、ドット周期。100msに設定したかったので、100/16 = 6。なのだけど、調整の結果、 int dot = 5; とした。(A1A breakerの出力が正常にdecode出来るように調整したって のが、舞台裏な訳ですよ。）まあ、タイミングの関係で、プラス・マイナス１カウントは 避けられない誤差となりますから、細かい事は、気にしない。
マーク/スペースの判定信号レベル(double th = 800.0;)は、データを見ながら適当に決めた。 他意はありません。
プログラム上重要な定数、FFTのポイント数は、#define DATAN 128 で、決め打。信号を 何処から取り出すかも、#define SOURCE で決め打。今は、ファイルから取っているけど、 リアルにサンプリングしたかったら、/dev/dspW を選んでおいて、main内の while(count--) を、while(1) にすればOK。エラーが無い限り、回り続けるでしょう。
#define MAXELE 10 は、１文字を構成するドットとダッシュの最大個数。訂正符号の "HH" が、ドット８個で一番多いと思うんだけど（和文知らないから嘘かも？）、余裕を 見て、10にしておいた。（実際はnull文字が必要なので9個まで)それを超える場合は、 諦めて、プログラムを強制終了させるようにした。
プログラムの動き
load_wave()で、128個のデータを取ってくる。それをFFTして、該当周波数の信号を抽出。 その信号をスレッショルドで判定して、マーク/スペースを決定。
前のマーク/スペース状態から変化したかを判定、(変更)。マークかスペースどちらかの カウンターを増加。
スペースのカウンターが、ドット分を超えていたら、ドットかダッシュを決定して、配列に 記録。
スペースが、３ドット分なら（プログラム上は、２ドット分以上と言う、緩やかな設定にしてますが)、 １文字を構成するドット・ダッシュが揃った事になる。 ドットは、"1" ダッシュは、"3"として記録してるので、ツー・ト・ト・ト なら、"3111" という文字列になる。
この文字列(符号)を数字に変換してから、lookup を使って検索し文字に変換する。モールスのテーブルは、{3111,"B"} と言うような、構造体をずらずら並べた配列になっている。 検索したい数字と、構造体内の数字が一致したら、文字を返せば良い。
配列の最後には、番兵として数字の０が入っているので、そこまで行っても、一致しなかったら それは、不明な符号なので、 ＜？？＞ を返す。 後は、それを表示。
ワードを検出したら、スペースを表示するのが正当だろうけど、今回は、分かり易いように、 ("_"後)改行するようにしてみた。
DUMPと言うマクロは、debug用に仕込んだもので、信号のon/off状態を表示します。
結果発表
モールス練習ソフト A1A Breaker を、Tone 600Hz 60Char/分に設定して、初画面のテキストを キーイングさせ、それを１分間録音した。その生録データを、オフラインで解析しました。 まずは、実行時間です。(on VMWARE)
[sakae@nil ~/CW]$ time ./a.out > /dev/null real 0m0.388s user 0m0.046s sys 0m0.323s
0.388/60 = 0.6% っつう事で、まあこんなものかと思いますね。次は、実際の実行結果です。
[sakae@nil ~/CW]$ ./a.out RY_ TIRED_ OF_ SITTING_ BY_ HER_ SISTER_ ON_ THE_ BANK,_ AND_ OF_ HAVING_ NOTHING_ TO_ DO:_
途中からの録音なので、出だしの部分が単語の途中からになっちゃってますが、正しくdecodingさ れている事を確認出来ました。
次は、内部状態も含めて、表示してみます。(debug mode)
[sakae@nil ~/CW]$ ./a.out debug 0011111100000011111111111111111100000011111100000000000R 0000000111111111111111111000000011111100000011111111111111111100000011111111111111111100000000000Y 000000000000000_ 000000000000000011111111111111111100000000000T 000000011111100000011111100000000000I 0000000111111100000011111111111111111100000011111100000000000R 00000011111100000000000E 000000011111111111111111100000011111100000011111100000000000D 000000000000000_ 0000000000000000111111111111111111000000111111111111111111000000111111111111111111100000000000O 000000011111100000011111100000011111111111111111100000011111100000000000F 000000000000000_ 000000000000000011111100000011111100000011111100000000000S 0000000111111100000011111100000000000I 000000011111111111111111100000000000T 000000011111111111111111100000000000T 000000011111100000011111100000000000I 0000000111111111111111111000000111111100000000000N 000000011111111111111111100000011111111111111111100000011111100000000000G 000000000000000_ 000000000000000011111111111111111100000011111100000011111100000011111100000000000B 0000000011111111111111111100000011111100000011111111111111111100000011111111111111111100000000000Y 000000000000000_ 0000000000000000111111000000111111000000111111000000011111100000000000H 000000011111100000000000E 000000011111100000011111111111111111100000011111100000000000R 000000000000000_ .....
ノイズも無い、スピードも一定と言う限定条件ながら、ちゃんと動いているようです。良かった、 良かった。
嗚呼、今気づいたんだけど、このままじゃ、表示文字列がバファリングされちゃって、 リアルタイムに文字が表示されないな。後で何とかしよう。
(9/6 追記）リアルタイム解析をやってみた。心配してたオーバーヘッドも無く、無事に decodeしてくれたよ。上のソースのままだと、やはりバファリングされちゃって、単語が 確定された時(改行が入るので）、づらっと単語が表示された。これは通な人みたいで面白い んだけど、新人にはやはり、ちときつい。fflush(stdout); を、printfの次の行に入れて あげたら新人向けとなった。