ここ数年の僕はプログラミングというと音と統計ばかりやっているのですが、特に音の分析や加工に関してよく使うものはMatlab、Python、Rで実装したので、次はCommon LispかJuliaで実装してみようかと思います。最終的にはCommon Lispでなんでもやるのが夢なのですが、今回はJuliaで。
基本的なことがら
音ファイルを扱えないと、その後の分析も何もできませんので、いの一番にやることはJuliaでWAVファイルの読み込みです。音に関するパッケージをPkg Serverで探してみると、使えそうなパッケージはいくつかあります。その中でもLibSndFileを使えるようにしてくれるパッケージが(現時点で丸井綜研が使っている)Julia 0.6のテストを通過してるので良さそうです。
インストール
インストールは簡単で、Juliaのコマンドラインから
Pkg.add("LibSndFile")
とやるだけ。初回に一度やれば、次からはこの部分は不要です。
使用準備
using LibSndFile
でパッケージの読み込みができます。初回はプリコンパイルなどで少し時間がかかるかもしれません。
音の読み込み
読み込みにはload
を使います。JuliaのコマンドラインREPLではいろいろな情報を表示してくれます。以下ではCDから抜粋した15秒ほどのWAVファイルを読み込んでいますが、1チャンネルあたり587,264サンプルあり、それが2チャンネルあること、全体で約13.3秒、サンプリング周波数は44,100Hzなどの情報に加えて、波形の絶対値?をプリティプリントまでしてくれます。
julia> x = load("gershwin_piano_concerto.wav")
587264-frame, 2-channel SampleBuf{FixedPointNumbers.Fixed{Int16,15}, 2}
13.316643990929705s sampled at 44100.0Hz
▅▆▆▆▆▇▆▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▆▆▇▇▇▇▇▇▇▇▇▇▇▇▇▆▆▇▆▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▆▇▇▇▇▇▇▇▇▇▇▇▆▇▆▆▇▆▆▆▆▅▄
▅▆▆▆▆▆▆▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▆▇▇▆▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▆▆▆▆▆▆▆▅▄
IJuliaノートブック上では、なんとグラフィカルな波形表示と再生機能までしてくれます。(下はそのスクリーンショットです)
格納先の変数はSampledSignals.SampleBuf型で、data
とsamplerate
の二つのフィールドを持っています。
julia> x.data # 各サンプル値(2chファイルなので2列です)
587264×2 Array{FixedPointNumbers.Fixed{Int16,15},2}:
-3.0e-5Q0f15 0.0Q0f15
0.0Q0f15 3.0e-5Q0f15
-3.0e-5Q0f15 3.0e-5Q0f15
-3.0e-5Q0f15 3.0e-5Q0f15
-3.0e-5Q0f15 6.0e-5Q0f15
⋮
6.0e-5Q0f15 -6.0e-5Q0f15
0.0Q0f15 3.0e-5Q0f15
-6.0e-5Q0f15 3.0e-5Q0f15
-6.0e-5Q0f15 6.0e-5Q0f15
-6.0e-5Q0f15 6.0e-5Q0f15
-6.0e-5Q0f15 3.0e-5Q0f15
julia> x.samperate # サンプリング周波数
44100.0
波形の表示
音の加工などをする前に波形を見てみましょう。
x1 = x.data[:, 1];
t = (0 : length(x1)-1) / x.samplerate;
using PyPlot
plot(t, x1);
ylim(-1, +1);
grid();
title("Music Excerpt")
xlabel("Time (s)")
ylabel("Amplitude")
ひとまずはここまで。Octaveを使っているのとそれほど変わりない感覚です。
応用例(周波数分析)
おまけとして、リコーダーの音の倍音構造を見てみようと思います。変数名がぐちゃぐちゃだけど許して。
y = load("recorder.flac");
fs = y.samplerate;
y1 = y[1:Int(fs*1.8)];
y2 = zeros(nextpow(2, length(y1)));
y2[1:length(y1)] = y1;
Y2 = abs.(fft(y2));
Y2 = (Y2 / length(y2)) .^ 2;
Y3 = Y2[1 : Int(length(Y2)/2) + 1];
Y3[2:end-1] = Y3[2:end-1] * 2.0;
Y4 = 10*log10.(Y3);
f = linspace(0, y.samplerate/2, length(Y3));
plot(f, Y4-maximum(Y4));
xlim(0, 15000);
ylim(-80, +3);
grid();
title("Power Spectrum");
xlabel("Frequency (Hz)");
ylabel("Power (dB)");
savefig("output3.png");
この音は第2・第4倍音が少なめですね。