類似楽曲検索システムを作ろう
もう1年以上かけて音声信号処理の勉強をしてきました（Pythonで音声信号処理）。ここらで具体的なアプリケーションとして類似楽曲検索の実験をしてみたのでレポートをまとめておきます。言語はPythonです。
前に
- 類似画像検索システムを作ろう（2009/10/3）
- Visual Wordsを用いた類似画像検索（2010/2/27）
という画像の類似検索に関するエントリを書きましたが、今回は画像ではなく音楽を対象に類似検索をやってみたいと思います！
今回作る類似楽曲検索システムは、従来からよくあるアーティスト名や曲名などテキストで検索するシステムや購買履歴をもとにオススメする協調フィルタリングベースのシステムとは異なります。WAVEファイルやMP3ファイルなどの音楽波形そのものを入力とするのが特徴です。たとえば、「具体的なアーティストや曲名は知らないけれど、この曲とメロディや雰囲気が似た曲がほしいな」なんていうときに便利に使えるシステムです。
この研究分野は、Music Information Retrieval (MIR) と呼ばれています。私自身、アーティストによって曲を買うことはほとんどないので、以前からこういうシステムがほしいなぁと思ってました。
とりあえず実験結果だけ見たい方は下のページにあります。1位がクエリの楽曲で似ている順に上位10件を出力しました。HTML5でオーディオ断片を埋め込んだのでFirefoxやChromeだと音声で確認できます。
どうでしょう？似てますかね？
概要
今回の実験で参考にした資料は、
- E. Pampalk: Computational Models of Music Similarity and their Application in Music Information Retrieval (PDF), Doctoral Thesis, Vienna University of Technology, Austria, March 2006.
- B. Logan and A. Salomon: A Content-based Music Similarity Function, Technical report, Compaq Cambridge Research Lab, 2001.
という2つの論文です。上の論文は、類似楽曲検索のサーベイ論文としてよくまとまっています。今回はPampalkさんの論文で紹介されているLoganさんの論文の手法を実際に試してみました。
Loganさんが提案したのは音声の局所特徴量を用いたアプローチでVisual Wordsを用いた類似画像検索（2010/2/27）と考え方がよく似ています。類似画像検索だとSIFTやSURFといった局所特徴量をクラスタリングしてVisual wordsを作りました。そして、各画像をVisual wordsのヒストグラムで表現し、ヒストグラム間の距離をヒストグラムインターセクションで求めました。
今回の音楽の場合はメル周波数ケプストラム係数（2012/2/25）、略してMFCCというスペクトル特徴量を使います。各楽曲から得られるMFCCベクトル集合をクラスタリングして楽曲をクラスタの集合で表現し、クラスタ集合間の距離をEarth Mover's Distance（2012/8/4）で求めます。EMDが近い楽曲が互いに類似した楽曲となるわけです。
コンテンツの表現方法や使用する距離の定義が多少違いますが、アプローチはよく似ています。というわけで論文にそって早速試してみます！
MP3ファイルの準備
類似楽曲検索なのでまず楽曲ファイルが必要です。最初、ドラゴンクエストやファイナルファンタジーのサウンドトラックで試してみたのですが、どうもよい結果が得られませんでした。いろいろ試したところインストだけの曲よりボーカルが入ってた方がよい結果になるみたいです。
私自身はボーカル入りのCDはあんまり持ってなかったので、妹のiTunesリポジトリをパクって借りて、500曲くらいのボーカル曲のリストを作りました。そんなわけで、検索結果の曲のすべてが私の趣味というわけではないのであしからず（笑）500曲のMP3ファイルは当然ながら公開できないのでお手持ちの曲でぜひ試してみてください。
iTunesはアルバム名やアーティスト名で階層構造になっていますが、ここでは簡単のためすべての曲をmp3というフォルダにまとめました。iTunesで買った曲などMP3以外のフォーマットがある場合は、MP3への変換が必要です。やり方はググってください。
MFCCの抽出
まず、MP3ファイルからスペクトル特徴量のメル周波数ケプストラム係数（MFCC）を抽出します。MFCCはスペクトルの概形を表すパラメータなので音色を表すと考えてよいと思います。MFCCに関しては、前に書いた記事を参考にしてください。
- メル周波数ケプストラム係数（MFCC）（2012/2/25）
- SPTKの使い方 (6) MFCCの抽出（2012/8/5）
以下が、指定したMP3ディレクトリ（mp3）にあるすべてのMP3ファイルからMFCCを抽出してMFCCディレクトリ（mfcc）に保存するPythonスクリプトです。音声フォーマットの変換にsoxとlame、波形の切り出し、MFCCの抽出にSPTKというツールを使っています。両方ともインストールしないと動きません。SPTKのインストールは
- SPTKの使い方 (1) インストール・波形描画・音声再生（2012/7/1）
という記事にまとめました。
#coding:utf-8 import os import sys # mp3_to_mfcc.py # usage: python mp3_to_mfcc.py [mp3dir] [mfccdir] [rawdir] # ディレクトリ内のMP3ファイルからMFCCを抽出する def mp3ToRaw(mp3File, rawFile): # mp3を16kHz, 32bitでリサンプリング os.system("lame --resample 16 -b 32 -a '%s' temp.mp3" % mp3File) # mp3をwavに変換 os.system("lame --decode temp.mp3 temp.wav") # wavをrawに変換 os.system("sox temp.wav %s" % rawFile) os.remove("temp.mp3") os.remove("temp.wav") def calcNumSample(rawFile): # 1サンプルはshort型（2byte）なのでファイルサイズを2で割る filesize = os.path.getsize("temp.raw") numsample = filesize / 2 return numsample def extractCenter(inFile, outFile, period): # 波形のサンプル数を求める numsample = calcNumSample(inFile) fs = 16000 center = numsample / 2 start = center - fs * period end = center + fs * period # period*2秒未満の場合は範囲を狭める if start < 0: start = 0 if end > numsample - 1: end = numsample - 1 # SPTKのbcutコマンドで切り出す os.system("bcut +s -s %d -e %d < '%s' > '%s'" \ % (start, end, "temp.raw", rawFile)) def calcMFCC(rawFile, mfccFile): # サンプリング周波数: 16kHz # フレーム長: 400サンプル # シフト幅 : 160サンプル # チャンネル数: 40 # MFCC: 19次元 + エネルギー os.system("x2x +sf < '%s' | frame -l 400 -p 160 | mfcc -l 400 -f 16 -n 40 -m 19 -E > '%s'" % (rawFile, mfccFile)) if __name__ == "__main__": if len(sys.argv) != 4: print "usage: python mp3_to_mfcc.py [mp3dir] [mfccdir] [rawdir]" sys.exit() mp3Dir = sys.argv[1] mfccDir = sys.argv[2] rawDir = sys.argv[3] if not os.path.exists(mfccDir): os.mkdir(mfccDir) if not os.path.exists(rawDir): os.mkdir(rawDir) for file in os.listdir(mp3Dir): if not file.endswith(".mp3"): continue mp3File = os.path.join(mp3Dir, file) mfccFile = os.path.join(mfccDir, file.replace(".mp3", ".mfc")) rawFile = os.path.join(rawDir, file.replace(".mp3", ".raw")) try: # MP3を変換 mp3ToRaw(mp3File, "temp.raw") # 中央の30秒だけ抽出してrawFileへ extractCenter("temp.raw", rawFile, 15) # MFCCを計算 calcMFCC(rawFile, mfccFile) print "%s => %s" % (mp3File, mfccFile) # 後片付け os.remove("temp.raw") except: continue
曲全体からMFCCを抽出すると処理が重すぎるので、Pampalkさんの論文に従って、曲の中心の30秒間だけを対象にしました。本当は、画像の局所特徴量SIFTのようにその曲の特徴をよく表す部分（サビ？）からとりたかったのだけどどうやればいいんだろう。
MFCCファイルは、SPTKのバイナリ形式になっています。前に書いたprint_mfcc.pyというスクリプトでダンプできます。愛は勝つ.mfcを20次元ずつダンプしてみます。
> python print_mfcc.py mfcc/愛は勝つ.mfc 20 -14.77 -8.78 -6.70 ... 3.45 7.37 1.44 2.80 22.29 -13.84 -9.71 -3.28 ... -0.82 6.72 2.70 1.92 22.85 -15.45 -12.54 -7.64 ... -5.64 3.45 2.31 6.03 22.67 -14.64 -10.72 -4.96 ... -5.24 3.64 0.71 5.49 22.78
各行が1フレームの20次元MFCCベクトル（19次元+パワー）です。サンプリング周波数16kHz、シフト幅160サンプル、30秒間の楽曲なので、30 x 16000 / 160 = 3000フレームあります。つまり、「愛は勝つ」のmfcファイルは20次元のMFCCベクトルが3000行もあります。1曲を表すベクトルとしては大きすぎですね・・・
シグネチャの計算
1曲を表す特徴ベクトルが、20次元ベクトル3000個（20x3000 = 60000個の浮動小数点数）では大きすぎるのでもう少し情報を圧縮します。こんなときに使えるのが前に紹介したベクトル量子化（2012/8/13）というアルゴリズムです。ようは3000個のベクトルをクラスタリングしてよく似たベクトルをクラスタにまとめてしまいます。そして各クラスタを正規分布と仮定してその平均ベクトルと分散共分散行列を求め、それをあらたな特徴量とします。式で表すと楽曲Pの特徴量は、
となります。mはクラスタの数です。
クラスタリングにはよく使われるk-meansを用いました。クラスタ数は上の論文に合わせて16としました。つまり、3000個のベクトル集合をたった16個のクラスタに分類します。各クラスタの平均ベクトルは20次元ベクトル、分散共分散行列は20x20次元の行列になります。つまり、
(20次元平均ベクトル + 20x20次元の分散共分散行列) x 16クラスタ = 6720個の浮動小数点数
まで情報を圧縮できます。実際は、各クラスタの平均と分散のほかに重み（そのクラスタに分類されたベクトルの数）の情報も加えるので正確には6736個の浮動小数点数になります。60000個に比べたらずいぶん減ります。下のスクリプトは、各局のMFCCファイルをシグネチャファイルに変換するスクリプトです。
#coding:utf-8 import os import struct import sys import numpy as np import scipy.cluster # mfcc_to_signature.py # usage: python mfcc_to_signature.py [mfccdir] [sigdir] # 各曲のMFCCをシグネチャに変換する def loadMFCC(mfccFile, m): """MFCCをロードする、mはMFCCの次元数""" mfcc = [] fp = open(mfccFile, "rb") while True: b = fp.read(4) if b == "": break val = struct.unpack("f", b)[0] mfcc.append(val) fp.close() # 各行がフレームのMFCC # numFrame行、m列の行列形式に変換 mfcc = np.array(mfcc) numFrame = len(mfcc) / m mfcc = mfcc.reshape(numFrame, m) return mfcc def vq(mfcc, k): """mfccのベクトル集合をk個のクラスタにベクトル量子化""" codebook, destortion = scipy.cluster.vq.kmeans(mfcc, k) code, dist = scipy.cluster.vq.vq(mfcc, codebook) return code if __name__ == "__main__": if len(sys.argv) != 3: print "usage: python mfcc_to_signature.py [mfccdir] [sigdir]" sys.exit() mfccDir = sys.argv[1] sigDir = sys.argv[2] if not os.path.exists(sigDir): os.mkdir(sigDir) for file in os.listdir(mfccDir): if not file.endswith(".mfc"): continue mfccFile = os.path.join(mfccDir, file) sigFile = os.path.join(sigDir, file.replace(".mfc", ".sig")) print mfccFile, "=>", sigFile fout = open(sigFile, "w") # MFCCをロード # 各行がフレームのMFCCベクトル mfcc = loadMFCC(mfccFile, 20) # MFCCをベクトル量子化してコードを求める code = vq(mfcc, 16) # 各クラスタのデータ数、平均ベクトル、 # 共分散行列を求めてシグネチャとする for k in range(16): # クラスタkのフレームのみ抽出 frames = np.array([mfcc[i] for i in range(len(mfcc)) if code[i] == k]) # MFCCの各次元の平均をとって平均ベクトルを求める m = np.apply_along_axis(np.mean, 0, frames) # 0は縦方向 # MFCCの各次元間での分散・共分散行列を求める sigma = np.cov(frames.T) # 重み（各クラスタのデータ数） w = len(frames) # このクラスタの特徴量をフラット形式で出力 # 1行が重み1個、平均ベクトル20個、分散・共分散行列400個の計421個の数値列 features = np.hstack((w, m, sigma.flatten())) features = [str(x) for x in features] fout.write(" ".join(features) + "\n") fout.close()
出力である楽曲のシグネチャファイルは、1行に (重み、20次元平均ベクトル、20x20次元の分散・共分散行列) の421個の数値が並びます。それが、16クラスタ分あるので16行のファイルです。
シグネチャ間の距離と楽曲間の距離
次に2つの楽曲のシグネチャ間の距離を求めます。これには、前に書いたEarth Mover's Distance（2012/8/4）を使います。
左側が楽曲Pのシグネチャで右側が楽曲Qのシグネチャです。楽曲Pと楽曲Qの距離を求めるのが課題です。EMDの詳しい解説は、Earth Mover's Distance（2012/8/4）を参照してください。先のページのEMD具体例では、各シグネチャが単純なベクトルなので、各特徴量間の距離にユークリッド距離が使えました。ですが、今回は各特徴量がベクトルではなく、正規分布になっています。というわけで単純なユークリッド距離が使えません・・・
正規分布間の距離がはかれる指標としてカルバック・ライブラー情報量（Kullback Leibler Divergence）というのがあります。Kullback-Leibler divergence - Wikipediaによると正規分布N1とN2（それぞれ平均ベクトル、分散共分散行列がある）の間のカルバック・ライブラー情報量の定義は、
です。うう、目がくらむぅ。kは平均ベクトルの次元数です。一般的にカルバック・ライブラー情報量は対称性が成り立ちません。つまり、と
で異なる値が出てきます。そのため、上の式をそのまま使うと楽曲Pからみた楽曲Qの距離と楽曲Qからみた楽曲Pの距離が違ってきます。これでは使いにくい。そのため、論文では対称性のあるカルバック・ライブラー情報量を使っています*1。
単純にひっくり返した距離も足し合わせているだけですね。これで対称性が保たれます。カルバック・ライブラー情報量は、Rには関数が用意されているのですが、SciPy/NumPyにはないので自作します。カルバック・ライブラー情報量の定義式において定数倍や定数項は、距離の比較では必要ないので省略しました。また、行列式の項は、対称性のあるカルバック・ライブラー情報量を使うと打ち消し合って消えるので省略しました。trは、traceの略で行列の対角成分の和です。スクリプトでは、diag()で対角成分を取り出して、sum()で合計しています。
import numpy as np import numpy.linalg def KLDiv(mu1, S1, mu2, S2): """正規分布間のカルバック・ライブラー情報量""" # 逆行列を計算 try: invS1 = np.linalg.inv(S1) except numpy.linalg.linalg.LinAlgError: raise; try: invS2 = np.linalg.inv(S2) except numpy.linalg.linalg.LinAlgError: raise; # KL Divergenceを計算 t1 = np.sum(np.diag(np.dot(invS2, S1))) t2 = (mu2 - mu1).transpose() t3 = mu2 - mu1 return t1 + np.dot(np.dot(t2, invS2), t3) def symKLDiv(mu1, S1, mu2, S2): """対称性のあるカルバック・ライブラー情報量""" return 0.5 * (KLDiv(mu1, S1, mu2, S2) + KLDiv(mu2, S2, mu1, S1))
これで正規分布間の距離が定義できたのでEMDを計算する関数も作ります。Earth Mover's Distance（2012/8/4）にも書きましたが、RにはEMDの輸送問題を解く関数があるのでrpy2を使ってRの関数をPythonから呼び出しました。
import rpy2.robjects as robjects # Rで輸送問題を解くライブラリ # Rのデフォルトパッケージではないのでインストールが必要 # Rでinstall.packages("lpSolve") robjects.r['library']('lpSolve') transport = robjects.r['lp.transport'] def calcEMD(sigFile1, sigFile2): # シグネチャをロード sig1 = loadSignature(sigFile1) sig2 = loadSignature(sigFile2) # 距離行列を計算 numFeatures = sig1.shape[0] # クラスタの数 dist = np.zeros(numFeatures * numFeatures) # 距離行列（フラット形式） for i in range(numFeatures): mu1 = sig1[i, 1:21].reshape(20, 1) # 縦ベクトル S1 = sig1[i, 21:421].reshape(20, 20) for j in range(numFeatures): mu2 = sig2[j, 1:21].reshape(20, 1) S2 = sig2[j, 21:421].reshape(20, 20) # 特徴量iと特徴量j間のKLダイバージェンスを計算 dist[i * numFeatures + j] = symKLDiv(mu1, S1, mu2, S2) # シグネチャの重み（0列目）を取得 w1 = sig1[:,0] w2 = sig2[:,0] # 重みと距離行列からEMDを計算 # transport()の引数を用意 costs = robjects.r['matrix'](robjects.FloatVector(dist), nrow=len(w1), ncol=len(w2), byrow=True) row_signs = ["<"] * len(w1) row_rhs = robjects.FloatVector(w1) col_signs = [">"] * len(w2) col_rhs = robjects.FloatVector(w2) t = transport(costs, "min", row_signs, row_rhs, col_signs, col_rhs) flow = t.rx2('solution') dist = dist.reshape(len(w1), len(w2)) flow = np.array(flow) work = np.sum(flow * dist) emd = work / np.sum(flow) return emd
これでようやく二つの楽曲間の距離を測れるようになりました。次は、与えたクエリと距離が近い楽曲を検索する部分を作ります。
類似楽曲検索
今回は検索対象の楽曲が500曲程度なので単純な線形探索を使いました。クエリの楽曲と500曲全部の間でEMDを計算して、距離が近い順にランキングを出力します。ついでにeyeD3というMP3ファイルからタグ情報を読み書きするライブラリを用いて、MP3ファイルからアーティスト名を読み取って出力してみました。また、検索結果をテキストで出力しただけでは検証しにくいので、オーディオタグを埋め込んだHTMLも合わせて出力しています。
先ほどの3つの関数（KLDiv、symKLDiv、calcEMD）は省略しているので下のスクリプトに追加してください。
#coding:utf-8 import os import sys import numpy as np import numpy.linalg import rpy2.robjects as robjects from collections import defaultdict # mir.py # usage: python mir.py [sig file] [sig dir] [html file] # sig file : クエリ楽曲のシグネチャファイル # sig dir : 検索対象のシグネチャファイルのディレクトリ # html file : 検索結果を出力するHTMLファイル # 引数で指定したシグネチャファイルに近い # 上位N件の楽曲を出力する def loadSignature(sigFile): """シグネチャファイルをロード""" mat = [] fp = open(sigFile, "r") for line in fp: line = line.rstrip() mat.append([float(x) for x in line.split()]) fp.close() return np.array(mat) def getArtist(mp3Path): """MP3ファイルからアーティストを取得""" import eyeD3 try: tag = eyeD3.Tag() tag.link(mp3Path) artist = tag.getArtist() except: artist = "None" # 空白のとき if artist == "": artist = "None" return artist def makeHTML(ranking, htmlFile, N=10): """ランキングをHTML形式で出力""" import codecs fout = codecs.open(htmlFile, "w", "utf-8") # HTMLヘッダを出力 fout.write('<!DOCTYPE html>\n') fout.write('<html lang="ja">\n') fout.write('<head><meta charset="UTF-8" /><title>%s</title></head>\n' % htmlFile) fout.write('<body>\n') fout.write('<table border="1">\n') fout.write(u'<thead><tr><th>ランク</th><th>EMD</th><th>タイトル</th>') fout.write(u'<th>アーティスト</th><th>音声</th></tr></thead>\n') fout.write(u'<tbody>\n') # ランキングを出力 rank = 1 for sigFile, emd in sorted(ranking.items(), key=lambda x:x[1], reverse=False)[:N]: prefix = sigFile.replace(".sig", "") # rawをwavに変換（HTMLプレーヤー用） rawPath = os.path.join("raw", prefix + ".raw") wavPath = os.path.join("wav", prefix + ".wav") if not os.path.exists("wav"): os.mkdir("wav") os.system('sox -r 16000 -e signed-integer -b 16 "%s" "%s"' % (rawPath, wavPath)) # アーティスト名を取得 mp3Path = os.path.join("mp3", prefix + ".mp3") artist = getArtist(mp3Path) # HTML出力 # HTML5のオーディオプレーヤーを埋め込む audio = '<audio src="%s" controls>' % wavPath fout.write("<tr><td>%d</td><td>%.2f</td><td>%s</td><td>%s</td><td>%s</td></tr>\n" % (rank, emd, prefix, artist, audio)) rank += 1 fout.write("</tbody>\n"); fout.write("</table>\n") fout.write("</body>\n") fout.write("</html>\n") fout.close() if __name__ == "__main__": if len(sys.argv) != 4: print "python mir.py [sig file] [sig dir] [html file]" sys.exit() targetSigPath = sys.argv[1] sigDir = sys.argv[2] htmlFile = sys.argv[3] ranking = defaultdict(float) # 全楽曲との間で距離を求める for sigFile in os.listdir(sigDir): sigPath = os.path.join(sigDir, sigFile) emd = calcEMD(targetSigPath, sigPath) if emd < 0: continue ranking[sigFile] = emd # ランキングをEMDの降順にソートして出力 N = 10 rank = 1 for sigFile, emd in sorted(ranking.items(), key=lambda x:x[1], reverse=False)[:N]: print "%d\t%.2f\t%s" % (rank, emd, sigFile) rank += 1 # EMDの昇順に上位10件をHTMLにして出力 makeHTML(ranking, htmlFile, N)
実験
では、いくつかの楽曲をクエリとして与えてさっそく結果をみてみましょう！
この実験は、やはり音が聞けないと大部分の人にはわからないと思うので、各楽曲の中心30秒間のWAVファイルを再生できるHTMLページも用意しました。引用の範囲なのでおそらく大丈夫だとは思いますが・・・著作権者様からの連絡があったら削除します。結果の2列目がクエリとの間のEMDです。1位がクエリの曲になります。1位にクエリと同じ曲がくるのはまあ当然ですよね。
python mir.py sig/愛は勝つ.sig sig ai_ha_katu.html 1 20.00 愛は勝つ.sig 2 51.51 WE GOTTA POWER.sig 3 51.67 時の河.sig 4 52.30 Ilusion.sig 5 52.51 君がいるだけで.sig 6 53.26 君がいるから.sig 7 53.50 心のPhotograph.sig 8 53.65 ルネッサンス情熱.sig 9 54.25 嘘.sig 10 55.17 謳う丘.sig
愛は勝つ（KAN）の検索結果です。小学生のときすごい流行りました（笑）検索結果は、男性の力強いボーカル中心ということでけっこう似ているように思うのですがいかがでしょう？ところどころ女性ボーカルが入ってますが何でだろう？
python mir.py sig/魔法の人.sig sig maho_no_hito.html 1 20.00 魔法の人.sig 2 64.30 変わらないもの.sig 3 65.87 優しい風.sig 4 66.34 Thanks.sig 5 67.06 旅の途中.sig 6 68.11 碧いうさぎ.sig 7 68.32 窓絵.sig 8 69.03 Garnet.sig 9 72.18 奇跡の海 (オリジナル・カラオケ).sig 10 72.19 甘えんぼ.sig
魔法の人（奥華子）の検索結果です。女性の優しいボーカルが検索結果として並んでいるように思います。「変わらないもの」と「Garnet」はクエリと同じ奥華子さんの曲みたいですね。聞いてみたところ雰囲気は何となく似ていると思います。
python mir.py sig/Prayer.sig sig prayer.html 1 20.00 Prayer.sig 2 81.31 Windancer.sig 3 85.38 Escape.sig 4 89.42 Appassionata.sig 5 89.42 First Day Of Spring.sig 6 90.56 Celebration.sig 7 91.46 Steps.sig 8 91.84 Moving.sig 9 97.62 Passacaglia.sig 10 98.97 陽だまりの歌.sig
Prayer（Secret Garden）の検索結果です。海外アーティストの英語の曲です。9位まではすべて同じSecret Gardenのアルバムからでした。同じアルバムということで非常に雰囲気が似ています。10位も優しい曲で同じような雰囲気を感じます。
python mir.py sig/君をのせて.sig sig kimi_wo_nosete.html 1 20.00 君をのせて.sig 2 80.20 魔法のぬくもり.sig 3 82.78 青空っていいな.sig 4 83.66 風のとおり道.sig 5 84.07 カントリー・ロード.sig 6 87.44 Garnet.sig 7 88.42 やさしさに包まれたなら.sig 8 88.91 わたしが不思議.sig 9 89.48 Hello My Friend.sig 10 89.95 つないだ手.sig
君をのせて（井上あずみ）の検索結果です。ジブリの「天空の城ラピュタ」のエンディング曲ですね。2位の「魔法のぬくもり」は同じ井上あずみさんの曲です。「風のとおり道」「カントリー・ロード」「やさしさに包まれたなら」は全部ジブリの曲じゃないか！？これはおどろき。やっぱり雰囲気は似ているのかな？
python mir.py sig/蒼い鳥.sig sig aoi_tori.html 1 20.00 蒼い鳥.sig 2 66.21 蒼い鳥（アレンジ版）.sig 3 83.96 青空っていいな.sig 4 85.03 春よ来い.sig 5 87.98 わたしが不思議.sig 6 88.12 遠い音楽.sig 7 90.25 Appassionata.sig 8 92.86 Dear You -Feel-.sig 9 93.21 やさしさに包まれたなら.sig 10 94.33 葬列.sig
蒼い鳥（今井麻美）の検索結果です。悲しげな感じがするきれいな曲です。2位に同じ「蒼い鳥」という曲が入りましたが、微妙に違うアレンジ版のようです。6位の「遠い音楽」も同じ今井さんの曲です。その他も何か悲しげな曲が多い感じです。
python mir.py sig/さくらんぼ.sig sig sakuranbo.html 1 20.00 さくらんぼ.sig 2 49.84 WILL.sig 3 50.58 Good Day.sig 4 51.12 扉の向こうへ.sig 5 51.75 夢想歌.sig 6 51.80 ショットガン・ラヴァーズ.sig 7 52.08 深く眠れ.sig 8 52.25 Vermillion.sig 9 52.46 西のそらへ.sig 10 52.87 大爆発No.1.sig
最後は、さくらんぼ（大塚愛）の検索結果です。元気な感じの女性のボーカルです。まあ検索結果も元気な感じの女性の曲が並んでいました。6位の「ショットガン・ラヴァーズ」は今話題の初音ミクの曲のようです。合成音声でも関係ないのかな？
おわりに
今回は、類似楽曲検索システムを実験してみました。今回作ったのは非常に基本的なシステムでMFCCという特徴量だけに基づいています。さらにメロディ、リズムの類似性を考慮するなど改良できる点はたくさんあると思います。また、大量の楽曲を対象に検索するには並列処理やより効率的なアルゴリズムも必要ですね。たぶん、実用化に関してはそっちの方が問題になりそうです・・・上であげたPampalkさんの論文はさらに追究するスタートポイントとしてよくまとまっていると思います。
今年の目標（2012/1/1）は、類似楽曲検索システムの実験をすることだったのでとりあえず達成！達成できたのは珍しいぞ。
*1:ただ、先の論文には誤りがあるようです。traceがぬけているため結果がスカラーになりません。
