workbench(2)
空のおくには何がある。 空のおくには星がある。 星のおくには何がある。 星のおくには星がある。 眼には見えない星がある。 みえない星はなんの星。 （金子みすず『みえない星』より)
みすずさんは、ACでいい加減うんざりだったんだけど、石黒正人さんが『ALMA電波望遠鏡』（ちくまプリマー新書）の 最後で、深遠な天の川の写真と共に載せておられたので、引いてみた。
この書、電波少年だった石黒さんが趣味が高じて、電波天文学に足を踏み入れ、その意義と、世界最高の 電波望遠鏡の話が出てくる。
星から出てくる電波を観測する事により、星の誕生の過程や状態が分かるそうな。電波と言ってもいろいろな 波長があり、観測出来る状態も変わるとか。ある波長の電波は空気やらで吸収されてしまうため、なるべく 空気の薄い所を求めた所、チリの高地の砂漠が選ばれたとか。
望遠鏡の解像度を上げる為には、パラボラアンテナの口径を大きくしなければならないけど、１基で巨大なのは は限度がある。そこで、複数のパラボラを組み合わせ、それぞれのパラボラからの出力を計算処理（相関器） する事により性能を上げているとか。
勿論、パラボラのヘッド部に付いているプリアンプは、熱雑音を極力小さくする為、絶対零度近くまで 冷やすとか。地上からの人工雑音から逃れるため、人里離れたチリの山奥ってのだから、大変だな。 人工雑音が無い理想の場所は、月の裏側ってのには、至極納得しました。
この望遠鏡のフルオープンは、２０１２年との事だけど、一部では観測が始まっていて測定例が出ていた。 横軸が電波の波長（周波数）、縦軸が信号強度って事で、これFFTの画面じゃないですか。見方によっては、 ガスクロマトグラフとも見えるな。
いずれにしても、最先端をやると全部自分達で設計して自作する必要があるってのが、アマチュア（無線）に 通じていて親近感が沸くよ。それにしても、お金が大層かかるので、アメリカ、ヨーロッパ、アジア（日本、 台湾）の共同運用ってのには驚かされる。
電波少年（JA2***)さんの長年に渡る情熱に乾杯。
my REPL ?
先週やってたのの続きで、wireを定義するのに文字列である証をpythonに示してあげるの面倒なんて 事に気がついた。だったら、自前のREPLを書いちゃえって発想が出てくるな。今あるrun()をちょいと 手直ししてあげればいいんか。readは、raw_inputすればいいんだな。printは、まあ 主目的が、グラフへ出すんで、いらんか。
本格的な例は、 (How to Write a (Lisp) Interpreter (in Python))を読めばいいんだな。
余談はさて置き、run()をどうするかだな。どうせならrun()を単独でも使いたいぞ。そうすれば、ipythonとか のヒストリー機能がそのまま使えるからね。ついでに、例外の処理も何とかしておこう。run()の中で、最後に 評価した値を返すようにしとけば、後で役に立つかも知れんな。(Schemeのbegin相当かな）
そんな欲張り精神で手直ししたのが、以下の部分だ。
def run(wire):
try:
exlude = set(['show()',])
cmds = wire.split('|')
for n, cmd in enumerate(cmds):
cmd = cmd.strip()
if (n == 0) or cmd in exlude:
do_cmd = (cmd if cmd[-1] == ')' else (cmd + '()' ))
else:
do_cmd = (cmd.replace('(', '(val,', 1) if cmd[-1] == ')'
else (cmd + '(val)' ))
val = eval(do_cmd)
return val
except Exception as e:
print("Error: {0}".format(e))
def wb(): #################### my REPL
while True:
run(raw_input('wb> '))
if __name__ == '__main__': wb()
普通にこのスクリプトを起動すると my REPL である、wb()が回り出すので、普通に使えば良い。 ここから抜ける方法は特に用意してないけど、Ctrl-dすれば、文句を言われつつもshellに復帰する。
wb> tt|dt
wb> ；；Ctrl-dしたよ
Traceback (most recent call last):
File "C:\homes\WORK\wb.py", line 106, in <module>
if __name__ == '__main__': wb()
File "C:\homes\WORK\wb.py", line 104, in wb
run(raw_input('wb> '))
EOFError: EOF when reading a line
>>>
>>> run('tt|fft|geta|df')
>>> x = run('tt|fft')
>>> x
array([ -1.28605937e-13 +0.00000000e+00j,
6.37438194e-14 +2.51167817e-13j,
-2.29666762e-13 +5.70902738e-14j, ...,
1.63840000e+03 -7.59037277e-12j,
-2.29666762e-13 -5.70902738e-14j, -8.03186872e-13 +6.12763016e-13j])
>>> wb()
wb>
もう、いちいち wire = なんて書かなくてもいい。どうしても実体配線図を残したかったら、分かり易い 変数名に定義しておいて、それを、run()の引数として渡せばOK。run()の返り値は、最後に評価したものに なってるので、使いたかったらご自由に。実験を続けるなら、wb()で、 my REPL へ落ちていけば良い。 これ、実験用のブレッドボードだな。
wb> dd|fft|fcc Error: name 'fcc' is not defined wb>
ちゃんとエラーチェックもやってて、ボードが壊れないように保護してくれてるな。
それらしく
周波数ドメインの表示が今までは、折れ線グラフだった。まあ、これでもいいんだけど、拡大して いった時に、スペクトルが三角になってしまって、ちょっとかっこ悪い。pylabで棒グラフってどう すればいいんだろう？ あの膨大な資料に当たるのは勘弁な。
ぐぐったら Matplotlib サンプルコードを公開されてる ぐうたらさん(とは、とても思えない)という方がおられた。有り難く拝借させて貰った。また、 pythonのサンプルコードも、実用的で 嬉しいですよ。
def df(a): # display freq domain
#plot(a[:SMPL/2])
#bar(arange(SMPL/2), a[:SMPL/2])
bar(arange(SMPL/8), a[:SMPL/8])
show()
棒グラフを本当に一本々書いてくれるので、マウスの反応についてこれませんでした。よって、良く使う 周波数エリアに限って表示させてます。FFTすると、サンプリングサイズの丁度半分の所で、折り返し になるんで、演算量的には有利なんですけどね。
ドットとダッシュ
今年はモールスさん生誕220周年とかだそうなので、記念にモールス符号の発生でも思った。けど、テーブルを 用意するのが面倒なので、ドットとダッシュで表す事にした。デコードはご自由に。
def dd(msg='2133321331213116', m=0.5, s=-0.5): # dot dash (msg=JPL)
wave = SMPL * [s]
dots = 0
for n in msg:
if n == '1' or n == '3':
dots += (1 + int(n))
elif n == '2' or n == '6':
dots += int(n)
else:
pass
t = int(SMPL / dots)
ref = [m] * 3 * t
st = 0
for n in msg:
if n == '1' or n == '3':
wave[st:st+t*int(n)] = ref[0:t*int(n)]
st += (t * (int(n) + 1))
elif n == '2' or n == '6':
st += (t * int(n))
else:
pass
return array(wave)
ドットを文字の１で、ダッシュを３で表す事にします。文字区切りは２で、単語区切りは６です。 あれ、文字区切りは３だよねと言うなかれ、これは仕様です。（と、きっぱり。ちょいとした手抜きです）
信号レベルは、マークが0.5で、スペースが-0.5をデフォルトとしました。あれ、このレベルってRS232C でも無いし、ECLレベルでも無いし、TTLレベルとも違うなと思った人、鋭い。普通は、m=1,s=0 ぐらいに しとくんでしょうけど、オシロスコープで波形を眺めた時、オートスケールのおかげで、波形が目一杯に 拡大されてしまい、縦線しか見えなかったので、こうしてあります。
このオートスケールも、人間の感性に訴える設定って難しいよね。m=1,s=0だと、 １が連続した部分、０が連続した部分が、表示の上部、下部に 張り付いてしまって、信号が遷移する縦線しか見えなくなっちゃうんだ。表示範囲をもう少し、広く してくれれば、いいのに。
レベルが0.5,-0.5ってのは、いかにも中途半端だな。ここは、Voh = 3.4, Vol = 0.4 ぐらいな、TTLの規格と してはぎりぎりぐらいにしておいた方が良かったかも。えー、Vohとかの意味知らないって？
かの昔、SN74XXシリーズと言う、デジタルICが有ってだな、その規格表に書いてあるんよ。最初のVは、電圧 と言う意味ね。他には、I で、電流とか、T で時間とかあるよ。真ん中のo は、出力側端子って事ね。 他には、i で、入力側端子ってのがある。最期の h やｌは、ハイレベル、ローレベルって意味だ。
メッセージは、デフォで、JPL にしといた。CQ誌に載ってた、火星探検衛星から繰り出されるローバーの キャタピラーに、JPLのモールスコードが埋め込まれているそうなんで、宇宙に夢を馳せて頂いてきた。
コードの前半は、メッセージのドット数を算出し、そこから、ダッシュ分のデータ(ref)を作っている。 このデータの冒頭付近に手を加えれば、リンギングとか、キークリックとかが実現出来るんで、暇な人は やってみると面白いぞ。 後半部分では、refから、ドット/ダッシュ分のデータを所定の位置、サイズでコピーして、求める信号に仕上げています。 ここでも、スライスの機能が活躍してますよ。
変調器
お次は、変調器だな。２つ例を挙げる。（まだ、それしか作ってないから） 何はなくとも、振幅変調だな。
def am(fs, fc=300, m=0.5): # amplitude modurator
for n in range(SMPL):
fs[n] = (1 + m * fs[n] ) * cos(2 * pi * fc * n / SMPL)
return fs
もう一つは、振幅変調から、キャリアーを抜いてしまうやつ。ダブリサイドバンド。ああ、書き間違っちゃった。 ダブルサイドバンドが正解。でも、だぶり でも通じそうだな。
def dsb(fs, fc=300): # double side band
for n in range(SMPL):
fs[n] = fs[n] * cos(2 * pi * fc * n / SMPL)
return fs
amとdsbで何処が違うん？ と目を凝らして見れば、キャリアーの振幅部分で、被変調波(fs)に、１を加えて いるかどうかしか違わない。（変調率はこの場合本質ではないので、乗算の単位元１とみなせば良い） すなわち、amの場合は、被変調波がゼロだった場合でも、キャリアーは常に出ている。それに対しdsbの方は、 キャリアーが消失していて、変調波の影響しかない。
かの昔、807の終段のタンク回路の電源側に変調トランスを入れた、AM送信機を作った。マイクからの信号が 無くても、常にプレートには電圧がかかっていて、電波は出ていた。変調トランスと陽極の直流電源が直列に なってるんで、トランスを介して音声が加われば、プレート電圧はそれに合わせて高くなったり低くなったり する。それで、電波が強くなったり弱くなったりして、音声信号がキャリアーに乗って出てくんだ。 こうやって、式にしてみると原理が良く分かって面白いな。
写真集
AKB48の写真集じゃなくてごめんなさい。これから、各種の波形を開陳します。変調波は、シングルトーンじゃ ありきたりなので、４声のコーラスを使う事にします。
wb> vo | dt
まずは、被変調波となる４音のコーラス音源
wb> vo | am | dt
それを変調して。。。 こんな電波が発射されるんだな
振幅変調後の周波数特性。酷試の問題集にある通りだな。
wb> vo | am | fft | nb | gett | df
振幅変調後の位相特性、途中に入っている nb は、ノイズブランカです。ええ、fftした結果、非常に 小さな値が出てくるんですが、それの位相も表示されちゃって本信号が埋もれちゃうんです。そんな訳で 小さな複素数は、0+0j に丸めちゃってます。
サイン波の位相は９０度、コサイン波は０度だな。上側波帯(USB)と下側波帯(LSB)で、位相が反転するなんて、 初めて知りましたよ。こんなニッチな事は教科書にも載ってないものねぇ。大体、位相は無視される事が 多く、かわ位相。
dsbすると、キャリアーが消えるよ。今流行の省エネだな。更に省エネしたかったら、どっちかの側波帯を 削ってSSBだな。そして、そこから更に省エネを目指すなら。。。
wb> dd | dt
ドット ダッシュ発生器のお世話になるしか。
方形波の集まりは、いろいろな周波数成分が混じってるな。側波帯を広げない為、40wpmなんてのはもっての外です。 1wpmぐらいで、のんびりと行きましょう。
上記の写真集は、emacsからmyREPLである、wb()を起動して撮ったけど、一度もハングしなかった。どうやら 普通にpythonを使うと、shellとemacs間で、タイムアウトしちゃうんだな。疑ってゴメンな、emacs君。後で、 python.elでも眺めておくかな。眺めてみたら、C-c C-c なんてキーバインドを発見。こいつは便利。
現在のワークベンチ一覧
from pylab import *
from numpy import *
from numpy.fft import *
from math import *
from cmath import phase
SMPL= 4096 # sampling points, must be 2**n
r2d = lambda r: 180 * r / pi # radian to degree
d2r = lambda d: d * pi / 180 # degree to radian
th = lambda z: r2d(phase(z)) # angle
getr = lambda a: array([x.real for x in a]) # real
geti = lambda a: array([x.imag for x in a]) # imag
geta = lambda a: array([abs(x) for x in a]) # power
gett = lambda a: array([th(x) for x in a]) # phase
def nb(a, th=1e-4): # round to 0+0j when less th
for n in range(SMPL):
a[n] = 0+0j if abs(a[n]) < th else a[n]
return a
def dt(a): # display time domain
plot(a)
show()
def df(a): # display freq domain
#plot(a[:SMPL/2])
#bar(arange(SMPL/2), a[:SMPL/2])
bar(arange(SMPL/8), a[:SMPL/8])
show()
def st(fs1=3, g=1.0): # single tone (for audio)
wave = []
for n in range(SMPL):
wave.append(g * cos(2 * pi * fs1 * n / SMPL) )
return array(wave)
def tt(fs1=3, fs2=5, g=0.8): # two tone (for audio)
wave = []
for n in range(SMPL):
wave.append(g * (cos(2 * pi * fs1 * n / SMPL) +
sin(2 * pi * fs2 * n / SMPL)) )
return array(wave)
def vo(fs1=2, fs2=3, fs3=4, fs4=5, g=0.7): # like voice
wave = []
for n in range(SMPL):
wave.append(g * (0.4 * cos(2 * pi * fs1 * n / SMPL) +
0.8 * sin(2 * pi * fs2 * n / SMPL) +
0.6 * cos(2 * pi * fs3 * n / SMPL) +
0.2 * sin(2 * pi * fs4 * n / SMPL)) )
return array(wave)
def dd(msg='2133321331213116', m=0.5, s=-0.5): # dot dash (msg=JPL)
wave = SMPL * [s]
dots = 0
for n in msg:
if n == '1' or n == '3':
dots += (1 + int(n))
elif n == '2' or n == '6':
dots += int(n)
else:
pass
t = int(SMPL / dots)
ref = [m] * 3 * t
st = 0
for n in msg:
if n == '1' or n == '3':
wave[st:st+t*int(n)] = ref[0:t*int(n)]
st += (t * (int(n) + 1))
elif n == '2' or n == '6':
st += (t * int(n))
else:
pass
return array(wave)
def am(fs, fc=300, m=0.5): # amplitude modurator
for n in range(SMPL):
fs[n] = (1 + m * fs[n] ) * cos(2 * pi * fc * n / SMPL)
return fs
def dsb(fs, fc=300): # double side band
for n in range(SMPL):
fs[n] = fs[n] * cos(2 * pi * fc * n / SMPL)
return fs
def run(wire):
try:
exlude = set(['show()',])
cmds = wire.split('|')
for n, cmd in enumerate(cmds):
cmd = cmd.strip()
if (n == 0) or cmd in exlude:
do_cmd = (cmd if cmd[-1] == ')' else (cmd + '()' ))
else:
do_cmd = (cmd.replace('(', '(val,', 1) if cmd[-1] == ')'
else (cmd + '(val)' ))
val = eval(do_cmd)
return val
except Exception as e:
print("Error: {0}".format(e))
def wb(): #################### my REPL
while True:
run(raw_input('wb> '))
if __name__ == '__main__': wb()