可視化と可聴化のためのMathematica
〜Wolfram言語による映像と音楽の融合〜

面白法人カヤック  岩淵 勇樹

Join the Conversation @WolframJapan

はじめに

Mathematicaは ListPlot, ArrayPlot, Graphics, ListPlay, Exportなど、画像処理や音声処理の関数が充実

それに加えて数式処理の強みを活かし、画像⇔音声などの変換にも強い

本講演では5つの事例をプログラムを通して紹介

お品書き

物智数

bion

数式曲線

解析信号

フラクタル音楽

自己紹介

岩淵 勇樹 (IWABUCHI Yu(u)ki)

金沢大学自然科学研究科修了

博士（工学）

面白法人カヤック HTMLファイ部・人事部

博士論文は「図形と音声の変換方法とその応用に関する研究」

数システム → 画像
（その1）

Mathematicaによる可視化

数学的対象をプロットしたりグラフィックス表示するのが簡単

研究の効率を上げることができる

物智数

2進数や10進数などの「記数法」のひとつ

発表者が2006年に発見

Hyperlink["物智数（Butchi number, 平面的2進数）","http://butchi.jp/documents/planar_binary/"]

岩淵勇樹 卒業論文「2進数の平面的表現に関する考察」 (2007)

2進数変換の再発明

binaryFromPolynomial[n_]:=Reverse[CoefficientList[FixedPoint[Expand[((#-PolynomialMod[#,2])/2) p+PolynomialMod[#,2]]&,n],p]]

binaryFromPolynomial[5]

binaryFromPolynomial[(p^2)+1]

多項式を用いた実装

2進数の「2」を「p」に置き換えて多項式表現（ex: 5 = (Subscript[101, 2])= (2^2)+1 = (p^2)+1）

数式処理ソフトならではの実装

2進数の平面拡張

↓

「左に繰り上がる」から「左と上に繰り上がる」にルール変更

ビットが平面的に広がる

物智数の実装

butchiNumberFromPolynomial[n_?PolynomialQ]:=butchiNumber[CoefficientList[FixedPoint[Expand[((#-PolynomialMod[#,2])/2) (p+q)+PolynomialMod[#,2]]&,n],{q,p}]]

Format[butchiNumber[li_]]:=Grid[Reverse[li,{1,2}],Dividers->-1->True]

butchiNumberFromPolynomial[4]
First[%]

Table[butchiNumberFromPolynomial[n],{n,10}]

物智数の実装

1→黒、0→白としてプロット

butchiNumberPlot[butchiNumber[li_]]:=ArrayPlot[Reverse[li,{1,2}]]

Table[butchiNumberPlot[butchiNumberFromPolynomial[n]],{n,10}]

大きな物智数

Manipulate[butchiNumberPlot[butchiNumberFromPolynomial[n]],{{n,10},1,1000},SaveDefinitions->True]

セルオートマトンを用いた実装

butchiNumberCellularAutomaton[n_]:=butchiNumber[FixedPoint[CellularAutomaton[{Mod[#[[1]],2]+Quotient[#[[2]],2]+Quotient[#[[3]],2]&,{},{{0,0},{-1,0},{0,-1}}},{#,0},{{{2}}}]&,{{n}}]]

butchiNumberCellularAutomaton[5]

お蔵入りとなったが、CやJavaなどよりは短いコードで書ける

数式処理と豊富な関数により、実装の選択肢が増える

スライド

対話的なプログラミングにより、簡潔な関数定義と1行レベルの命令で実験可能

即座に可視化結果を確認できるので試行錯誤しやすい

数システム → 画像
（その2）

bion

音楽における「倍音」をイメージした数システム

螺旋状に1ずつ加算される

放射状に2倍される

bionの実装

bion[n_]:=Log2[n]*(E^(2π I Log[2,n]))

Graphics[Table[Text[k,ReIm[bion[k]]],{k,1,64}]]

複素数でグラフィックス処理

x座標を実部、y座標を虚部として扱う

リスト操作や幾何学的変換が{x, y}のリストに比べて簡便

個人的にはほとんどこの方式を用いている

ReIm[3+4I]

Complex@@{3,4}

(((-1))^((4/5)Range[0,5]));
ReIm/@%;
Graphics[Line[%]]

参考

Hyperlink["複素数の魔法 - Qiita","https://qiita.com/butchi_y/items/a4a7e992870d75e5ee7a"]

bionの作品化

研究の一環でありつつ、CG作品としても出力

Exportによるグラフィックス出力はラスター形式もベクトル形式もサポート

n=2048;
gBion=Graphics[Table[{Hue[0.57,(((1-Log[2,k]/Log[2,n]))^2),1],Disk[ReIm[bion[k]],0.5(k^(-0.25))]},{k,n,1,-1}]]

SetDirectory[NotebookDirectory[]];

Export["bion.eps",gBion]

waon

bionの関連作品

RGB, HSB, CMYKなど表色系も充実

waon[n_]:=(E^(2π I Log[2,n]))

n=256;
gWaon=Graphics[Table[{Hue[0.58+Log[2,k],1-(((Log[2,k]/Log[2,n]))^2),1],Disk[ReIm[waon[k]],0.3(k^(-0.2))]},{k,n,1,-2}]]

Export["waon.eps",gWaon]

映像 → 数式

数式曲線

Wolfram|Alphaの“Hatsune Miku curve”がネットで話題に

アニメーション版の『Bad Apple Curve!!』を制作 (2015)

モノクロのアニメーション版ができたので今度はカラー版も (2018)

参考

Hyperlink["Making Formulas… for Everything—From Pi to the Pink Panther to Sir Isaac Newton—Wolfram Blog","http://blog.wolfram.com/2013/05/17/making-formulas-for-everything-from-pi-to-the-pink-panther-to-sir-isaac-newton/"]

でびどる曲線！

Hyperlink["https://twitter.com/butchi_y/status/1072374226588459008"]

Special Thanks: コチンPa!製作委員会

フーリエ記述子

曲線の座標をフーリエ変換し、低周波成分のみを拾う

今回はP型ではなくG型のフーリエ記述子を採用

参考

岩淵勇樹・秋田純一 「非周期的な画素配置を持つCMOSイメージセンサの基礎検討」, 電子情報通信学会技術研究報告 (2007)

でびどる曲線！ 制作過程

実際使ったコードを整理し、一連の流れを抽出

img=

binaryImageAira[image_]:=ColorNegate[Binarize[ColorDistance[image,],0.5]]

binImg=binaryImageAira[img]

gcImg=GeodesicClosing[binImg,1]

mcImg=MorphologicalComponents[gcImg];
%//Colorize

layerLi=Image/@Table[Unitize[mcImg-n],{n,1,Max[mcImg]}]

selLayerLi=Select[layerLi,Total[Flatten[1-ImageData[#]]]>200&]

findBoundary[img_]:=Module[{imgData=ImageData[Binarize[ImagePad[ImagePad[img,-2],2,1]]],
d=(({{0, 1}, {1, 1}, {1, 0}, {1, -1}, {0, -1}, {-1, -1}, {-1, 0}, {-1, 1}})),
vec=0,
nextPointLi,offset,boundary},
boundary=Position[imgData,0,2,1];
While[Length[boundary]==1||First[boundary]≠Last[boundary],
nextPointLi=Select[RotateLeft[d,vec-3],Extract[imgData,Last[boundary]+#]==0&];
If[nextPointLi=={},Break[]];
offset=First[nextPointLi];
AppendTo[boundary,Last[boundary]+offset];
vec=First[First[Position[d,offset]]]-1;
];
Return[Most[Reverse/@boundary]];
]

lineLi=Cases[findBoundary/@selLayerLi,Except[{}]]
Graphics[Line/@lineLi]

fft[li_]:=Fourier[Complex@@#&/@li,FourierParameters->{-1,1}]
fftLi=fft/@lineLi

atan[0,b_]:=(π/2)
atan[0.,b_]:=N[(π/2)]
atan[a_,b_]:=ArcTan[(b/a)]+({{\[Piecewise], ({{0, a≥0}, {π, True}})}})
(* reSin・imSinの係数と位相だけのリスト *)
reSinComponent[c1_,c2_]:=With[{p=-Im[c1]+Im[c2],q=Re[c1]+Re[c2]},{(Sqrt[(p^2)+(q^2)]),atan[p,q]}]
imSinComponent[c1_,c2_]:=With[{p=Re[c1]-Re[c2],q=Im[c1]+Im[c2]},{(Sqrt[(p^2)+(q^2)]),atan[p,q]}]
(* 係数・位相リストのリスト *)
reComponents[components_]:=With[{n=((Length[components])/2)},Join[{Re[First[components]]},Table[reSinComponent[components[[k+1]],components[[-k]]],{k,1,n}]]]
imComponents[components_]:=With[{n=((Length[components])/2)},Join
[{Im[First[components]]},Table[imSinComponent[components[[k+1]],components[[-k]]],{k,1,n}]]]

reComp=reComponents/@fftLi;
imComp=imComponents/@fftLi;
Short[First[reComp]]
Short[First[imComp]]

sinFromComponent[{k_,θ_},n_]:=k Sin[n t+θ]
reExpressionList[n_]:=First[#]+Sum[sinFromComponent[PadRight[#,n,{{0,0}}][[m]],m-1],{m,2,n}]&/@reComp;
imExpressionList[n_]:=First[#]+Sum[sinFromComponent[PadRight[#,n,{{0,0}}][[m]],m-1],{m,2,n}]&/@imComp;

reExprLi=reExpressionList[10];
imExprLi=imExpressionList[10];
First[reExprLi]
First[imExprLi]

ParametricPlot[{reExprLi,360-imExprLi}\[Transpose],{t,0,2π},PlotRange->{{0,640},{0,360}},AspectRatio->Automatic]

Manipulate[ParametricPlot[Evaluate[{reExpressionList[n],360-imExpressionList[n]}\[Transpose]],{t,0,2π},PlotRange->{{0,640},{0,360}}],{{n,50},2,150,1},SaveDefinitions->True]

音声 → 曲線

解析信号

音声信号を複素信号化してプロット

周期信号に対して閉曲線が描かれる

解析信号フィルタ

負の周波数成分をカットするフィルタをFFT後のスペクトルに対して適用する

参考文献

岩淵勇樹，秋田純一，北川章夫 「閉曲線を利用した音色の操作法」 芸術科学会論文誌 (2011)

岩淵勇樹 博士論文 「図形と音声の変換方法とその応用に関する研究」 (2012)

解析信号の実演

step[n_Integer?Positive]:=1+Sign[Sin[π Range[1,-1,-2/n]//Most]]

step[10]

analyticSignalTransform[li_List]:=InverseFourier[Fourier[li]*step[Length[li]]]

snd=ExampleData[{"Sound","Flute"}]

sndData=snd[[1]][[1]]
sampleRate=snd[[1]][[2]]
len=Length[sndData[[1]]]

asLi=analyticSignalTransform[sndData[[1]]]

ListLinePlot[ReIm/@asLi,PlotRange->All,AspectRatio->Automatic]

asLi//Length

ListLinePlot[ReIm/@Take[asLi,{20000,20500}],PlotRange->All,AspectRatio->Automatic]

解析信号のアニメーション

fps=30;
totalSecond=N[len/sampleRate];
totalFrame=Floor[len*fps/sampleRate];
samplePerFrame=Floor[sampleRate/fps];
ListAnimate[Table[ListPlot[Take[{Re[asLi],Im[asLi]}\[Transpose],{frame*samplePerFrame+1,(frame+1)*samplePerFrame}],Joined->True,PlotRange->{{-1,1},{-1,1}},AspectRatio->Automatic,Axes->False,PlotStyle->Black],{frame,0,totalFrame-1}],fps]

さまざまな解析信号

楽器や声質など、音色に特有の形が現れる

参考

Hyperlink["解析信号ギャラリー", "http://butchi.jp/documents/analytic_signal/"]

図形 → 音楽

フラクタル音楽

フラクタル図形を可聴化するアート活動

TwitterとInstagramで約1年間毎日投稿、計459作品

Hyperlink["カーペットフラクタル(@carpet_fractal) • Instagram","https://www.instagram.com/carpet_fractal/"]

フラクタル音楽制作の流れ

フラクタル図形の描画

音声信号に変換

解析信号で可視化

参考

Hyperlink["フラクタル音楽 〜可視化と可聴化の世界〜","https://www.slideshare.net/yuukiiwabuchi9/ss-66759132"]

岩淵勇樹 「カーペットフラクタルの可聴化と音楽生成」, 情報処理学会第105回音楽情報科学研究会 (2014)

Hyperlink["http://butchi.jp/documents/mus105/"]

フラクタル生成から可聴化まで

iterate[x_]:=ArrayFlatten[x/.{0->(({{0, 0, 0}, {0, 0, 0}, {0, 0, 0}})),1->(({{1, 1, 1}, {1, 0, 1}, {1, 1, 1}}))}]

iterate[1]//Grid

iterate[iterate[1]];
Grid[%]

Nest[iterate,1,2];
Grid[%]

Nest[iterate,1,2]//ArrayPlot

Clear[n];
Manipulate[ArrayPlot[Nest[iterate,1,n]],{n,1,6,1},SaveDefinitions->True]

carpet=Nest[iterate,1,5];
ArrayPlot[carpet,Frame->False,PixelConstrained->1]

ListPlay[Flatten[carpet]]

randomIterate[x_]:=ArrayFlatten[x/.{0->(({{0, 0, 0}, {0, 0, 0}, {0, 0, 0}})),1->iterationTbl}]

SeedRandom[0]

iterationTbl=RandomInteger[{0,1},{3,3}];
Grid[iterationTbl]

randomCarpet=Nest[randomIterate,1,5];
ArrayPlot[randomCarpet,Frame->False,PixelConstrained->1]

ListPlay[Flatten[randomCarpet]]

shiftCarpet=BitXor[carpet,RotateLeft[carpet,{20,45}]];
ArrayPlot[shiftCarpet,PixelConstrained->1,Frame->False,ColorFunction->(Blend[{LCHColor[0.2,1,1.8],White},#]&)]

ListPlay[Flatten[shiftCarpet]]

まとめ

「数システム → 画像」「映像 → 数式」「音声 → 曲線」「図形 → 音楽」の変換について解説

数式を介して画像・音声間を柔軟に行き来可能

Wolfram Mathematicaは研究にも芸術活動にも最適

ご清聴ありがとうございました