Workbench(4)

ワオー６EQUJ5

昔々、有名な科学雑誌ネイチャーにある論文が掲載された。それによると、ある知性体が別の星に住む知性体にみずからの存在を知らせたい場合は、電波放射と言う手段を使うはずだと。

そしてその電波の周波数は、宇宙に最も豊富に存在する元素に関わるだろう。その元素とは水素。水素原子が1420MHzの電波を励磁する事はよく知られた事実だ。そして、少ない電力で遠方まで電波が届くようにするには、帯域を絞る必要がある。

この推論を信じて、とある大学の人が電波望遠鏡を使って観測を続けた所、エルビスプレスリーが没した時に、ついに上記の信号を捕らえた。それも、ただの一度だけ。

ちなみに、上記信号は、36進法によるSメーターの値だそうです。関係ないけど、おいらの持ってる FT-817のSメーターは、０からFまでの１６段階だったりします。

上記の信号を読み解くと、弱い信号がどーんと大きくなって、すぐにフェードアウトしてるって事です。当時、この1420MHzで帯域わずか10KHz内と言う信号は、地球人が発するはずは無く、すわ ETからの呼びかけと大騒ぎになったのでした。しかし、ただ一度だけってのが、ハムをやる人なら分かるけど、こんなに潔いETは居ないだろうと。。。何度も、信号が送られてくるはずだろうからまた、受信できるのではなかろうか？

NASAが協力を申し出た。国に予算を付けてもらって、観測を続けよう。でも、反対する仕分け人が居た。そんな無駄な事に血税は使えませんとね。日本にもこういう人が居て、科学者の総スカンを喰らってましたね。誰かは言わないけど。。最近、TOP500で一躍有名になったあれです。

民間が援助に乗り出した、HP社を作った、ヒューレッドとパッカードさん。この会社、ハムの好き者ならスペアナとかでネットワークアナライザーでお世話になってるはず。そして、マイクロソフトの共同創立者のポール・アレン。

NASAが達成できなかったSETIが、こうして始まった。

『ETよ、こちらからかけ直す。いつになるかはわからないけど。』

そして、イタリアの有名な物理学者、エンリコ・フェルミが、『みんな、どこにいるのか？』と問いている。

さあ、SETI@homeへ！

命名

前回、pipeに流れる信号をファイルに落とす為の手続き名を決めた。（ほとんど、閃きだけど）それと逆の事をやる手続き名は、決めかねていた。で、あれから、ずっと考えていて、やっと決めたよ。まあ、生まれてくる前と言うか、決めてから仕込むと言うのが、私の流儀なんですけど。

ファイルのデータを土管もとえパイプに放り込む訳だから、Lispから名前を拝借してきて、throw かなあ。でも、この名前って、ごみの不法投棄っぽくて、気にいらんな。こういう時は英語の辞書を引く。chargeで投入。dropで投入するってのが出てるな。でも、なんかピンとこないな。

名前一つで、血を血で洗うすさまじい抗争が続いているrubyを思い出しちゃったよ。（mapとcollect派の抗争が、来週練馬であるんかな。おいらもmap派と徒党を組たかったけど、訳有りで行けなくなった。まあ、最近はrubyなんて下敷きのごとく扱われていて、レールが敷いてあれば、細かい事はどうでもいいよって人がほとんどなのかな。やっぱり、一年に一度の祭りだから、丁々発止とやらなくちゃ、面白くないぞ）

ここは、もう素直にunix界隈から、dump/restore なんてのを貰ってくるか。そうすると、saveonは釣り合いが悪いんでこれは避けたいな。ああ、派閥の論理が炸裂する世界だ事。

ええい、こんな事に時間を費やすのは無駄。えいっと決めるよ。entry って単語から、最初の２文字を貰ってこよう。en、をIMEで変換してみると、円、得ん、獲ん、とか出てきたけど、雰囲気的には、縁かなあ。ファイルとパイプの縁を取り持ちますってね。

def saveon(a, fn='zmywb.dat'):                # data save
    with open(fn, 'w') as f:
        for n in range(SMPL):
            f.write(str(a[n]))
            f.write('\n')
    return array(a)

def en(fn='zmywb.dat'):                       # entry into pipe
    a = []
    with open(fn) as f:
        dats = f.readlines()
    jp = True if dats[0].find('j') > 0 else False
    for n in range(SMPL):
        try:
            v = complex(dats[n]) if jp else float(dats[n])
            a.append(v)
        except:
            v = 0.0+0.0j if jp else 0.0
            a.append(v)
    return array(a)

名前が決まれば（決めれば）、コードなんて、サクッと書けちゃうよ。saveonの方は、パイプの途中にも挿入出来るように、スルー端子を付けておいた。(tapのように、pipeの途中に挿入可）

enの方は、一応、パイプの先端に置く事を狙った。日本の伝統的なパイプ技法、ながし素麺も、そうめんが流れてこなければ始まりませんからね。データが多すぎたら、多い分は捨てるし、少なすぎたら、zeroを補填するようにした。

wb> vo | fft | saveon
wb> en | disp

MSDOSの偽物パイプより進化してて、時空を超えてパイプを接続出来ますから、一行に収まらない長いパイプも接続可能です。更に、データは、普通のファイルですから、editorで修正する事も可能。ちなみに、複素数データは

(-1.10800257858e-13+0j)
(-1.29358329604e-13+2.2226664953e-13j)
(573.44-4.98282838601e-13j)
(-1.95464477937e-12-1146.88j)
(860.16-8.83182416089e-13j)
   :

こんな風になってました。これで、胸の閊えも取れましたから、少し新しい事をやってみます。今度は変調の反対操作である復調です。やっぱり、直交性は大事って事で。。。

変復調

前回、周波数変調とかのコードを書いたけど、改めて調べてみたら、どうも嘘っぽいぞ。無線機の変調・復調についてなんてページを見ると、へぇーとなってしまう。まだまだ修行が足りんな。本格的にやるなら、ちゃんとした本を読破しないといけないね。どんなのが良いのだろうか？

はじめて学ぶディジタル・フィルタと高速フーリエ変換このあたりが適当かなあ？誰もが悩む、負の周波数なんて言うコラムも載っているし。。。フィルター設計プログラム（ソース付き）が公開されてたので、先走って頂いてきてしまった。いつもお世話になります、CQ出版さん。

包絡線検波

毎度おなじみの百科事典から復調を引いてみる。すると、まあ、いろんな復調方式がある事が分かる。

一番馴染みがあるのは、ダイオードを使った、包絡線検波だな。早速やってみる。ダイオードって、正負の信号のうち正の信号しか通さないやつだな。半端もとえ半波整流回路と全波整流回路があるんだ。

電源回路だと、この整流回路の後ろに平滑回路が有ったな。なんだ、電源の整流回路と一緒ぢゃん。ならば、全波整流回路の方が効率よかったな。全波整流って、負の信号部分も反転して通しちゃうんだから、デジタル的には、絶対値を取るのと一緒だ。絶対値を取るやつは、workbench内に既にあるぞ。そう、getaだ。

wb> vo | am(fc=100) | geta | fft | disp

整流した波形のFFTを取って表示してみた。確かに、低い周波数領域に、被変調波が現れているな。後は、これをローパスフィルターで取り出せばいいんだ。えっと、今出来るローパスフィルターは、前回バラックで組んだ、周波数領域でのローパスフィルターだな。ちゃんとしたモジュールにしておくか。

def lff(a, fc=20):                            # lf module
    v = lf(fft(a), fc)
    return array(getr(ifft(v)))

本当は一行で書けるけど、debugの為に分けておいた。こうしとけば、debug用のプローブも当て易いしね。

wb> vo | tap | am | geta | lff | disp

信号レベルがシフトしてるのと振幅が小さくなったのを除けば、ちゃんと復調出来る事が確認できましたよ。

で、次は。。。

プロダクト検波

いよいよプロダクト検波です。SSBやらDSB信号の復調に使われるやつね。プロダクトって言うぐらいだから、掛け算器が必要なのね。用意しておいて良かった！もう球の7360はディスコンだしね。

wb> vo | dsb(fc=50) | mul(st(fs1=50)) | fft | disp

スケール内にスペクトラムが収まるよう、搬送波周波数を下げてみました。mulの所が検波器、変調した時と同じ周波数を注入しています。

スペクトラムを見ると、低い周波数領域と、搬送波の２倍のエリアにスペクトラムが表れました。ちゃんと、三角数学通りになりましたよ。

フィルターを通して高い周波数をろ過してみました。ちゃんと元の波形（の半分の振幅）が再現出来ましたよ。

ちょっと悪戯して、復調用のキャリアー周波数をずらしてみました。あのモガモガと言うか、モコモコ（どこかの議員か？）というか、復調後の波形は、大いに崩れていましたよ。

CWの再生

で、CW信号の再生はどうすんだっけ？ BFOでビート取って聞くのが一般的なんで、プロダクト検波でいいのかな？

wb> mul(st(fs1=300),dd(m=1,s=0)) | tap | mul(st(fs1=300)) | lff(fc=100) | disp

パイプの初段はCW送信機、tap は、一応送信波形のモニターです。実世界では、アンテナを経由して電波発射。地表と電離層の間を往復しながら電波は旅を続けると。受信部ではダイレクトに検波ですよ。その後、フィルターを通して可聴信号を取り出し、音を鳴らす。（今は、イヤホンの代わりにオシロスコープを繋いでいますけど）

なお、ローパスフィルターの帯域を広くしているのは、ドット（の高い周波数成分）を正しく再生するためです。

つう事で、上のパイプは、トランシーバーを模したものになっています。これって、ひょっとしたらエレクラフトのKX1 ぐらいのプログラム表現になるんでしょうかね？

証明

三角数学の有名な定理に、加法定理と言うのがある。大事な定理なので受験生諸君は家宝にしとくと、きっと果報が訪れる事じゃろう。で、この加法定理から導かれるものに、積和公式と言うのがある。

            cos(A-B) + cos(A+B)
cosAcosB = --------------------
                     2

cosAとcosBを掛けると、cosのA-B とcosのA+Bが出てくると言うんだ。もし、Aが３でBが５なら、３引く５で、－２と、３足す５で８の周波数が出てくるはず。あれれ？マイナスの周波数だよ。本当に出てくるんか？

wb> mul(st(fs1=3), st(fs1=5)) | fft | geta | df

２と８が出てきた。掛け算は、A掛けるB とやっても、B掛けるAとやっても同じになるんで、上の場合は Aが５でBが３と置いても、何ら問題は無い訳だ。

すなわち、周波数領域で考えたら、周波数０の所で、マイナスの周波数は、折り返してるって思ってよさそうって事だな。こんなんで、積和公式の証明になってるんですかね？単に実験しただけじゃん！

積和公式には、cosとsinの掛け算はどうなるとか、別の式が有るけど、家にある唯一の数学の本では、そんなのは忘れてしまっても良い事になってます。三角で一番重要なのは、sinの二乗とcosの二乗を加えると１になるよと言う式。これは３平方の定理（ピタゴラスの定理）の言い換えに過ぎないと言う事情があるからです。

数学の人って、基本の所からどんどん発展させて定理を作って行くなあ。想像力につくづく脱帽しますよ。

で、脱帽ついでにつらつらと本を繰って行くと、デルタ関数なんてのが出てくる。xがzeroの時、デルタ関数の値は無限大。xがzero以外の時はデルタ関数の値はzero。無限な区間に渡ってデルタ関数を積分すると、結果は１になると言うんを発明しちゃったらしい。よく、こんなの想像するよ。もうついて行けんです。

でも、このデルタ関数、フーリエ変換と深ーーーい関係があるそうですから、以外だなあ。ああ、インパルスとお友達かも。こんなのに付き合ってたらどうしようも無くなってしまいそうなので、もう少し現実的なやつを。

この記事の冒頭で出てきた、Sメータの振れ具合のデータ。６EQUJ5 こんな説明が出てましたよ。

波形にして眺めてみましょうかね。ひょっとしたら、ETからのメッセージが潜んでいないかな。

まず、Sメータの振れを、10進数に直してファイルに落としておきますか。ファイル名は、s.log ぐらいでいいかな。

早速、分析してみます。

wb> en(fn='s.log') | disp

wb> en(fn='s.log') | fft | disp

時間領域では、何て事ないけど、周波数領域では綺麗なやつが表れました。はて、ETはこれに何を託したのだろうか？謎は謎を生みますよ。