2012-02-03
■[Arduino] 空中 Quine (空中に自分自身を出力するプログラム)
何を言っているかわからないと思いますので、まずは動画を見てください。
つまり、LED の残像として空中に自分自身を「出力」するプログラム (＋回路) です。こんな風に空中に字を描くのを POV (Persistance of Vision = 残像) とか Light writing とか言うらしいですね。やることは LED 光らせるだけなのに、やってみるとかなり楽しい。
というわけで、最近 Arduino とかいうマイコンで遊んでます。アナログ電気回路の知識は中学生レベル *1 ですが、USB で PC に繋いで IDE の upload ボタンを押すだけで C プログラムをマイコン上で実行できるので、まあそれなりに遊べます。
最初に作るのは当然 Quine なわけですが、シリアル通信で Quine 出力するのは既出だったので、ちょっと凝ってみました。もっと読みやすいコードにしたかったのですが、メモリが 2 KB しかないのでゴルフを余儀なくされた。フォントは @hirekoke さんに作ってもらいました。
回路は 2 〜 8 番ピンについて、ピン 〜 抵抗 〜 LED 〜 GND とつなぐだけです。Quine のコードの最後に書いてある通り *2 。自作シールドの作り方は先人の知恵を参考にさせてもらいました。
ところでこの Quine 、ちょっと費用がかかってます。工具から買ったから 2 万円くらい。Amazon アフィリエイトのリンクを置いときますので、よければ買って Quine してみてください。
次は Arduino から別の Arduino へ自己複製するものつくろうと思っていたら、これも既出だった。この業界は意外に Quine されてるなあ。
/* Arduino LED Quine */ /* (c) @mametter, @hirekoke, 2012 */ unsigned int x,y,a[6],b[6],c,d,n,N=128; extern char *p,T[]; void w(int v){ if(!v){for(x=1; y=digitalRead(13),x+!y; x=y); return;} for(x=0,b[0]=1; x<7; x++) for(y=c=d=0; y<6; y++) c/=N,a[y]=(c+=b[y]*((T[v*7+x+412]-1)%92)+a[y])%N, d/=N,b[y]=(d+=b[y]*91)%N; for(x=0; x<6; x++,delay(4)) for(y=b[x]=0; y<7; y++) digitalWrite(8-y,a[x]&1?HIGH:LOW),a[x]/=2; } void setup(){ for(x=2; x<14; x++) pinMode(x,x<10?OUTPUT:INPUT); for(p=T,w(0); *p!=36; p++) w(*p%126); for(p=T,w(34),n=12; *p; w(34),w(n=0),w(34+25*!*p)) while(n<53) w(*p++),n++; } void loop() { } char *p,T[]="/* Arduino LED Quine */~/* (c) @mametter," " @hirekoke, 2012 */~unsigned int x,y,a[6],b[6],c,d,n," "N=128;~extern char *p,T[];~void w(int v){~ if(!v){for" "(x=1; y=digitalRead(13),x+!y; x=y); return;}~ for(x=0" ",b[0]=1; x<7; x++) for(y=c=d=0; y<6; y++)~ c/=N,a[y]" "=(c+=b[y]*((T[v*7+x+412]-1)%92)+a[y])%N,~ d/=N,b[y]=" "(d+=b[y]*91)%N;~ for(x=0; x<6; x++,delay(4)) for(y=b[" "x]=0; y<7; y++)~ digitalWrite(8-y,a[x]&1?HIGH:LOW),a" "[x]/=2;~}~void setup(){~ for(x=2; x<14; x++) pinMode(" "x,x<10?OUTPUT:INPUT);~ for(p=T,w(0); *p!=36; p++) w(*" "p%126);~ for(p=T,w(34),n=12; *p; w(34),w(n=0),w(34+25" "*!*p))~ while(n<53) w(*p++),n++;~}~void loop() { }~c" "har *p,T[]=$]]]]]]]~#<_]]]l}mV_]]+sxH6^].f.r4^]E$*rQh" "aSBbl?H]qVd_]]]ZgZ__]],Cy^]]]D`'B+^]-$}5|]]>Qh]]]]g3m" "4|]]eX`]]]]<x>E#^]^eGQ8_]9yye]]]ro*]m_]V[5FT^]RY]1p]]" "W?3_x`]dK8.2]]k?|.ka]1s6FT^]K0:j8_]6qH]]]]doL]]]](jba" "^]]v#*SO]]BdAs]]]~CqC%^]|B$ZB+^zMs&9]](t6FT^]cFXW4^]R" "$Oba^]x(#QO_]Iw{[G_]Mmj0T^]GZzB'b]BB'a_]]]U9*{b]Nee+N" "_]RjaRd]]xpOo(b]usc,'b]w6nS7_]~i&ge_]/2AX@j]qHYbm_]_<" "5o5^]8cJ/G_]zV9*{b]&=MHKa]>p~Ukb]50e+N_]dsEvia]`Nljsa" "]Q@&a_]]cCDWd]]#yE_]]]Pc*[]]]k-8b]]]O5)[]]]SY{Cx^]Qwf" "10]]`RMMO]]>*VG'b]g>+ci^]^i$+_]]ZP(dx^]'jWw8]]iZ`_]]]" "{~9P_]]v9GcH]]IHxe]]]V-Qw8]]$=Uw8]]:vf10]]?4oba^]gH@B" "B_],Z}h@]]FU;p4^]*a.3]]]}fYDe]],|U[?]]L%U=o^]-~v^H]]Y" "x46w^]ViFtH]]vsYa^]]2wD_]]]4m9t]]]YaA-]]] " " ,----------------------------. " " | Arduino| " " | 2 3 4 5 6 7 8 9 13 Vcc GND | " " `----------------------------' " " | | | | | | | | | | | " " R R R R R R R R R B | R: resister " " | | | | | | | | | | | L: LED (common cathode) " " L L L L L L L L `--' | B: push button " " | | | | | | | | | " " `-+-+-+-+-+-+-+---------' " ;
2012-01-22
■[Ruby] 音声の波形からピッチを検出するアルゴリズム
去年のクリスマスに公開したカラオケ機能つき Quine の仕組みについて。
ref: 声の高さで操作するゲームを作ってみた
で解説されている内容と同一です。おわり。
で終わるのもつまらないので、簡単に解説します。でも思いだしながら書いているので嘘書いてたらごめんなさい。動画には図とかあるので、やはりそっち見た方がいいと思うけど。
「ピッチ検出なんて FFT するだけでしょ」と思ってる人は素人で、音叉みたいにきれいな正弦波を測りたいならともかく、声や楽器の音など倍音を含んだ音では誤判定が起きまくるようです。偉そうなこと言ってる私も素人です。そこで、Wikipedia の Pitch detection algorithm で挙げられている、MPM アルゴリズムを調べて実装してみました。以下の論文。
ref: P. McLeod and G. Wyvill. A smarter way to find pitch
波形から自己相関関数を計算する
自己相関関数とは、波形を時間軸方向にずらしてみて、元の波形とどのくらい一致しているか、という指標。その定義通りに計算すると非効率ですが、不思議な力によって、FFT 2 回で効率的に計算できるそうです。ちなみに Ruby で FFT は
F=proc{|n,t,a|n/=2;r=(0..n-1);n<1?(a):(F[n,t*2,r.map{|i|a[i]+a[i+n]}].zip(F[n,t*2,r.map{|i|(a[i]-a[i+n])*Complex.polar(1,t*i)}]).flatten)}
というように書けます (参考) 。定数オーダの効率はよくないですが、カラオケ判定程度ならぎりぎり実時間でできるようです。(8000 bps で 1024 点ずつ処理。1 秒あたりおよそ 8 * 2 = 16 回くらい FFT してる)
MPM アルゴリズムでは、自己相関関数を基にした Normalized Square Difference Function (NSDF) とかいう関数を使います。-1 から 1 に正規化されているし、インクリメンタルに計算できるし、何かと都合がいい関数のようです。
nsdf = F[N*2,-T,F[N*2,T,w+[0.0]*N].map{|v|v.abs2}].map{|v|v.real/N/2}[0,N].map{|v|v=2*v/m;m-=w[i]**2+w[-i+=1]**2;v}
NSDF のピークを見つける
波形のピッチを求めるには、元の波形とずらした波形の一致度が高いところ、つまり NSDF の値が大きいところを求めればいいわけです。しかし、誤検知を防ぐためのバッドノウハウのところのようで、なかなか dirty になっています。
まず波形を時間軸方向になめていって、値が負から正に入ってから、正から負に落ちるまでの間で、もっとも値の大きなところ (ピーク) を列挙していきます。単純なピーク (微分したら 0 になる地点) ではないことに注意。
さらに、各ピークの前後の値を見て、放物線補完を使ってより正確にピークの時刻・大きさを特定する。
そうして列挙されたピークの中で、最大ピークの大きさの 8 割以上の大きさを持つピークで、一番最初に現れるものを採用します。このピークが現れる時刻 (つまり波長) の逆数が、求めるピッチの周波数。
また、ピークの値の大きさは 0 から 1 に正規化されていて、ピッチ検出の確度として使えます (大きいほど確実) 。
図がないとわからないでしょうね。先の動画や論文を見てください。ソースコードで言うとこの部分。k が波長、v がピークの大きさになります。
nsdf.each_cons(3){|x,y,z|j+=1;x>0&&(0>=y&&(e<<[k*11,h];g<h&&g=h;h=0.0);y-x>0&&z-y<=0&&(d=(x-z)/t=2*(x-2*y+z);(h<c=y-t*d*d/4)&&(k,h=j+d,c)))};k,v=e.find{|i,v|v>g*0.9}
ピッチから音階に直す
音の波長から音階に直すのは、log で計算できます。440Hz がオクターブ 4 のラの音なので
Math.log(440/k,2) * 12
の mod 12 が 0 ならラ、1 ならラ# 、2 ならシ、……、11 ならソ#、という感じです。
まとめ
音声からピッチを検出する MPM アルゴリズムについて、せっかく調べたのでメモがてら簡単に書いてみました。
でも、図なしで説明してもわからないですね。動画の方がいくぶん分かりやすいです。詳しく知りたければ論文みましょう。門外漢にも分かりやすくて self-contained で、4 ページと短いアルゴリズム論文で、すばらしいです。
しかしこんな解説よりデモ動画の方が必要なんだけど、歌声とか公開したくないのだった。誰か歌って。
2012-01-01
■[Ruby] あけましておめでとうございます 2012 - なめらかライフゲーム
あけましておめでとうございます。日本のプログラマには古来より「正月はフラクタル」という習わしがあります。正月はフラクタルに触れて心穏やかに過ごそうというものです。
というわけで動画を作りました。何かわかるでしょうか。
答えは、ライフゲームのグライダーとブロックの衝突です。セルをそのまま四角で表示するのではなく、アニメーションでなめらかにつないで表示してみました。
グライダー銃 (上) と宇宙船 (下) をグライダーから作る様子。
銀河 (左) とヒキガエル (中央) とパルサー (右) 。かなりきもい。
まあ、ライフゲームはフラクタルかというと違う気がしますが気にしない。一応フラクタルらしいドラゴン曲線も。辰年だし。
実装を簡単にいうと
- 普通のマスの代わりに電荷みたいなものを配置して、ポテンシャル場を計算する
- ポテンシャル場に対して marching squares で等高線を描く
という感じですが、いずれ詳細をきちんと書きたい。ソースコードも意外とぐちゃぐちゃになったので、整理してから公開したい。
あと、参考にしたもの。
2011-12-24
■[Ruby] Merry Quine-mas 2011
今夜は W. V. Quine の命日です。
eval$s=%q(eval((%w(A="+_Sing_do,_re,_mi_to_the_mike!*z+_Do!*{+_Re!*|+_Mi!__OK
!_Let's_sing_Silent_Night!*~+_Sxi~lent_|night!_~Hxo~lY_|night!*{+_All_is_ycal
m,_zall_is_~bright,*x+_Round_You_zViyrxgin_~Moxther_~and_|Child!*x+_HolY_ziny
fant_xso_~tenxder_~and_|mild,*{+_Sleep_in_}heav{enylY_zpea|ce,*z+_Slee~p_|in_
~heav}en{lY_zpeace.**";T=-Math::PI/N= 2**10;F=proc{|n,t,a|n/=2;r=(0.
.n-1);n<1?(a):(F[n,t*2,r.map{|i| a[i]+a[i+n]}].zip(F[n,t*2,r
.map{|i|(a[i]-a[i+n])*Comple x.polar(1,t*i)}]).flat
ten)};u=open("/dev/ds p","rb"){|s|loop{u||=
[1,3,4,6,8,10,11][( A.slice!(/([!-w]+)([x
-~])?/);$><<$1.tr(" Y*_+","y"<<10<<32<<35
);$2||(c=1614776052 204941692509352;s=["e
val$s=%q("+$s+")",r= "#" *2,"Yusuke","Endoh,",
2011,r]*32.chr;s[i=34 9]. ord>31&&"5i=bFV0IS3
HS4GS5FT6>-.U7<-0T8:,2 T: 7,2V:5-0Z;0f>g<j8*
".each_byte{|d|s[i,0]=" \e [%dm"%[0,41,47][c%
3];c/=3;i+=d-38};puts(s); exi t);$2.ord-121)];i=j
=k=m=0;w=s.read(N).unpack("C *").map{|v|v-=128.0;m+=v*v*2
;v};g=h=0.0;e=[];F[N*2,-T,F[N*2 ,T,w+[0.0]*N].map{|v|v.abs2
}].map{|v|v.real/N/2}[0,N].map{|v |v=2*v/m;m-=w[i]**2+w[-i+=1
]**2;v}.each_cons(3){|x,y,z|j+=1;x>0 &&(0>=y&&(e<<[k*11,h];g<h&&g
=h;h=0.0);y-x>0&&z-y<=0&&(d=(x-z)/t=2*(x-2*y+z);(h<c=y-t*d*d/4)&&(k,h=j+d,c))
)};k,v=e.find{|i,v|v>g*0.9};k&&k>0&&v>0.9&&99<k&&k<2E3&&(((u-Math.log(25/k,2)
*12)%12).abs<2.0&&u=$a)}})*"").gsub(/\e\[\d+m/,""))) ## Yusuke Endoh, 2011 ##
カラオケ機能のついた Quine です。/dev/dsp とマイクが必須です。
実行したらまずド、次にレ、最後にミの音を出してみて、音を確認してください。そのあと、「きよしこの夜」をゆっくり歌ってください。音に合わせて歌詞が出るはず (バッファリングのせいで結構遅れます) 。
音が外れたり認識できなかったりすると歌詞が進みませんので、そこから歌いなおしてください。最後まで進んだら Quine になります。
動画は後日。
2011-11-24
■[Ruby] Ruby で FFT (高速フーリエ変換) を書いてみた
ref: 【ニコニコ動画】ミクをPCの再生音に合わせて自動で踊らせてみた
↑に触発されて波の処理をしたくなったので、Ruby で FFT (高速フーリエ変換) を書いてみました。
FFT とは、波の形を見て周波数とかを見抜く魔法のことです。数式とか考えたくないので、とにかく Ruby で書いてみました *1 。
def fft(a) n = a.size return a if n == 1 w = Complex.polar(1, -2 * Math::PI / n) a1 = fft((0 .. n / 2 - 1).map {|i| a[i] + a[i + n / 2] }) a2 = fft((0 .. n / 2 - 1).map {|i| (a[i] - a[i + n / 2]) * (w ** i) }) a1.zip(a2).flatten end
これだけです。短いですよね。
実用上はここから in-place 化とか scramble とか定数オーダの最適化が待っているらしいですが、Ruby で書いてる時点でそんな方向には興味がありません。*2
そんなことよりちょっと動かしてみましょう。まずは入力用の波形を作る。さっきの fft の実装では、サンプル数が 2 の累乗でないとダメです。
N = 64 a = (0...N).map do |n| v = Math.sin(2 * 2 * Math::PI * n / N) * 2 v + Math.cos(5 * 2 * Math::PI * n / N) end
2 Hz の sin 波と 5 Hz の cos 波を合成した波形。2 Hz の振幅は 2 倍。図示して確認すると
a.map do |v| s = [" "] * 20 min, max = [(-v * 3 + 10).round, 10].sort s[min..max] = ["#"] * (max - min) s end.transpose.each do |l| puts l.join end
#
###
## #####
#### #####
###### #######
### ####### #######
################ #########
################ #########
################ ########### #
################# ###########
################ #########
################ #########
####### ### #######
###### #######
#### #####
## #####
###
#
なんかよくわかりませんが、できてるのかな？ これを見て、2 Hz はなんとなくわかりますが、5 Hz が混ざってることはちょっとわからないですよね。さっきの FFT にかけてみます。
fft(a)[0, N/2].each_with_index do |v, n| next if n == 0 puts "%2d Hz: %.3f" % [n, v.abs / (N / 2)] end
1 Hz: 0.000 2 Hz: 2.000 3 Hz: 0.000 4 Hz: 0.000 5 Hz: 1.000 6 Hz: 0.000 7 Hz: 0.000 8 Hz: 0.000 9 Hz: 0.000 10 Hz: 0.000 11 Hz: 0.000 12 Hz: 0.000 13 Hz: 0.000 14 Hz: 0.000 15 Hz: 0.000 16 Hz: 0.000 17 Hz: 0.000 18 Hz: 0.000 19 Hz: 0.000 20 Hz: 0.000 21 Hz: 0.000 22 Hz: 0.000 23 Hz: 0.000 24 Hz: 0.000 25 Hz: 0.000 26 Hz: 0.000 27 Hz: 0.000 28 Hz: 0.000 29 Hz: 0.000 30 Hz: 0.000 31 Hz: 0.000
ばーん。2 Hz と 5 Hz を見事見破ってくれました。fft の返り値の見かたのポイントは、
- 意味のある情報は配列の前半だけ
- 0 番目の値も無視する
- 値は複素数になっているので大きさ (v.abs) をとる
- (N / 2) で割った値がその周波数の強さ
- ちなみに位相の情報は偏角 (v.angle) として入っているらしい (あんま見てない)
という感じです。*3
*1：「FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法」 の「リスト1.2.1-1. 再帰呼び出しによる FFT」の移植です。
*2：まじめに Ruby で FFT とかしたい人は NArray の FFTW3 とか使うといいんだと思います。ぼくはやったことありません。
