2017年7月3日月曜日

Karplus-Strong合成 初歩的なギター音の物理シミュレーション

シンセサイザの合成方式に,物理モデル音源という方式があるのをご存知でしょうか?(というより,シンセサイザの合成方式,いくつ挙げることができますか?)

名前から想像がつくと思いますが,手本にする音源(楽器である場合がほとんど)の発音原理を物理的にシミュレーションする方式です.この方式を初めてシンセサイザとして実装したのがヤマハです.



論理的に考えれば,この方式によってかなりリアルな音が得られそうな気がしますが,実際言うほどリアルな音は得られず(業界では分析が足りなかったと言われています),その後同社が物理モデル音源を前面に推した商品を出すことは無くなりました.しかし,エレクトーン等で管楽器の息遣いも含んだ質感のサウンドを演出したりするために使われるなど,物理モデル音源は今や電子楽器ではあって当たり前の技術となりました.

Rolandが数年前発表したV-Pianoは,ピアノの弦の振動からハンマー反響板に至るまで全ての要素をシミュレーションするだけでなく,パラメータ調整によって様々な物理的状態のピアノを作り出すことができます.


こうした物理モデル音源の名器達は,人間の優れた英知とそれを実装するに足る半導体の進歩によって可能になったものです.こうした物理的な分析に基づく音源合成は,1983年に発表されたKarplus-Strong合成が最初です.この合成法は,考案者の名前が元になっていて,KarplusさんとStrongさんによって考案されたことがわかります.Strongさんは何やら強そうです.

Karplus-Strong合成は,撥弦のシミュレーションです.「撥弦」とは,弦をはじく演奏法の事です.恐らく,この演奏法を使う最も知名度のある楽器はギターで,Karplus-Strong合成自体,撥弦の中でも「エレクトリック・ギター」のシミュレーションと紹介されることもあります.物理的にシミュレーションするというと,非常に難しそうな印象を受けますが,Karplus-Strong合成は,非常に簡単に作ることができます.まず,弦の振動について考えてみます.

皆さん既にご存知のように,音というのは振動現象です.音源の振動が波として空気を伝わって(例えば)耳まで届いて私たちは音を聞く訳です.波が伝わっていくもののことを「媒質」と言います.今の音の説明の場合,空気が媒質となっている訳です.ところで,この音の正体たる振動を伝える媒質は,空気に限らないことが知られています.例えば水や木の板等,振動を伝える媒質は空気以外にもたくさんあります.例えば,こういうものも媒質になります.

これ,学校にありませんでしたか?ウェーブマシーンと呼ばれる実験装置です.魚の骨のように,たくさんの棒が一本の軸で繋がっていて,棒の間を波が伝わっていきます.ここで重要なのは,「繋がっている」ということです.物凄く雑に言ってしまえば,それが何であれ,繋がっていれば,よほど凝集されていない限り,振動を伝える媒質として機能します.ですから,こんな風に「繋がった」物も波を伝える媒質になる訳です.


もう何が言いたいか分かりますね.私たちは弦の振動について考えるとき,大雑把に「弦が振動する」と考えてしまうかもしれません.でも,弦をはじく時のことを考えてみると,弦を指やピックで弾くとき,はじかれているのは,弦全体の局所的な部分です.そのときの振動が,弦という媒質を伝わって,やがて弦全体が振動するようになるのです.

この時,具体的にどのような現象が起きるのかを知るために,弦の振動をシミュレーションで図示してみましょう.

媒質を伝わる波の様子は,波動方程式


\[\frac{\partial^2u}{{\partial}t^2}=c^2\frac{\partial^2u}{{\partial}x^2}\]

で表現できます.この方程式をFDTD法で計算します.

尚,このシミュレーションは,神戸市立工業高等専門学校長谷芳樹准教授がwebで研究成果の一部を公開している

FDTD Simulation Movie & Demo

を参考にしました.上記ページで紹介されているのは2次元以上の波動方程式ですが,媒質の違いによって速度が変化する様子等が分かりやすく実装されているので非常に勉強になります.私のような現代のゆとり世代にも馴染みやすいツールで実装されているのも非常にありがたいです.

とりあえずやってみます.尚,上記の式では,減衰が考慮されていません.弦に限ったことではありませんが,現実的な環境で物が振動する時,何かが(地球上では大抵空気)抵抗となって,振動が時間の経過とともに減衰していきます.今回のシミュレーションでは,無理やり減衰をつけました.厳密な物理的解析に基づく結果ではありませんが,何が起きているか視覚的に伝われば幸いです.


何やら弦が振動してるっぽい動きをしてますね.初期状態(シミュレーション開始時点で弦がどういう形をしていたか)によって振動の仕方が様々に変化します.



ただ,これではよく分からないので,初期状態での変位(波の盛り上がり)をもう少し局所的にしてみます.

お分かりいただけますでしょうか?中央付近で発生した波が弦の左右へ伝わっていきます.山が両端まで到達したら,進行方向を逆にして跳ね返ってきます.跳ね返った波はいずれ再び中央付近へ戻ることになります.そして中央付近へ戻ったらまた両端まで伝わっていく・・・これを何度も繰り返しているのが弦の振動です.これを弦全体ではなく中央付近のみに注目して考えると,同じ波が両端の反射によって何度も到達することになります.この,同じ波が何度も到達する様子をシミュレーションしたのが,Karplus-Strong合成です.

MaxでKarplus-Strong合成をやってみましょう.

まず,撥弦を作ります.撥弦という動作は,弦にほんの一瞬だけ力学的刺激を与える行為だと考えることができます.ですから,パルスでこそないものの,パルスのような短い時間で変化する信号が適しています.


3 msだけ持続するホワイトノイズを撥弦とします.この時点でButtonをクリックすると,そのままですが,非常に時間の短いノイズが聞こえると思います.

次に撥弦によって弦に与えた刺激が,両端の跳ね返りによって何度も到達する様子をシミュレーションします.これは,ディレイを使います.

Maxのディレイのためのオブジェクトtapin~とtapout~は,tapin~のアーギュメントによって,ディレイの最大時間(ディレイのために確保するメモリ容量)を変えられます.今回はそんなに長いディレイタイムは必要としないので,100(ms)としておきます.「何度も到達する」訳なので,ディレイにフィードバックもかけます.とりあえず,ゲインは0.99にしておきます.

この画像のようにディレイタイムを10.msにして,Buttonをクリックするとどうでしょうか?先ほどのノイズがディレイとフィードバックによって何度も繰り返されているのがわかると思います.ここで更にディレイタイムを短くすると,徐々にピッチがはっきりと感じられると思います.勘の良い方はお気づきと思いますが,3msのノイズが周期的に繰り返されることによって,音源がノイズであったにもかかわらずピッチが生じているのです.これがKarplus-Strong合成です.

最後に,KSliderもつけて,少しばかり演奏インフェースっぽくしてみましょう.MIDIのNote Numberから周波数へはmtofで変換することができます.周波数からmsへの変換は,1000msに何回振動するかというのが周波数の定義なのですから,1000を周波数で割り算すれば求めることができます.



いかがでしょうか?物理モデリングといってもかなりシンプルに実装できることがわかります.初歩的な物理モデリングといえば,Karplus-Strong合成が最も有名ですが,実はこの考え方を応用したかなりたくさんの種類の物理モデリング音源が提案されています.色々な場所でサンプルとして見かけることができますから,ぜひ探してみてください.

※番外編

tapout~で設定するディレイタイムはsignalとして入力することで,「連続的」にディレイタイムを変化させることが出来ます.これを利用することで,チョーキングやビブラートのような効果を作ることが出来ます.


2017年6月27日火曜日

働きたくないので働かないアリのシミュレーション(もどき)

突然ですが,働きたくないんです.ないんです働きたく.

働いている人なら大体共感していただけると思うのですが,「俺は仕事大好きでさ,もう生きがい!!」なんて人はそんなにはいないんじゃないかと思います.皆毎日行きたくない職場へ行って,したくない仕事をしている面が少なからずあると思います.でも,私の働きたくなさはそんなもんではありません.よく,人よりちょっと顔が大きい人や,大きなほくろが顔にある人,目が大きい人なんかがいますね.あんな感じで,私は人よりも「仕事をしたくない」という感情が強いようです.この事は働き出してから気が付いたことで,理由は分かりません.これは,今やっている業種に不満があるとか,そういうことに原因があるわけではありません.「仕事」という営みをしたくないんです.だから,例えば今の職場を退職して,もっと「自分に向いている!!天職だ!!」と思う職業に就いたとしても,同じように仕事をしたくないと思うのでしょう.

そんなことを職場のデスクで考えていたら,以前のNHK 教育テレビ「サイエンスZERO」で働きアリについての研究が取り上げられていたのを思い出しました.

働きアリの中に一部働かないアリが発生するのは何故なのか?

そう,あるアリの一族の働きアリの中に,全く働かない,怠け者のアリがいるのです.これ,シフト制で一部のアリが交代で休憩しているのかと思いきや,そうではない.必ず決まったアリが働いていないのです.更に,驚くべきは,この一族から働かないアリだけを取り除いた時,全てのアリが働くようになりそう(働かないアリはいなくなった訳だから)ですが,そうではない.それまで働きアリであったアリの中から,働かないアリが発生するのです.

この事はそこそこ昔から知られていたことでしたが,その理由がついに判明したという研究です.

働かないアリも集団維持に必要 北大研究者が興味深い研究成果

この研究結果の概要を簡単に書くと,以下のようになります.

働きアリが全員働くと,全ての働きアリが一斉に疲れてしまって,最終的に早くリソースが底をついてしまう.一方,一部の働かないアリがいると,働きアリが疲れた時や,想定外の新しい仕事が現れたときに,その働かないアリが働くことで,リソースを継続することが出来る.

働かないアリは怠けていたのではなく,緊急時に備えて待機していたのです.

私も働きたくない人間として,この精神を見習わなければなりません.そこで,上記の原理を計算機シミュレーションで作ってみます.上記研究でも計算機シミュレーションによって検証が行われたようですが,そのシミュレーションと本記事のシミュレーションとは全く関係ありません.私が雑念として考えたモデルなので,アルゴリズムが正しいかどうかも分かりません.尚,実際の働きアリの原理は上記よりももう少し多層的で,各アリに固有の閾値があり,これによって働き始めるタイミングが決まるようです.が,そういったことを全て実装すると大変なので,話をシンプルにして実装します.

以下のような条件を設定して,時間を進め,アリの集団で何が起こるかを観察します.

今回題材とするのは,上記研究でも取り上げていた,アリの卵掃除です.アリの卵は,非常に汚れに弱く,汚れが蓄積していくと卵自体が死んでしまいます.そこで,働きアリの唾液で卵を嘗め回し,常に卵を清潔に保つというのは,働きアリの重要な仕事のひとつです.

このタスクをアリの集団が行う状況は,以下の点を実装すれば良いでしょう.

①卵には,汚れた卵と汚れていない卵がある.
②アリには働きアリと働かないアリがいる
③卵は時間が経過するほど汚れている時間が蓄積していく
④汚れの蓄積時間がある一定を超えると卵は死んでしまう
⑤汚れの蓄積時間がある一定を超えると他の汚れていない卵も汚してしまう
⑥死んだ卵も汚れが蓄積していく
⑦働きアリは生きている卵のみ掃除する
⑧働きアリの「掃除」は,卵の汚れの蓄積時間を減らす働きをする
⑨働きアリの掃除のスピードはランダムに速かったり遅かったりする
⑩働きアリは掃除すればするほど疲れがたまっていく
⑪疲れきると働きアリは,全く働かなくなり,二度と復活しない(これは示唆に富んだ設定です.)
⑫疲れきった働きアリが出ると,入れ替わりに働かないアリが働き始める

これらを実装します.MATLABで実装すると以下のようになります.構造体やMATLAB特有の手法を使うことでもっとシンプルになる気がしますが,MATLABの構造体の定義が少し面倒なので,あまりそういった手法は使っていません.その代わり,配列の定義以外のところは,他の言語でも楽に移植できるコードになっていると思います.(移植する必要性が全くありませんが)

clear;
close all;
numegg=100;%卵の数
egglife=1;%卵の寿命(汚れが蓄積すると死ぬ)
infecttime=1;%汚れ状態が伝染する時間
numant=100;%アリの数
workant=100;%働きアリ
tired=5;%アリが疲れて働かなくなる時間
egg_islife=ones(numegg,1);%卵が生きているかどうか(1:生 0:死)
egg_isdirty=zeros(numegg,1);%卵が汚れ状態かどうか (1:汚れ状態 0:汚れ状態でない)
egg_isdirty(1:30)=1;%汚れ状態は初期状態では30個
egg_dirtytime=zeros(numegg,1);%卵の汚れの蓄積時間
ant_life=ones(numant,1)*tired;%アリの疲労度(0になると疲れて働かなくなる)
ant_iswork=zeros(numant,1);%アリが働いているか(1:働いている 0:働いていない)
ant_lovework=zeros(numant,1);%働きアリかどうか(1:働きアリ 0:働かないアリ)
for k=1:workant,
    ant_lovework(k)=1;
    ant_iswork(k)=1;
end
delta=0.01;%刻み時間
tmax=4;
len=tmax/delta;
t=zeros(len,1);
answer=zeros(len,1);%汚れている卵の数
deadegg=zeros(len,1);%死んだ卵の数
liveant=zeros(len,1);%疲れていないアリの数(働きアリ・働かないアリ両方含む)
for k=1:len,
    deltat=tmax*delta;
    t(k)=k*delta;
    for n=1:numegg,
        if egg_isdirty(n)==1
            %汚れが蓄積する
            egg_dirtytime(n)=egg_dirtytime(n)+deltat;
            %汚れの蓄積時間が閾値を超えると,汚れていない卵に汚れが伝染する
            if egg_dirtytime(n)>=infecttime
                for m=1:numegg,
                    if egg_isdirty(m)==0
                        egg_isdirty(m)=1;
                        break;
                    end
                end
            end
       %汚れの蓄積時間が閾値を超えると卵は死ぬ
            if egg_dirtytime(n)>=egglife
                egg_islife(n)=0;
            end
        end
    end
    for n=1:numant,
        if ant_iswork(n)==1&&ant_life(n)>0
            ant_life(n)=ant_life(n)-deltat;
            for m=1:numegg,
                if (egg_isdirty(m)==1)&&(egg_islife(m)==1)
           %働きアリが卵を掃除する(掃除するスピードはランダムである)
                    egg_dirtytime(m)=egg_dirtytime(m)-deltat*rand(1,1)*0.1;
                    if egg_dirtytime(m)<=0
                        egg_isdirty(m)=0;
                        egg_dirtytime(m)=0;
                    end
                    break;
                end
            end
        end
    %疲れきった働きアリと入れ替わりに働かないアリが働き始める
        if ant_life(n)<0&&ant_lovework(n)==1
            for m=1:numant,
                if ant_life(m)>=0&&ant_lovework(m)~=1
                    ant_iswork(m)=1;
                end
            end
        end
    end
    for n=1:numegg,
        answer(k)=answer(k)+egg_isdirty(n);
        deadegg(k)=deadegg(k)+(1-egg_islife(n));
    end
    for n=1:numant,
        if ant_life(n)>0
            liveant(k)=liveant(k)+1;
        end
    end
end
hold on;
plot(t,answer);
plot(t,deadegg,'r');
plot(t,liveant,'g');
xlim([0 4]);
 legend('汚れている卵','死んだ卵','疲れていないアリ');

このシミュレーションの結果が以下のようになります.今回は,卵の数100個,そのうち初期状態で汚れている卵を30個,アリ全体を100匹でシミュレーションを行いました.

まず,全員が働きアリの場合.


横軸は時間tです.青が汚れた卵の数.赤が死んだ卵の数,緑がまだ働けるアリの数です.ここで注目したいのは,赤い線が,上限の100に達するタイミング.t=1.6あたりで,全ての卵が死んでしまいました.

次は,働かないアリを20匹にした場合.

卵が全滅するまでの時間がt=3まで延びました.雑なモデルと実装にしてはよくこんなに明白に差が出たもんだと自画自賛したいですね.

いかがでしょうか?この結果を見れば,働かないアリの必要性は議論の余地がないですね.働かないアリはいざというときの最後の砦として,自分の力を温存しているのです.そしてそれは,集団をより長く維持するのに貢献しています.

これは,人間社会でも同じではないでしょうか?今働いている人が疲れきってしまったとき,今働いている人だけでは手に負えない課題が発生した時.そんな時こそ働いていない人の出番です.人類が直面する未曾有の危機に,立ち向かう.そんな人間になりたいものです.

・・・というのは全て嘘です.とにかく働きたくありません.明日も客先の現場で仕事なので朝早く起きなければならず,とても憂鬱です.と書いて気が付きましたが,私がそんなにも仕事をしたくない理由は,朝決まった時間に起きなければならないからかもしれません.

2017年6月8日木曜日

boidの色々な動き 〜自分だけの群衆を作り出せ!〜

前回の記事でboidという,初歩的な群衆シミュレーションのアルゴリズムを紹介しました.bloggerを使っている方はご存知だと思いますが,bloggerの管理画面で,各記事が何回閲覧されたかというのが表示されます.その数が1日2日でこれまでの記事に匹敵するくらいの数になっていて驚いています.Twitterで組み込みやプログラミング分野に非常に影響力のある方に紹介していただいたお陰かもしれません.ありがとうございます.

せっかくですので,もう何回か,boidについて突っ込んだ記事を書いてみたいと思います.(前回の前半で,私のしようもないディズニー話がなければ,一つの記事で済んでいたかもしれませんね.)

前回の記事の続きという形で書きますので,前回の記述内容等を当記事で詳しく引用することは行いません.

1, 条件②の定数について

boidアルゴリズムにおける条件②について,多くのサンプルにはなかった定数での割り算を行っていたのを覚えていますか?これは,そうしないアルゴリズムでは,各個体の動きが不自然になるからでした.

②の部分を以下のように変えます.

vx2-=(x[n]-x[k])/2;
vy2-=(y[n]-y[k])/2;

するとこんな感じ.上が前回のもの.下が今回のものです.



どうでしょう?群れてはいますけどなんか変な感じは否めません.(個人の見解です)

ここで,各条件の計算式はそのままに,個体数(len)を1500,半径(r)を2,近さの閾値(dist)を10に設定すると,このようになります.



いいですね.昔スーパーファミコンに「スーパードンキーコング 謎のクレミス島」というゲームがあって,虫の大群が追いかけてくるというステージがあったのですが,それを思い出す質感です.

2, 速度算出時の重み係数について

各条件の計算結果を結合して速度を算出する時に,重み係数をかけていたのを覚えていますか?この係数の値を変えることで,動きの印象が大きく変わります.今回はProcessingのマウス位置検出機能を使って,重み係数を毎フレームで以下のように計算します.

float r1 = (float)mouseX/Width;
float r2 = abs((float)mouseX-(float)mouseY)/Width;
float r3 = (float)mouseY/Height;

マウス位置で群れの動きに変化が生まれます.


2017年6月4日日曜日

Processingでboid 〜初歩的な群衆シミュレーション〜

実は私,マニアとまでは行かないまでも,ディズニーは結構好きです.東京ディズニーランドなんかは,子供の頃は毎年長期休みのたびに家族で行っていましたし,大人になってからも何度か一人で行った程度には好きです.
ディズニーのアニメーション作品も,子供の頃から結構よく見ていました.親が結構ディズニー好きだったので,普段ビデオソフトは絶対買わずに,テレビ放送時に録画してすませるのに,ディズニーのアニメーションだけはビデオを買ってきて見せられたものでした.見終わった後で母親に感想を聞かれ,その時に母親が望む内容の感想を答えられないと(どのシーンを見てこのように思った.この物語の教訓を受けてこれから○○して生きていこうと思った.というような感想)急に不機嫌になり,育て方を間違えた的なことを言われるので,私は家族でディズニーのビデオを見るのがとても嫌でした.その代わり,家で留守番している時に一人で見るのは好きでした.
私の地元では,子供が見るアニメーション映画といえばほとんどディズニーくらいでした.後はドラえもんやクレヨンしんちゃんといった,テレビアニメーションの劇場版くらいで,スタジオジブリなどは,もののけ姫の頃までほとんど流行っていなかったと記憶しています.なので子供は皆当然のように公開されるディズニー映画を見ていました.
その頃公開されていた作品群は,ディズニーの長編アニメーション史で言えば,「ディズニー・ルネッサンス」と言われる時代の作品群になります.「白雪姫」から始まり,「ふしぎの国のアリス」や「ピーター・パン」等のヒット作を連発していたディズニー社は,創業者ウォルト・ディズニーがディズニーランドの建設・運営事業に熱中し始めた頃から徐々に陰りが出始めます.更に,ウォルト・ディズニーの死後,実写ハリウッド映画の台頭に押されディズニー社の人気は風前の灯火となりました.その状況を打破したのが1989年公開の「リトル・マーメイド」という作品で,これ以後の1990年代,爆発的なヒット作を飛ばし続け,現在に続く不動の人気を築いていきました.この時期が「ディズニー・ルネッサンス」で,私や私の幼馴染達がディズニー映画に親しんでいた時期です.
この時期の作品の特徴は,段階的にピクサー社の3DCGによる表現が取り入れられたことです.「リトル・マーメイド」で初めてピクサー社のプログラムを取り入れた(この時使ったのは3DCGではなく,セル画の彩色支援ソフトウェア)ディズニーは,次々とコンピュータによる表現を作品に取り入れるようになります.


©︎Disney

有名なのが「美女と野獣」の舞踏会のシーン.ダンスホールの背景を3DCGで描き,ダイナミックに視点移動することで,これまでのアニメーションでは出すのが難しかった奥行きのある映像表現を行い,世界中のアニメータを驚かせました.また「アラジン」では,飛ぶだけでなく,感情を持って動く魔法の絨毯の造形の全てが3DCGによって行われました.




©︎Disney

このように,ディズニー社を不動の地位にまで押し上げた作品群において,3DCGが大きな役割を果たしてきたと言えるでしょう.そして,この時期の3DCGを取り入れた表現において最も実験的だったと思われるのが1996年「ノートルダムの鐘」におけるラスト部分のこのシーン.

©︎Disney

この作品はご存知ヴィクトル・ユゴーの小説を原作とした作品ですが,小説とは違ってハッピーエンドとなっており,(少しネタバレ)悪者は死んで,奇形児の孤児として差別や偏見に晒されてきた主人公はついに町の人々に受け入れられます.町に平和が訪れたことを喜びながら人々がパリの街を行進し,やがて視点がだんだん遠のいていき,街全体が映ったところで物語が終了します.(このシーンだけ見たい方は,ぶっちゃけYouTubeにありますのでご覧ください.

この壮大なシーンは,群衆の動きを計算モデルによってシミュレーションし,自動生成することによって作られました.大勢の人々の動きを一人一人描き分け,一つのシーンとして纏め上げる,手書きであれば膨大な作業量が必要とされますが,それをコンピュータによる自動生成で行ってしまう,これまでの作品の奥行きやダイナミズムの表現とは違った角度のCGの活用法に,人々は改めて驚いたのでした.

群衆の動きをシミュレーションするアルゴリズムは,この作品の10年ほど前にCraig Reynoldsによって提唱されたBoidというアルゴリズムが有名です.Boidは人ではなく,鳥の群れ(っぽい動き)をシミュレーションするアルゴリズムですが,非常に自然な動きをするため,改良されて多くのアニメーションで使われています.言わば群衆シミュレーションの古典とも言えるアルゴリズム.現在では,各言語でクラスとしてそのまま使えるコードがあったりして,わざわざ実装する気にもならないようなアルゴリズムですが,改めて実際どうなっているのか復習したいと思います.

言語はProcessing.再利用しやすい実装を目指して色々と活用したいわけではなく ,どういう処理手順で実現されているか知りたいので,draw()関数にダダダダダッとアルゴリズムを書いていきます.

まずは初期化です.個体をランダムな位置に,ランダムな初速度で配置します.色も適当につけておきます.


int len=50; //個体数
int Width=960,Height=540; //解像度
float[] x=new float[len]; //個体のx座標
float[] y=new float[len]; //個体のy座標
float[] r=new float[len]; //個体の半径
float[] dx=new float[len]; //個体のx方向の移動量
float[] dy=new float[len]; //個体のy方向の移動量
// 色
float[] red=new float[len];
float[] green=new float[len];
float[] blue=new float[len];
void setup(){
    int k;
    for(k=0;k<len;k++){
        r[k]=10;
        x[k]=random(r[k]/2,Width-r[k]/2);
        y[k]=random(r[k]/2,Height-r[k]/2);
        red[k]=random(0,255);
        green[k]=random(0,255);
        blue[k]=random(0,255);
        dx[k]=random(-2,2);
        dy[k]=random(-2,2);
    }
    size(960,540);
    background(0,0,0);
}




 

いよいよdraw()関数の中にboidを書いていきます.boidは,各個体がランダム(ランダムな位置や速さ)に動いている時に,以下の条件を加えると実現できます.

① それぞれの個体が,全体の分布の中心に向かって進む
② 互いに近づきすぎたら反対へ動く
③ 他の個体の速度に合わせる

現在の各個体の位置や速度から,①〜③についての移動速度を別々に算出し,最後にそれを合計して,その個体の移動速度とします.

 float vx1=0,vy1=0,vx2=0,vy2=0,vx3=0,vy3=0; //それぞれの条件の速度
 float ax=0,ay=0; //位置の平均値
 float avx=0;avy=0; //速度の平均値

まずは①.自分以外の位置の平均を求め,自分自身が平均に向かう方向の速度を求めます.
        //① それぞれの個体が,全体の分布の中心に向かって進む
        //自分以外の個体の位置の平均を求める
        for(n=0;n<len;n++){
            if(k!=n){
                ax+=x[n];
                ay+=y[n];
            }
        }
        ax/=len-1;
        ay/=len-1;
        //自分を平均値に近づける方向を速度とする(適当に重みをつけて)
        vx1=(ax-x[k])/400;
        vy1=(ay-y[k])/400;
            }
        }
ここで最後の/400の400は定数です.この値を変えると,動きの印象が変わってきます.

次は②.自分以外の個体との距離を求め,近づいていたら,反対方向の速度を加えます.

Post to 
LiveJournal
        //② 互いに近づきすぎたら反対へ動く
        for(n=0;n<len;n++){
            if(k!=n){
                float dist=sqrt((x[k]-x[n])*(x[k]-x[n])+(y[k]-y[n])*(y[k]-y[n])); //自分とそれ以外の距離
                //距離が近かったら反対方向の速度成分を加える
                if(dist<30){
                    vx2-=(x[n]-x[k])/25;
                    vy2-=(y[n]-y[k])/25;
                }
            }
        }
<30 pre="" vx2-="(x[n]-x[k])/25;" vy2-="(y[n]-y[k])/25;"> ここでも速度を定数で割り算しています.ネットで見かけるサンプルでは,割り算していないものが多いのですが,割り算しないと,見えない壁にぶつかって跳ね返ったような動きが目立つのですが,割り算するとより自然な動きになります.

最後に③.今度は速度の平均を求め,自分自身の速度との差に応じて,速度を算出します.


        //③ 他の個体の速度に合わせる
        //自分以外の個体の現在の速度の平均を求める
        for(n=0;n<len;n++){
            if(k!=n){
                avx+=dx[n];
                avy+=dy[n];
            }
        }
        avx/=len-1;
        avy/=len-1;
        //自分の速度との差を速度成分に加える
        vx3=(avx-dx[k])/2;
        vy3=(avy-dy[k])/2;


これで①〜③の計算ができました.これらの値を現在の速度に加えていきます.この時,①〜③の値それぞれに異なった重みをつけることで,動きをより自然にします.この値も,変化させると印象が変わります.

また,速度が想像以上に早かった場合,速度の正規化を行います.

     //計算結果を各個体の速度に加える(重み付き)
        dx[k]+=vx1*r1+vx2*r2+vx3*r3;
        dy[k]+=vy1*r1+vy2*r2+vy3*r3;
        
        //速度の絶対値が4以上だったら正規化
        float velocity=sqrt(dx[k]*dx[k]+dy[k]*dy[k]);
        if(velocity>4){
            dx[k]=(dx[k]/velocity)*4;
            dy[k]=(dy[k]/velocity)*4;
        }
        //自分の速度との差を速度成分に加える
        vx3=(avx-dx[k])/2;
        vy3=(avy-dy[k])/2;

以下がソースの全体です.

int len=50; //個体数
int Width=960,Height=540; //解像度
float[] x=new float[len]; //個体のx座標
float[] y=new float[len]; //個体のy座標
float[] r=new float[len]; //個体の半径
float[] dx=new float[len]; //個体のx方向の移動量
float[] dy=new float[len]; //個体のy方向の移動量
// 色
float[] red=new float[len];
float[] green=new float[len];
float[] blue=new float[len];
void setup()
{
    int k;
    for(k=0;k<len;k++){
        r[k]=10;
        x[k]=random(r[k]/2,Width-r[k]/2);
        y[k]=random(r[k]/2,Height-r[k]/2);
        red[k]=random(0,255);
        green[k]=random(0,255);
        blue[k]=random(0,255);
        dx[k]=random(-2,2);
        dy[k]=random(-2,2);
    }
    size(960,540);
    background(0,0,0);
}

void draw()
{
    int k;
    int n;
    //速度を足し合わせる重み係数
    float r1 = 1;
    float r2 = 0.8;
    float r3 = 0.1;
    background(0,0,0);
    for(k=0;k<len;k++){
        float vx1=0,vy1=0,vx2=0,vy2=0,vx3=0,vy3=0; //それぞれの条件の速度
        float ax=0,ay=0; //位置の平均値
        float avx=0;avy=0; //速度の平均値
        noStroke();
        fill(red[k],green[k],blue[k]);
        ellipse(x[k],y[k],r[k],r[k]);
        
        //① それぞれの個体が,全体の分布の中心に向かって進む
        //自分以外の個体の位置の平均を求める
        for(n=0;n<len;n++){
            if(k!=n){
                ax+=x[n];
                ay+=y[n];
            }
        }
        ax/=len-1;
        ay/=len-1;
        //自分を平均値に近づける方向を速度とする(適当に重みをつけて)
        vx1=(ax-x[k])/400;
        vy1=(ay-y[k])/400;
        
        //② 互いに近づきすぎたら反対へ動く
        for(n=0;n<len;n++){
            if(k!=n){
                float dist=sqrt((x[k]-x[n])*(x[k]-x[n])+(y[k]-y[n])*(y[k]-y[n])); //自分とそれ以外の距離
                //距離が近かったら反対方向の速度成分を加える
                if(dist<30){
                    vx2-=(x[n]-x[k])/25;
                    vy2-=(y[n]-y[k])/25;
                }
            }
        }
        
        //③ 他の個体の速度に合わせる
        //自分以外の個体の現在の速度の平均を求める
        for(n=0;n<len;n++){
            if(k!=n){
                avx+=dx[n];
                avy+=dy[n];
            }
        }
        avx/=len-1;
        avy/=len-1;
        //自分の速度との差を速度成分に加える
        vx3=(avx-dx[k])/2;
        vy3=(avy-dy[k])/2;
        
        //計算結果を各個体の速度に加える(重み付き)
        dx[k]+=vx1*r1+vx2*r2+vx3*r3;
        dy[k]+=vy1*r1+vy2*r2+vy3*r3;
        
        //速度の絶対値が4以上だったら正規化
        float velocity=sqrt(dx[k]*dx[k]+dy[k]*dy[k]);
        if(velocity>4){
            dx[k]=(dx[k]/velocity)*4;
            dy[k]=(dy[k]/velocity)*4;
        }
        
        //求めた速度で移動
        x[k]+=dx[k];
        y[k]+=dy[k];
        
        //端まで行ったら跳ね返る
        if(x[k]-r[k]<=0){
            x[k]=r[k];
            dx[k]*=-1;
        }
        if(y[k]-r[k]<=0){
            y[k]=r[k];
            dy[k]*=-1;
        }
        if(x[k]+r[k]>=Width){
            x[k]=Width-r[k];
            dx[k]*=-1;
        }
        if(y[k]+r[k]>=Height){
            y[k]=Height-r[k];
            dy[k]*=-1;
        }
    }
}

こうして作ったboidの動きがこちらです.
単純な割に,色々変えると面白い動きになるので,遊びがいがありそうです.






2017年5月3日水曜日

Audio Sculptで音を彫刻

ずっとクレジットカードを持たない人生を送ってきたのですが,iPhoneを機種変更した時に,ソフトバンクカードというプリペイド式のVISAカードとしても使えるTポイントカードを貰ったので,クレジットカードでないと決済できない商品も買えるようになりました.

というわけで,ずっと憧れだった音響分析・合成ソフトウェア "Audio Sculpt" を買いました.

このソフトを販売しているのは,フランスのIRCAMという研究所で,日本語ではフランス国立音楽音響研究所と呼ばれる組織です.ここでは,電子音楽やメディアアート等の作品制作や,それに付随する技術の研究開発が行われていて,フランスの現代音楽の拠点とも言える施設です.電子音楽やメディアアートの世界で有名なソフトウェアMaxがIRCAMで作られたことはあまりにも有名です.
IRCAMを拠点として,フランスでは「スペクトル楽派」と呼ばれる現代音楽の潮流ができました,スペクトル楽派とは,音響学や音響信号処理の知見を音楽として再解釈して楽曲を構築する流派で,創始者の一人と言われるトリスタン・ミュライユという作曲家が有名です.(スペクトル楽派やその音楽は知らなくても,トリスタン・ミュライユという作曲家の存在を知っている人は多いと思います.)

スペクトル楽派では,音響の知識を楽曲として昇華させるので,作曲の際には音についての綿密な分析,解釈,合成が不可欠です,そういった事を音楽家が行うためのソフトウェアがAudio Sculptです.Sculptとは彫刻という意味.音をまさに彫刻のように加工する事ができるソフトウェアです.
買ったばかりなのでまだ簡単な事しかわかりませんが,あまりwebで紹介している例がないので,簡単な使い勝手を紹介します.

Audio Sculptは,IRCAM Forumというサイトで販売されています.ここでは,IRCAMの研究開発の成果がソフトウェアやライブラリとして有料・無料で配布されています.Audio Sculptは有料ソフトウェアで,€180.00です.


Audio Sculptで音声ファイルを開くとこのような画面になります.


この段階では再生ができるだけです.▶︎ボタンを押せば音声が再生されます.


Analysis -> Sonogram Analysis

を選んで,分析の各パラメータを設定して,"Do Analysis"をクリックすると,






波形の下にスペクトログラムが表示されます.Audio Sculptでは,このスペクトログラムをまさに彫刻のように加工して,音を再合成する事ができます.

例えば,"Selectio  Tool"モード(画面左上のボタン群の二番目)を選択して,スペクトログラム内のある範囲を選択し,Edit -> Cut を行うと,選択範囲をごっそり消す事ができます,


ここで,SVPと書かれた再生ボタンを押すと,今行った編集が適用された音声を聞く事ができます.ただし,ここで聞けるのは,あくまでもプレビュー的なものなので,上段の波形や,画面右のスペクトル表示はオリジナルのままです.

こんな風にフリーハンドで範囲を選択する事もできます.


また,範囲を選択した後,"Arrow Tool"モードを選択する事で,選択範囲を他の場所に持っていく事もできます.勿論,この操作も音として反映されます.



更に,"Filter Tool"モードにして,範囲選択をし,ctrlキーを押しながらマウスを上下させると,選択範囲を強調・抑制ができます.

グレーが強調,赤が抑制した部分です.




ライン上にフィルタを定義する事もできます.同様に強調,抑制が可能です.


一通りいじりましたね.


一通り編集が終わったら,Processing -> Process Treatments を選択し,各パラメータを設定し,Processを押すと,今編集した内容が適用されたオーディオファイルが書き出されます.




いかがでしょうか?今やった内容が適用された音声になったのではないかと思います.

これが,Audio Sculptのオーソドックスな使い方です.これらの操作自体は,一つ一つはこのソフトを使わなくても簡単にできる内容のものです.しかし,これほどインタフェースとしてまとめられているソフトウェアは他にはないと思います.しかも,これはAudio Sculptの機能のほんの一部でしかありません.

また面白い機能を見つけたら,本ブログで報告したいと思います.

2017年4月25日火曜日

機械学習を学習してみるか 〜単純パーセプトロンで顔文字の感情推定〜

最近やたら流行ってますね.人工知能.これからの世界は人工知能が牛耳って,人間の大半は人工知能に代替され職を失うようです.そうしたら働かなくて済むので1日も早くそうなって欲しいですが,本当に情報科学のあらゆる分野で「人工知能」というワードが聞かれるようになりました.どうやら,情報科学に携わるものは,人工知能が分からなければ,この先最新の研究内容に触れることはできなさそうです.
というわけで,ずっと避けてきた,人工知能分野の勉強を少しずつ始めてみたいと思います.人工知能といってもいろいろあると思いますが,まずは基本,機械学習から.

ニューラルネットワークの基本,単純パーセプロトロンをやってみます.ちなみに,この記事は,単純パーセプトロンとは何か,どういう仕組みなのかを勉強し,それを記事として記録するものではありません.そういった学習が終わった後,私なりに「こういう風に使えばいいのか?」というのをやってみるという内容です.ですから,機械学習初心者のみならず,ベテランの方もお読みいただいて,「ここはこうすればいいのだ」等の指摘をいただけると非常に励みになります.

単純パーセプトロンというのは皆さん知っていますね.こういうやつです.



青グループと赤グループを教師データとして学習しておけば,両グループの境界線を勝手に求めてくれます.一度学習してしまえば,次に新しいデータが入力されてもどちらのグループに属するのか判定することができます.

ここまでは,本を読むかGoogleで検索すれば簡単に概要およびプログラムを見つけることができます.ここからが本題.折角覚えた単純パーセプトロン.何かしらに使ってみたくなるのが人情です.そこで,今回はこんなことをやってみることにしました.今回の言語はPythonを使いました.

「いくつかの顔文字とその感情(今回は 喜と哀)を学習し分類する.その後未知の顔文字の感情を推定する」

今回は喜と哀を分析します.どうでもいいことですが,人間の基本的な感情としてよく言われる「喜怒哀楽」は単純パーセプトロンで推定することはできません.何故か?分からなくてなおかつ気になる方はもう少し単純パーセプトロンについて調べてみましょう.

まず,学習データとなる顔文字を用意します.ここでは,「顔文字」とは丸括弧の中に3文字の記号が入ったタイプを想定します.今回はこれらを使います.

(^_^)(^∀^)(^o^)(^ω^)(;_;)(>_<)(◞‸◟)(゚Д゚)

左から

喜 喜 喜 喜 哀 哀 哀 哀 

とします.今考えると,一番右は哀じゃない気がしますが,例題なので大目に見ます .

次に,それぞれの顔文字に使われている文字に数字を割り当てます.とりあえず,喜に使われそうな文字は小さい数字,哀に使われそうな文字は大きい数字にしました.こんな感じで0から12まで割り当てました.

(^  _ ^)(^∀^)(^o^)(^ω^)( ;  _;) (> _ <) (◞  ‸  ◟)  (゚ Д゚)
 0 7       1     6      2   12     10 11   3 9 8    5 4

ここまで決めたら,データとしても顔文字を作ります.また,喜を1,哀を-1として教師データも作ります,

N=8
ETA=0.1
#顔データ
face=[[0,7,0],[0,1,0],[0,6,0],[0,2,0],[12,7,12],[10,7,11],[8,9,3],[5,4,5]]
#教師データ
t=[1,1,1,1,-1,-1,-1,-1]

今回は,顔データとして3次元のベクトルを8個定義しましたので,それらを二分する境界面(3次元なので面になる)

ax+by+cz+d=0

を単純パーセプトロンによって求めます.機械学習なんて大層な名前の技術の割に,結構簡単にできます,詳しいメカニズムは書籍や他のブログ等を当たってみてください.

w=[1.0,1.0,1.0,1.0] #ax+by+cz+d=0の各係数を[a,b,c,d]とする
correct=0
while correct<N:
    correct=0
    for i in range(N):
        if np.dot(w,[face[i][0],face[i][1],face[i][2],1])*t[i]>0:
            correct+=1
        else:
            w+=ETA*np.array([face[i][0],face[i][1],face[i][2],1])*t[i]w=[1.0,1.0,1.0,1.0]

それによって求めた境界面の式がこちら.

-1.0x+0.1y+-0.5z+0.9=0

そして,境界面と上記顔文字(ベクトル)との距離を求めると以下のようになりました.

(^_^): 1.4253932902
(^^): 0.890870806375
(^o^): 1.33630620956
(^ω^): 0.979957887012
(;_;): -14.6102812245
(>_<): -12.3831042086
(◞‸◟): -6.85970520909
(゚Д゚): -5.52339899952

ちゃんと喜と哀がプラス側とマイナス側に分かれています.喜側の距離が小さい気がしますが、喜側は目の部分が「^」しか使っていないため,特徴の学習が甘くなっているのではないかと考えられます.

それでは,未知の顔文字の感情推定.

上記文字をランダムに使って顔文字を生成.それぞれの感情を推定します.数値が正なら喜,負なら哀です.

(<‸;): -13.5412362569
(∀o^): 0.445435403187
(◟◟Д): -3.38530906422
(゚∀◟): -4.89978943506
(>Д;): -13.0958008537
(゚o>): -7.57240185419
(>Д^): -7.75057601546
(Д◟゚): -4.72161527379
(^ω∀): 0.534522483825
(<◟ω): -9.62140470885
(^Д∀): 0.7126966451

なるほど.納得できるようなできないような結果.まあこんな感じでしょう.

今回,学習によって求めた境界面と顔文字(ベクトル)との距離を求めましたが,距離が小さいということは,喜と哀の感情が曖昧,距離が大きければはっきりと感情が出ていると考えられるのでしょう.今回,(多分)上記理由で,喜側の距離があまり大きくなりませんでしたが,顔文字の定義の仕方やデータの作り方を工夫すればもっとはっきりと分かる結果になると思います.