Workbench(3)
防衛省で、球形飛行体 なんてのを研究してるらしい。市ヶ谷のあそこですかね。昔はおいらも あの辺へよく通ったものだ。 あの辺は、怪電波が出てるらしくて、一緒になった元東大の教授様が、電波腕時計があのあたりへ行くと 狂うとおっしゃっていたっけ。怪しさ満載の市ヶ谷・四ツ谷エリアですよ。
だけど、あんな飛行物体が秋葉原で材料を仕入れて作れちゃうって、怖いな。あそこに灯油と燐でも混ぜた たいまつを吊るして、お墓の回りを旋回させたら、涼しくなるぞーーー。防衛省さんよ、今年の夏は、節電 グッズとして、貸し出したらどうでっしゃろ？
ネーミングは、これで涼しさ満点、クールスポット ぐらいでいいかな。防衛省のお膝元には、お岩さん とか三島さんとかがおられるはずですから、これらの方とタイアップすれば、もう万全でっせ。
こちらは、噂によると例の爆発現場を撮影した飛行物体だそうです。カメラの代わりに、O-104でも搭載して、 ばら撒くとかは、考えないで下さい。
どうせなら、かわいいやつがいいです。上と同じ会社で作っているそうです。おいらも一羽欲しい。餌は 何でしょうね？ 電気を食べているんでしょうか？小さくする事が得意な日本の技術を持ってすれば、蚊 ぐらいのサイズに出来るかな？ 京で世界一になったけど、マイクロマシーンとかMEMS分野でも、小さいので 世界一を目指してくらはい。
AMPM
前回は、振幅変調器とかを作ってみた。まだまだ電子工作は続くよ。今度は周波数変調をやってみるか。 周波数変調って、キャリアーがふらふらするやつね。昔はVFOを自作したけど、なかなか安定しなくて 苦労した事を覚えている。VFOが温度計になっちゃりしてね。そう、温度で発振周波数が勝手にドリフト しちゃうんだ。
で、そのドリフトを、音声の大小によって能動的に変えてしまえば、周波数変調器になるんだ。 電子工作的には、次のようになるかな。
def fm(fs, fc=300, m=0.5): # freq. modulation
for n in range(SMPL):
fs[n] = cos(2 * pi * (fc + m * fs[n]) * n / SMPL)
return fs
FMと言うとステレオを思い出す。おいらまっとうなクラッシックファンで、よくダビングしたっけ。ステレオ は、左右のスピーカーから別々な音が出てくるけど、どうなってるんだ？ この際だからちょいと調べてみた。
FMステレオ方式に よると、AM-AM方式と和差方式とスイッチング方式があるんだな。面白い。
AM,FMと来たら、次はPMだな。PMが揃えば、AMPMでコンビニが一丁出来上がり。簡単だから、ここでは いちいち表示しないけど、アマチュアの世界でPMって使われているんかいな？ こちら経由で なじみのあるのをを 見ると、確かに許されてますな。ファクシミリがそうか。最近は聞かないけど、昔はやってる人が居たような。。。。
CW送信機
各種の変調器が出来たんなら、A1Aの信号も発生したいよね。A1Aの送信機を作ってもいいんだけど、ちょいと 汎用的に考えてみる。パドルのon/offによって、搬送波が出たり出なかったりすればいいわけだ。
見方を変えると、常に発振してる搬送波とパドル信号のAND（論理積）を取ってあげれば良い。言い方を 変えると、搬送波をパドル信号でゲーティングするって事だ。そうすると、２入力のAND回路か。はて、 パイプの途中に信号を注入するんだな。
きっとA1Aクラブの人は、これを『ラブ注入』って言うんだろうな。
途中に注入するには、AND回路の一方の入力端子に、信号の発生手順を教えてあげれば良い。 勿論もう一方は、パイプの 上流から流れてくる信号ね。AND回路は、掛け算で実現出来るな。
mix = lambda a,b: array((a + b)) # add array mul = lambda a,b: array((a * b)) # product array
ついでなので、足し算も定義しておいた。 一応、array間の演算動作を、Pythonのリストと比べてみる。まずはリストから
>>> a = [1,2,3,4]
>>> b = [10,11,12,14]
>>> a + b
[1, 2, 3, 4, 10, 11, 12, 14]
>>> a * b
Traceback (most recent call last):
File "<pyshell#6>", line 1, in <module>
a * b
TypeError: can't multiply sequence by non-int of type 'list'
>>> a + 3
Traceback (most recent call last):
File "<pyshell#14>", line 1, in <module>
a + 3
TypeError: can only concatenate list (not "int") to list
>>> a * 3
[1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4]
一杯、エラーを喰らっちゃったよ。まあ、一部にLispっぽいのが混じっているね。それに対して、リストをarrayに変換して おくと。。。
>>> x = array(a) >>> y = array(b) >>> x + y array([11, 13, 15, 18]) >>> x * y array([10, 22, 36, 56]) >>> x * 5 array([ 5, 10, 15, 20]) >>> x + 3 array([4, 5, 6, 7])
行列の出来る、numpyって事ですね。（事前に、numpyをimportしておく事）
wb> st(fs1=300) | mul(dd(m=1,s=0)) | dt
シングルトーン発振器から高い周波数（キャリアー）を得て、それとキーイング情報を掛け算する。 気分は、グリッドバイアス・キーイングだな。
１は何回掛け算しても、元の値と同じと言う、掛け算の 単位元です。ゼロを掛けると元がどんな値でもゼロになるって言う、数学では当たり前の事実を利用して います。単位元なんて言葉は、先の東大教授様が訳されたSICPを読んだ時に出てきました。それ以来、 単位元、それから式年遷宮とGCの繋がりを忘れる事はありませんでしたよ。
この２０年毎に繰り返される、お伊勢様の神社移転は、技術の伝承を保障するために編み出されたしきたり と言うのが、もっぱらの通説になっているそうです。で、GCが何で関連してくるかと言うと、STOP and COPY によるGCなんですね。言われてみれば、その通り。例のGC本には書いてあるのかな？ ああ、無駄話だったわい。
盗み見
パイプの途中に信号を注入する方法は分かった。そうなると、パイプの途中にどんな信号が流れて いるか、盗み見 したくなる。（今までの波形観測機、dt,df は、パイプの終端でしか機能しなかった）
unixのパイプの先生に聞いてみると、tee じゃよと。ついでに、パイプの先端にファイルからデータを 注入する機能や、データをファイルに落とす機能が無いぞと指摘されちゃったよ。嗚呼、やぶ蛇に なっちゃったよ。
そこで、先生には、ファイルに落とす演算子（unixのshellでは、有名な、＞ 記号ね）は、名前を 考えてあります。その名前は、セーブオンですって、コンビニ繋がりで咄嗟に思い出したのを 口走ってしまったよ。saveon('hoge.txt') なら、語感がいいからね。でも、もう一つ有名な、＜ に変わる やつは、いい名前が思い付かない。ここは、名前が決まらないと実装しないと言う、あの人に習って おこう。それより、盗み見の方が、需要だよ。
あの人に習って、まず名前を決めよう。teeでもいいけど、アマチュア（無線）的な名前の方が親しみが 有っていいだろう。機能的には、パイプに挿入して信号をモニターするんだから、方向性結合器、サーキュレーター になるんかな。
トランシーバーとかの基板を見ると、要所要所に信号モニター用のテストポイントが立っていたりする。いわゆる tpか。まあ、それでもいいんだけど、ちょっと悪ぶって、tapなんてのはどうだろう。
どうせtapするなら、信号は自動検出したいね。実数が流れていたら、ああ、タイムドメインで表示するんか。 複素数なら、周波数ドメインで表示だ。周波数ドメインだと、振幅特性と位相特性があるな。どちらを 表示させる。SWを付けて選ばせるか？ 自動検出までやってくれるなら、両方同時に表示させちゃえ。 そうすれば、考える事いらんわな。
まあ、自動検出って考えた時点で、信号の推測を放棄しちゃってる訳 だから、そんなの今のお前には10年早いわと言われれば、返す言葉は有りませんが。。。デジタルオシロとかが 出てきて、信号を推測しながら観測するって風潮が失われたのは、今の電子業界に取って、痛し痒しです。 ついでに書いておくと、放射能を急に測り出して、大きい小さいって一喜一憂してるけど、3.11前の値を 各地で測ったデータって有るんでしょうかね？
def tap(a): # tapping
figure()
if True in iscomplex(a):
#subplots_adjust(hspace=0.001)
subplot(211)
f = geta(a)
bar(arange(SMPL/8), f[:SMPL/8])
subplot(212)
p = gett(nb(a))
bar(arange(SMPL/8), p[:SMPL/8])
else:
plot(a)
return array(a)
def disp(a): # put on pipe end
tap(a)
show()
iscomplexで複素数を判定する時、最初は、a[0] を見てたんだ。3.12+0j みたいな虚部が ゼロの値が出てきて、Trueにならないのね。しょうが ないので、上記のようにしましたよ。でも今考えたら、a[0]って、直流成分の所じゃん。直流分に 虚部があったらまずいんで、積極的に0jにしてるんかな？
それから、今はコメントにしてある部分を生かすと、上側のグラフ （の表示部分）と下側のグラフのくっついてくれます。
自分用のメモとして、subplot(abc) の、aは、縦方向のグラフ数、bは横方向のグラフ数、cは描画対象の グラフ番号を表している。
dispは、パイプの終端に入れるターミネーターです。これを入れておかないと、SWRが無限大になってしまい、 保護機能が働いて、結果が表示されないままパイプが終了してしまいます。終端は必ずターミネート しましょう。
ろ過器
次は、福島方面で活躍してる、ろ過器を取り上げます。不要なのを除去して必要なものだけ取り出すアレです。 無線機の中には、水晶で出来たのがいいですぜと、メーカーからの甘い誘いにそそのかされて、高価なやつを いっぱい入れてる人もいるでしょう。高いのを買うと運が向上して、いろいろな賞状がホイホイと手に 入るんでしょうか。 だとしたら、一本３０万円の開運印鑑と一緒かな？ お年寄りは注意して、甘い言葉には騙されないように しましょう。
世の中には、いろいろなタイプのろ過器が出回っています。たとえば、 料理本に出てた、 スムースとか、 FIRさん、 カルマン方式、 等。ろ過すれば、もっと見つかるかも知れないよ。
で、おいらは、金が無いので、手を動かします。 まずは、ノイズ交じりの信号だな。ノイズをmixすればいいんだな。
wb> st | tap | mix(tt(fs1=39,fs2=51)) | tap | fft | disp
入力信号をモニターし、ノイズ相当（周波数が高い繰り返し波形ですが）を加えてからモニターし、最後に周波数領域に変換して、スペクトラムを 観測。まあ、ここまでは、理論通りだな。
それじゃ、フィルターを作るか。周波数領域でフィルターする事を意図します。一応ローパスフィルターをば。
def lf(a, fc=20): # low pass filter
for n in range(SMPL):
a[n] = 0+0j if fc <= n < (SMPL - fc) else a[n]
return array(a)
カットオフ周波数２０のローパスフィルターです。２０以上の周波数データをゼロにして、無い事にして います。
wb> st | mix(tt(fs1=39,fs2=51)) | tap | fft | lf | ifft | getr | disp
シングルトーンの波形にノイズを混ぜ合わせ、その信号をfftして周波数領域へ持って行き、そこでフィルターを適用。 後は、ifftで時間領域へ戻して実部（の振幅成分）を 取り出すって事をやってます。
綺麗にノイズが消えて、正弦波に戻りました。そんじゃ、純正？なノイズを加えてみます。ノイズ発生器 は、まだ作ってなかったのでバラックで組んだ。本当は放射線の崩壊による方がいいんだけど。 ソフトウェアデザイン誌（SD)でもガイガーカウンターが取り上げられるぐらいだから、すっかりメジャーに なっちゃね。
def rnd(g=1): # random data
import random
wave = []
for n in range(SMPL):
wave.append(g * (random.random() - 0.5))
return array(wave)
wb> st | mix(rnd()) | tap | fft | lf | ifft | getr | disp
これ、ちょっと歪んでいるな。THDはどれぐらいだろう？
wb> st | mix(rnd()) | tap | fft | lf(fc=5) | ifft | getr | disp
フィルターの帯域を狭めてやると、ノイズの影響は軽減されて、元の波形に近づきますよ。
今日の実験台
だんだんとジャンク品が溜まってくるなあ。捨てるにはもったいないんで、まだそのままにしてるよ。
そうそう、変調器の出来具合も載せておく。シングルトーンを被変調波にした 周波数変調波と 位相変調波です。 上の周波数が512で打ち切りになってるのは、表示範囲を狭めているから。位相変調の位相が 拡大されているのは、オートスケールのおかげです。
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
mix = lambda a,b: array((a + b)) # add array
mul = lambda a,b: array((a * b)) # product array
def fm(fs, fc=300, m=0.5): # freq. modulation
for n in range(SMPL):
fs[n] = cos(2 * pi * (fc + m * fs[n]) * n / SMPL)
return array(fs)
def pm(fs, fc=300, m=0.5): # phase modulation
for n in range(SMPL):
fs[n] = cos(m * fs[n] + 2 * pi * fc * n / SMPL)
return array(fs)
def lf(a, fc=20): # low pass filter
for n in range(SMPL):
a[n] = 0+0j if fc <= n < (SMPL - fc) else a[n]
return array(a)
def tap(a): # tapping
figure()
if True in iscomplex(a):
#subplots_adjust(hspace=0.001)
subplot(211)
f = geta(a)
bar(arange(SMPL/8), f[:SMPL/8])
subplot(212)
p = gett(nb(a))
bar(arange(SMPL/8), p[:SMPL/8])
else:
plot(a)
return array(a)
def disp(a): # put on pipe end
tap(a)
show()
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
figure()
plot(a)
show()
def df(a): # display freq domain
figure()
#plot(a[:SMPL/2])
#bar(arange(SMPL/2), a[:SMPL/2])
bar(arange(SMPL/8), a[:SMPL/8])
show()
def rnd(g=1): # random data
import random
wave = []
for n in range(SMPL):
wave.append(g * (random.random() - 0.5))
return array(wave)
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', '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()