この記事はJulia Advent Calendar 2018 の24日目です。
僕は日頃、おもにオーディオ信号を扱うのにJuliaを使用しています。今日はそのうち楽器音の分析のために作成したプログラムを紹介します。
楽器音の分析
楽器の音を聞いたとき、トランペットとピアノとでは同じ大きさ・同じ高さで鳴らされた音であっても異なる楽器だと区別できます。これは音色が異なるからです。音色は「音の大きさと高さをそろえても、なお異なる聞こえ方をするとき、その違いに相当するもの」というような定義がされています。その音色に対応する要因としては、含まれる成分音のレベルや時間変化などが考えられます。
楽器音に含まれている振動を(フーリエ変換 などを利用して)周波数ごとに取り出すと「それぞれの周波数において、どのくらいの強さの音が含まれているのか」が分かります。これを周波数スペクトル (縦軸をどのようにとるかで振幅スペクトル、パワースペクトル などとも)と言います。成分音の中で低い周波数に突出しているのが基音で、ピッチ感のある音については基音周波数の整数倍の周波数に倍音 (あるいは高調波成分)があることがほとんどです。多くの場合は基音の周波数に対応したピッチを感じます(基音が最も強いエネルギーを持っている必要はありません)。打楽器のようにピッチ感が曖昧な音は、整数倍の音列になっていないことが原因のようです。
同じような音が続く(定常的な音の)場合は周波数分析をして倍音 がどのように含まれているかを見れば良いのですが、楽器の音の場合は弾き始めの音の立ち上がりや減衰部分など、様々な時間変化があります。そのとき、時間に応じて倍音 の構成なども変わる場合があります。そのため、時間波形の概形を見るためにエンベロープ (包絡線)を計算したりします。
時間波形全体の変化ではなく周波数成分の時間変化を見たいという時もあります。そこで登場するのがスペクトログラムです。主に短時間フーリエ変換 を使用し、その名の通り短い時間範囲をフーリエ変換 してその範囲のスペクトルを求めることを信号全体にわたって行うものです。下の図にその様子を描きました。横軸が時間、縦軸がパワー、手前から奥に向かって周波数としています。各時刻において、一定の時間幅についての周波数分析を行っており、その1フレームだけを取り出すと右側に取り出したようにスペクトルになっています。
フーリエ変換 では分析にある程度の時間区間 を必要とし、その時間の長さが分析できる周波数分解能に関わってきます。たとえば、4 Hzごとの波の振幅が知りたいときには、1 ÷ 4 = 0.25 秒の時間が必要です。一方で周波数成分の0.01秒単位の時間変化を見ようと思った場合には、 1 ÷ 0.01 = 100 Hzの周波数幅でしか分析ができません。
そういった弱点はありますが、以下ではスペクトログラムを見ていきましょう。
Juliaでのスペクトログラムの計算
これ以降のコードは次の環境で実行したものです。
Julia Version 1.0.3
Commit 099e826241 (2018-12-18 01:34 UTC)
Platform Info:
OS: macOS (x86_64-apple-darwin14.5.0)
CPU: Intel(R) Core(TM) i7-7920HQ CPU @ 3.10GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-6.0.0 (ORCJIT, skylake)
まずは必要になるパッケージ群をusing
します。
using FFTW, DSP, SampledSignals
using FileIO: load
import LibSndFile
using Plots
pyplot();
using Printf
LibSndFileを使ってWAVに保存された音ファイルを読み込みます。file
には拡張子まで入っていますが、拡張子を除いた部分をfilename
に入れて、これ以降で使用できるようにします。セミ コロンなしでload
すると、音の概形が表示されます。コマンドライン では実効値が、IJuliaでは下のような再生可能なインターフェイス が表示されます(下は画像ですので再生できません)。
file = "ching.wav" ;
filename = splitext(file)[1 ];
snd = load(" $( filename) .wav" )
音オブジエクトには標本化周波数(samplerate
)と音のサンプルデータ(data
)が入っています。それぞれを取り出して、時間波形をプロットしてみます。0.2秒付近で立ち上がり、そのあと急速に減衰していく様子が見えます。一定間隔で振幅が波打って増減している部分がありますが、うなりが発生しているのでしょうか。
fs = snd.samplerate;
x = vec(snd.data[: ,1 ]);
t = range(0 , stop= length(x)- 1 ) / fs;
plot(t, x, ylim= (- 1 , + 1 ),
xlabel= "Time (sec)" , ylabel= "Amplitude" , legend= false )
スペクトログラムを作成してみましょう。以下のようにDSP .jlのspectrogram
を使用すると簡単です。引数は
分析対象の信号
分析窓長
オーパ ーラップ量(ここでは窓の移動量(hopsize
)を定義しておいてwinsize-hopsize
と指定しています)
FFT 長...窓あたりのフーリエ変換 をどのくらいの長さで行うか(長く指定するとゼロパッドが行われます)
標本化周波数
窓の種類(今回はハニング窓)
の6つです。
winsize = 4096 ;
hopsize = 256 ;
S = spectrogram(x, winsize, winsize- hopsize,
nfft= winsize* 2 , fs= fs, window= hanning);
spectrogram
からの出力(S
)にはpower
、freq
、time
が入っています。power
は、【freq
の長さ】×【time
の長さ】の2次元行列で、各時刻・各周波数でのパワーが格納されています。$L=10\log_{10}\left(I\right)$を使ってデシベル に変換し、ヒートマップ表示してみます。
heatmap(S.time, S.freq, 10 * log10.(S.power),
xlabel= "Time (sec)" , ylabel= "Frequency (Hz)" )
これが音についての一般的なスペクトログラムで、横軸に時間、縦軸に周波数、色でパワーを表します。明るい部分に強いパワーがあるので、たとえば4000 Hz付近に横長に伸びる明るい線が基音を示しています。また、突出した成分音が整数次でない場所に出ていることも分かります(体鳴楽器 に多く見られる特徴です)。
アニメーションを作る
上記のスペクトログラムを、一枚ずつ次々と表示していけばアニメーションさせることもできます。このほうが時間変化が直感的に把握しやすいという人もいるかもしれません。
以下で表示したい周波数範囲のみを取り出す処理をしています。楽器によって音域が違うので、分析対象の楽器にあわせてこの部分を調整します。今回は0 Hzから分析できる上限であるナイキスト周波数 (標本化周波数の半分)までを分析するように指定しています。
Spow = S.power ./ maximum(S.power);
fL = 0 ;
fH = fs/ 2 ;
Sf = S.freq[fL .<= S.freq .<= fH];
Spow = Spow[fL .<= S.freq .<= fH, : ];
St = S.time;
アニメーションを作るのには@animate
を使います。@animate
つきのfor
文の中で一枚一枚をプロットして、anim
に格納された画像をアニメGIFファイルとして保存します。gif()
の引数にfps=15
とあるのは1秒あたりのフレーム数(frames per second)で、あまり大きくない整数値であればOKのようです。
(実行すると結構時間がかかります! 自分の環境(MacBook Pro 2017, 3.1 GHz Intel Core i7 )では2分ほどかかりました)
@time anim = @animate for fr= 1 : length(St)
timetext = @sprintf (" %.2f " , St[fr]);
plot(Sf, 10 * log10.(Spow[: ,fr]),
xlim= (fL, fH), ylim= (- 105 , + 5 ),
title= " $( filename) , t= $( timetext) sec" ,
xlabel= "Frequency (Hz)" , ylabel= "Power (dB)" ,
legend= false , size= (720 , 480 ));
end
gif(anim, "tmp_ $( filename) .gif" , fps= 15 );
run(`ffmpeg -loglevel quiet -y -r 24 -i tmp_ $( filename) .gif -r 24 -vcodec libx264 -pix_fmt yuv420p $( filename) .mp4` );
なお、run
命令を使って外部プログラムのffmpeg を呼び出し、アニメーションGIFをMP4に変換しています。GIFだけでいい場合はrun
の行を省略してください。完成したものを以下にアップロードしました。
VIDEO
蛇足
今回はDSP .jlのspectrogram
を使用しましたが、スペクトログラムのプログラム(短時間フーリエ変換 )を自作しようとすると以下のような感じになるでしょうか。ただし、簡易版なのでspectrogram
のように機能豊富ではなく、必要最低限のものしか入っていません。参考までに。
function spect (x, fs, window_size= 1024 , hop_size= 256 )
num_frames = ceil(Int , (length(x) + window_size- 1 ) / hop_size);
xx = [x ; zeros(2 * window_size)];
N = nextpow(2 , window_size);
zeropad = zeros(N - window_size);
y = zeros(Complex , round(Int , N/ 2 + 1 ), num_frames);
for n= 1 : num_frames
xxx = xx[hop_size* (n- 1 )+ 1 : hop_size* (n- 1 )+ window_size];
y[: ,n] = rfft([xxx ; zeropad]) / N;
end
t = ((0 : num_frames- 1 ) * hop_size) / fs;
f = range(0 , stop= fs/ 2 , length= round(Int , N/ 2 + 1 ));
return (y, t, f);
end