Go to Qiita Advent Calendar 2024 Top
LoginSignup
29
15

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

AizuAdvent Calendar 2016

Day 4
@stringamp

Rubyを使ったWAVファイルのBPM解析

Last updated at Posted at 2016-12-04

この記事は Aizu Advent Calendar 2016 4日目の記事です。

前の人は、@hnjkさん、次の人は @misoton665 さんです。

はじめに

この記事では汎用的な音楽ファイル形式であるWAVファイルについて、BPM解析を行うプログラムを作成します。
今回は手軽に書けるのと、最近あまり注力してRubyを書いてないので、リハビリという名目でRubyを使って書きました。特に技術的な理由はありません。

$ ruby -v
-> ruby 2.2.3p173

作成したプログラムへの指摘は大歓迎です。

WAVファイル形式について

BPMを解析するプログラムだとかアルゴリズムを調べる前に、まず今回解析の対象とするWAVファイルの形式について、基礎知識から調べてみました。

WAVまたはWAVE（ウェーブ、ウェブ） (RIFF waveform Audio Format) は、マイクロソフトとIBMにより開発された音声データ記述のためのフォーマットである。RIFFの一種。主としてWindowsで使われるファイル形式である。ファイルに格納した場合の拡張子は、.wav。
WAV - Wikipedia, https://ja.wikipedia.org/wiki/WAV

マイクロソフトとIBM開発だったんですね…。ここでRIFFとは？と思いさらに調べると、

Resource Interchange File Format（RIFF、「資源交換用ファイル形式」の意味）は、タグ付きのデータを格納するための汎用メタファイル形式である。
Resource Interchange File Format - Wikipedia, https://ja.wikipedia.org/wiki/Resource_Interchange_File_Format

とのこと。上のWikiの内容を簡単にまとめると、RIFFに基づいているファイルは小さなチャンクから構成されており、それぞれINFOチャンク、DATAチャンクなどから構成されている模様。つまりRIFFの一種であるWAVファイルも同じようなチャンクから成る構成をしているということでしょう。
では具体的にそのチャンクの内容とはどうなっているのでしょう。
WAVフォーマットのファイル内チャンクの構成について調べたところ、以下のようなサイトを見つけました。

このサイトの説明によると、WAVファイルはRIFFに基づき最初の数バイトはINFO情報（後に続くDATAチャンクの基本の情報）で、更にその後に続くnバイトがDATAチャンク（ファイル本体の情報）だそうです。
ということはそれぞれのチャンクの情報が取り出せれば取り敢えず解析は出来そうですね。

BPM解析手法

http://hp.vector.co.jp/authors/VA046927/tempo/tempo.html
理論的にはこのサイトでかなり詳しく説明されているようなのでこの場で込み入った説明はしません。サイトから引用すると、WAVファイルのBPMを求めるまでの手順は、

  • STEP 1. WAVファイルを一定時間（フレーム）ごとに区切る。
  • STEP 2. フレームごとの音量を求める。
  • STEP 3. 隣り合うフレームの音量の増加量を求める。
  • STEP 4. 増加量の時間変化の周波数成分を求める。
  • STEP 5. 周波数成分のピークを検出する。
  • STEP 6. ピークの周波数からテンポを計算する。

と幾つかのSTEPに分けられるそうです。なるほど？
BPMを求める際、ここではBPMは60から240の間にあることを仮定して実行しているようなので、この記事でも同じ仮定の上でプログラムを書いていきます。
さて、上の処理は大きく分けて前半と後半に分けることが出来そうです。具体的には、WAVファイルに対してデータ本体に処理を行う

STEP 1. WAVファイルを一定時間(以下フレーム)ごとに区切る。
STEP 2. フレームごとの音量を求める。
STEP 3. 隣り合うフレームの音量の増加量を求める。

このパートと、実際に得られた結果から

STEP 4. 増加量の時間変化の周波数成分を求める。
STEP 5. 周波数成分のピークを検出する。
STEP 6. ピークの周波数からテンポを計算する。

以上の項目を導出して行くパートです。
実際に次の実装の節では前半後半と二つに分けて実装して行きます。

実装

前半

では実際に上で挙げたBPM解析の手法に基づいてプログラムを書いて行きましょう。一からWAVを解析するようなプログラムを書いている程時間がなかったので、Wavファイルを扱えるwav-fileと言うgemを使用します。早速導入。

gem install wav-file

このライブラリでは、WAV形式のINFOチャンク部分とDATAチャンク部分を分けて出力出来るようだったので、適当な曲を拾って来て実際にINFOチャンクの情報を出力してみました。

require 'wav-file'

f = open("./Take-Her-Heart.wav")
format = WavFile::readFormat(f)
f.close

puts format

これだけでWAVファイルのINFOチャンクから情報を読み取り、出力してくれるそうです。実際に実行してみると、

スクリーンショット 2016-12-03 10.22.50.png

ちゃんと出力されました。すごい簡単。
動作が確認できた所で、実際にSTEPごとに実装して行ってみます。まず最初は、

STEP 1. WAVファイルを一定時間(以下フレーム)ごとに区切る。

との事でした。wav-fileを開発している方のブログ記事によると、WAVファイルのDATAチャンクバイナリから、WAVの波形データを配列 wavs として取得するには、

f = open("input.wav")
format = WavFile::readFormat(f)
dataChunk = WavFile::readDataChunk(f)
f.close

dataChunk = WavFile::readDataChunk(f)
bit = 's*' if format.bitPerSample == 16 # int16_t
bit = 'c*' if format.bitPerSample == 8 # signed char
wavs = dataChunk.data.unpack(bit) # read binary

と書くそうなので、DATAチャンクから波形データを取り出し、波形データを一定フレーム数（ここでは512にしています）に区切るプログラムは以下のように書けるはずです。

FRAME_LEN = 512

def bit_per_sample(format)
  # ここでは16bitか8bitしか対象にして居ないので三項演算子で表現
  format.bitPerSample == 16 ? 's*' : 'c*'
end

def get_wav_array(data_chunk, format)
  data_chunk.data.unpack(bit_per_sample(format))
end

f = open("./tempo_120.wav") # BPM120のメトロノーム音が入ったWAVファイル
data_chunk = WavFile::readDataChunk(f)
format = WavFile::readFormat(f)
f.close()

get_wav_array(data_chunk, format)
      .take(@wavs.size - @wavs.size % FRAME_LEN)
      .each_slice(FRAME_LEN).to_a

ここでは、Array#takeメソッドを使い、波形データ末尾にあるフレーム内データ数512に満たない余ったデータを切り捨て、
Enumerable#each_sliceとto_aメソッドでフレーム毎の配列の配列へと変換しています。

[[フレーム1の波形データ配列], [フレーム2の波形データ配列], ..., [フレームnの波形データ配列]]

# Enumerable#each_sliceメソッドは便利なメソッドなのですが、Arrayの高階関数として呼び出す事が出来ても戻り値の型がEnumerable型なので、メソッドチェーンを繋げていく場合はいちいちto_aしないとダメで結構面倒…。良い方法がないかな。

次に行きましょう。

STEP 2. フレームごとの音量を求める。

直前で挙げたように、現在WAVの波形データをフレーム毎に区切った配列

[[フレーム1の波形データ配列], [フレーム2の波形データ配列], ..., [フレームnの波形データ配列]]

が手に入っている状況ですので、ここでそれぞれのフレームに対して、そのフレーム内の波形データの二乗平均平方根（Root Mean Square, RMS）を計算する事でフレームごとの音量を求めたいと思います。

# 二乗平均平方根についての説明はコチラ。

これは、先ほどの get_wav_array(data_chunk, format) に続けて、

get_wav_array(data_chunk, format).take(@wavs.size - @wavs.size % FRAME_LEN)
      .each_slice(FRAME_LEN).to_a
      .map{|arr| Math.sqrt(arr.map{|elem| elem ** 2}.inject(:+) / arr.size)}

と書く事が出来るでしょう。最後に追加したArray#mapの中身の

Math.sqrt(arr.map{|elem| elem ** 2}.inject(:+) / arr.size)

この部分が、

この数式の右辺と対応しています。
これで

[フレーム1のRMS値, フレーム2のRMS値, フレーム3のRMS値, ..., フレームnのRMS値]

が手に入りました。
どんどん行きます。次にやるべきことは、

STEP 3. 隣り合うフレームの音量の増加量を求める。

とのことでした。つまり、先ほどの配列を、

[フレーム1と2のRMS値差, フレーム2と3のRMS値差, ..., フレームn-1とnのRMS値差]

という形に変換すれば良い訳です。ここで、変換前と変換後で配列の要素数が1つ減少している点に注意です。1つ前のstepで手に入った配列を、 diff_arr としておきましょう。すると、

diff_arr[0..-2].zip(diff_arr[1..-1]).map{|f,x| f - x < 0 ? 0 : f - x}

と書けます。上のプログラムは単純に、今得られている配列 diff_arr の1番目の要素からn-1番目の要素から、同じく diff_arr の2番目の要素からn番目の要素を引き算しているだけです。ここで、知りたいのは増加量であり、減少量については無視するので、引き算の結果がマイナスの場合は増加量を 0 としています。

# Rubyにはvector演算を行うメソッドが存在しないので、Array#zipで二つの配列をまとめ、mapでまとめられた値に対する演算を書く方式が主流のようです。
# しかしここまでメソッドチェーンを繋げた以上、キリの良いところまで繋げたかった…！

長々と書きましたが、一旦ここまでをまとめると以下になります。

diff_arr = @wavs.take(@wavs.size - @wavs.size % FRAME_LEN)
      .each_slice(FRAME_LEN).to_a
      .map{|arr| Math.sqrt(arr.map{|elem| elem ** 2}.inject(:+) / arr.size)}

diff_arr[0..-2].zip(diff_arr[1..-1]).map{|f,x| f - x < 0 ? 0 : f - x}

後半

ここから少し理論的に難しくなっていきます。まず、

STEP 4. 増加量の時間変化の周波数成分を求める

これを考えます。サイトによると、n番目のフレームの音量の増加量をD(n)、フレームのサンプリング周波数をsとすると、各BPMのマッチ度Rbpmは以下で表されるとのこと。
スクリーンショット 2016-12-03 12.08.27.png

この数式をRubyに落とし込みましょう。
まず最初に今回対象としているBPMの幅は60~240なので、BPMをキーとして、値がそのBPMのマッチ度となるようなハッシュを作成しておきます。

(60..240).inject({}){|acc, bpm|
      acc[bpm] = calc_bpm_match(data, bpm)
      acc
}

ハッシュを作成する過程で呼び出されているcalc_bpm_match関数は、渡されたbpmとdata（前半で取得したフレーム毎のRMS差配列）のマッチ度を求める関数です。
以下で順番に書いて行ってみます。

さて、最初の式
スクリーンショット 2016-12-03 12.21.19.png
こちらに注目。
これは簡単で、Ruby上では以下のように表すことができるでしょう。

Math.sqrt(a_bpm ** 2 + b_bpm ** 2)

では、次に a_bpm と b_bpm について考えてみます。
数式は、

スクリーンショット 2016-12-03 12.21.25.png スクリーンショット 2016-12-03 12.21.31.png

となっていますが、二つの式に大きな差異は無く、cosとsinのみが違っている定義になっています。
ここで、cosとsinの偏角部分の数式が少々複雑になっているので、 a_bpm と b_bpm を定義する前に、lambda式を定義しておくことにします。

phase_cos = lambda{|m| Math.cos(2 * Math::PI * f_bpm * m / SAMPLE_F_PER_FRAME)}
phase_sin = lambda{|m| Math.sin(2 * Math::PI * f_bpm * m / SAMPLE_F_PER_FRAME)}

これらはそれぞれ必要な時に .apply(m) もしくは .(m) によって呼び出して利用します。
さて、それでは a_bpm と b_bpm の定義に移りましょう。Ruby上の数式に直すと以下のようになると考えられます。

a_bpm = (0..data.size-1)
      .map{|m| phase_cos.(m)}
      .zip(data)
      .map{|x,y| x * y}
      .inject(:+) / data.size

a_bpm = (0..data.size-1)
      .map{|m| phase_sin.(m)}
      .zip(data)
      .map{|x,y| x * y}
      .inject(:+) / data.size

ここで、先ほども出たようにベクトル同士の乗算が出てきますが、ここでもベクトル同士の演算を行う際にはzipしてmapして演算しています。しかしこれを書くのが二回目となると流石にいい方法が無いか気になってきましたので、twitterで呟いてみると…

なんとzipWithをRubyで実装している方がいらっしゃる事を教えてもらいました。

この記事の内容によると、Enumerableモジュールに自前のzip_withメソッドを追加する事で実現しているようです。

module Enumerable
  def zip_with(*others, &block)
    zip(*others).map &block
  end
end

ここで定義したzip_withを利用すると、上で定義した a_bpm と b_bpm は、以下のように書き換えられます。

a_bpm = (0..data.size-1)
      .map{|m| phase_cos.(m)}
      .zip_with(data){|x,y| x * y}
      .inject(:+) / data.size

b_bpm = (0..data.size-1)
      .map{|m| phase_sin.(m)}
      .zip_with(data){|x,y| x * y}
      .inject(:+) / data.size

大してコード量が少なくなったわけでもないですが、個人的に可読性も上がって満足。ついでにSTEP3のzip+mapも書き換えちゃいましょう。

diff_arr = @wavs.take(@wavs.size - @wavs.size % FRAME_LEN)
      .each_slice(FRAME_LEN).to_a
      .map{|arr| Math.sqrt(arr.map{|elem| elem ** 2}.inject(:+) / arr.size)}

diff_arr[0..-2].zip_with(diff_arr[1..-1]){|f,x| f - x < 0 ? 0 : f - x}

以上を全て組み合わせ、最終的なcalc_bpm_match関数は以下のようになります。

def calc_bpm_match(data, bpm)
      f_bpm = bpm / 60.0

      phase_cos = lambda{|m| Math.cos(2 * Math::PI * f_bpm * m / SAMPLE_F_PER_FRAME)}
      phase_sin = lambda{|m| Math.sin(2 * Math::PI * f_bpm * m / SAMPLE_F_PER_FRAME)}

      a_bpm = (0..data.size-1)
        .map{|m| phase_cos.(m)}
        .zip_with(data){|x,y| x * y}
        .inject(:+) / data.size
      
      b_bpm = (0..data.size-1)
        .map{|m| phase_sin.(m)}
        .zip_with(data){|x,y| x * y}
        .inject(:+) / data.size

      Math.sqrt(a_bpm ** 2 + b_bpm ** 2)
end

この辺でコード量が増え、結構まとまりが無い雰囲気が漂って来たので、次のステップに移る前にclassにまとめてしまいます。

class BPMAnalyzer
  FRAME_LEN = 512
  SAMPLE_F_PER_FRAME = 44100.0 / FRAME_LEN

  def initialize(file_name)
    f = open(file_name)
    @format = WavFile::readFormat(f)
    @data_chunk = WavFile::readDataChunk(f)
    @wavs = get_wav_array(@data_chunk, @format)
    f.close
    @res = nil
  end

  def run
    diff_arr = @wavs.take(@wavs.size - @wavs.size % FRAME_LEN)
      .each_slice(FRAME_LEN).to_a
      .map{|arr| Math.sqrt(arr.map{|elem| elem ** 2}.inject(:+) / arr.size)}

    diff_arr[0..-2].zip(diff_arr[1..-1]).map{|f,x| f - x < 0 ? 0 : f - x}

    @res = calc_match(diff_arr)
  end

  def to_s
    "BPM,Match rate\n" + @res.map{|k, v| "#{k},#{v}"}.join("\n") if @res != nil
  end

private
  def bit_per_sample(format)
    format.bitPerSample == 16 ? 's*' : 'c*'
  end

  def get_wav_array(data_chunk, format)
    data_chunk.data.unpack(bit_per_sample(format))
  end

  def calc_bpm_match(data, bpm)
      f_bpm = bpm / 60.0

      phase_cos = lambda{|m| Math.cos(2 * Math::PI * f_bpm * m / SAMPLE_F_PER_FRAME)}
      phase_sin = lambda{|m| Math.sin(2 * Math::PI * f_bpm * m / SAMPLE_F_PER_FRAME)}

      a_bpm = (0..data.size-1)
        .map{|m| phase_cos.(m)}
        .zip_with(data){|x,y| x * y}
        .inject(:+) / data.size

      b_bpm = (0..data.size-1)
        .map{|m| phase_sin.(m)}
        .zip_with(data){|x,y| x * y}
        .inject(:+) / data.size

      Math.sqrt(a_bpm ** 2 + b_bpm ** 2)
  end

  def calc_match(data)
      (60..240).inject({}){|acc, bpm|
        acc[bpm] = calc_bpm_match(data, bpm)
        acc
      }
  end
end

to_s メソッドは単純にBPM解析した結果を標準出力するだけです。Rなどの統計ソフトで使いやすいように csv の形式にしてあります。
さて、オブジェクトの形に直したところで次のSTEPに移りましょう。ここまでくればあと一息っぽいです。

STEP 5. 周波数成分のピークを検出する。
STEP 6. ピークの周波数からテンポを計算する。

この項目は内容がほぼほぼ被っているので一気にやってしまいましょう。先ほど作成した60から240までのBPMと、そのマッチ度を利用してピークを検出します。コードは以下のように書きました。

# add public method to BPMAnalyzer class
def get_max_rate
  "BPM,Match rate\n" + @res.max{|k, v| k[1] <=> v[1]}.join(",") if @res != nil
end

analyzer = BPMAnalyzer.new(ARGV[0])
analyzer.run
puts analyzer.get_max_rate

この時点での出力は、
スクリーンショット 2016-12-03 15.35.49.png
となりました。どうやら正しくBPMが検出されているようです。

試しに幾つかの曲を用意して解析してみます。

スクリーンショット 2016-12-03 16.12.03.png

成功。音圧の違いで検出がしやすいのかと考え、静かめの曲で再実行

スクリーンショット 2016-12-03 15.55.31.png

大人しめの、激しい音圧の上下が無い曲でも無事正しく検出できました。

まとめ

BPM検出プログラムをRubyで、なるべく手続き的にならないように書いてみました。良い復習になったと思います。
プログラムにもまだまだ改善の余地がありそうなのと、時間があれば今回の結果を利用してBPMが途中で変わる曲とかも判別できたら面白いかな、とか考えてます。音響処理関係のプログラミングに触れる良い機会になりました。
最後まで見て頂きありがとうございました。

GitHubリポジトリ

参考

http://sky.geocities.jp/kmaedam/directx9/waveform.html
http://hp.vector.co.jp/authors/VA046927/tempo/tempo.html
http://ism1000ch.hatenablog.com/entry/2014/07/08/164124
http://ackiesound.ifdef.jp/doc/tempo/main.html

AIZU ADVENT CALENDAR 2016 明日12/5は @misoton665 さんです。よろしくお願いします。

29
15
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
29
15

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?