楽器音の分析とスペクトログラム (Julia Advent Calendar 2018)

この記事は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

# WAVファイルの読み込みに使用
using FileIO: load
import LibSndFile

# グラフの作成はPyPlotを使う
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を使用すると簡単です。引数は

  1. 分析対象の信号
  2. 分析窓長
  3. オーパーラップ量(ここでは窓の移動量(hopsize)を定義しておいてwinsize-hopsizeと指定しています)
  4. FFT長...窓あたりのフーリエ変換をどのくらいの長さで行うか(長く指定するとゼロパッドが行われます)
  5. 標本化周波数
  6. 窓の種類(今回はハニング窓)

の6つです。

# compute spectrogram
winsize = 4096;
hopsize = 256;
S = spectrogram(x, winsize, winsize-hopsize,
    nfft=winsize*2, fs=fs, window=hanning);

spectrogramからの出力(S)にはpowerfreqtimeが入っています。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から分析できる上限であるナイキスト周波数(標本化周波数の半分)までを分析するように指定しています。

# normalize the power (maximum will become 0 dBFS)
Spow = S.power ./ maximum(S.power);

# extract only the frequency range we need
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分ほどかかりました)

# create animation
@time anim = @animate for fr=1:length(St)
  timetext = @sprintf("%.2f", St[fr]);

  # plot spectrum
  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の行を省略してください。完成したものを以下にアップロードしました。


蛇足

今回は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