音ファイルの読み込みと単純な可視化

ここ数年の僕はプログラミングというと音と統計ばかりやっているのですが、特に音の分析や加工に関してよく使うものは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];   # チャンネル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");   # FLACも読み込める
fs = y.samplerate;
y1 = y[1:Int(fs*1.8)];   # 冒頭から1.8秒分だけ使う（約524Hzの吹奏音が入っています）

y2 = zeros(nextpow(2, length(y1)));   # y1より長い、2の累乗の長さの零ベクトルを作る
y2[1:length(y1)] = y1;   # その先頭にy1を入れる（つまりzero-padding）
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));   # 最大値が0 dBになるようにしてグラフ描画
xlim(0, 15000);   # X軸・Y軸の範囲指定
ylim(-80, +3);
grid();   # グリッドを描画
title("Power Spectrum");   # タイトルや軸のラベルを付ける
xlabel("Frequency (Hz)");
ylabel("Power (dB)");

savefig("output3.png");   # PNGファイルへの保存

この音は第2・第4倍音が少なめですね。