pythonのライブラリlibrosaのクロマグラム(chromagram)を用いた和音認識

[@namaozi](https://twitter.com/namaozi)です。
この記事は[eeic Advent Calendar 2016](http://qiita.com/advent-calendar/2016/eeic) 9日目の記事 **その2** です。

いきなり **その2** から見始めた方はよろしければ [**その1**](https://qiita.com/namaozi/items/dec1575cd455c746f597) を先にお読みください。

前回はスペクトログラムを作って音楽信号を可視化しました。

今回はクロマグラムを用いて

[このサンプル音源](http://twitsound.jp/musics/ts14XqIQU)

の和音認識にチャレンジしてみましょう。


# この記事でやりたいこと

* クロマグラムを使って和音認識してみる
    * なるべく数学・プログラミング分からない人でも読めるようにしたい
    * できるだけ定性的で音楽的な話を入れていきたい


## おことわり

###ライブラリ、使います。
[その1]()ではスペクトログラムをnumpyだけを使って実装したのですが、クロマグラムの和音認識は自分の実装だと<font color="blue">精度がゴミofゴミすぎたのでライブラリ様に頼ることにします</font> 、情けない話だ…(´･_･｀)

…というわけでこの記事では高度なpythonの音声信号解析ライブラリの[**librosa**](https://github.com/librosa/librosa)を使ってクロマグラムを取得し、和音認識してみましょう。実際にやってみたい方はまずリンクからlibrosaをインストールしてください。

###python2.7を使ってます
大変申し訳ありませんがpython2.7をこのサンプルでは使ってます、、、じきに3.0対応したいと思います

###言葉の問題
音楽用語とかいろいろ出てきますが、分からないところは基本スルーして大丈夫だと思います。興味のある方は、リンクを貼るようにしますのでそちらを参照ください。
とりあえず**コードは和音で、名前がついた和音がコード**ということさえわかっていればOKです。


## 和音(コード)認識とは

**コンピューターにコード進行を耳コピさせるタスク**です。

実は和音を推定するタスクは非常に難しい問題で、現在でも盛んに研究が続けられています。年々精度が改善されているものの、いまだに訓練された人間のように柔軟に和音を聴き分ける段階には至っていません。

詳しく知りたい方は[この辺の記事](http://hil.t.u-tokyo.ac.jp/mklab/publications/Kameoka2015IEICE06published.pdf)を読んでみてください。今回扱うクロマグラムのことも書かれています。最新手法については[こちらのブログ](http://xiao-ming.digick.jp/mir/mirex2016chord/)が詳しいです。

本当は**隠れマルコフモデル**(Hidden Markov Model, HMM)を用いた実装が主流なのですが、実装が大変なので今回は簡易版としてクロマグラムのみを用いてコード認識にチャレンジしたいと思います。

## クロマグラムとは

[クロマグラム](https://en.wikipedia.org/wiki/Chroma_feature)は、全周波数帯域のパワーを[C, Db, D, Eb, E, F, Gb, G, Ab, A, Bb, B]の12音階に落とし込み、ある区間の時間における音の成分を可視化したものです。
英語のwikiの画像で説明します。

![chromagram](https://upload.wikimedia.org/wikipedia/commons/2/25/ChromaFeatureCmajorScaleScoreAudioColor.png)

(a)はCメジャースケール、いわゆるピアノの白鍵の音が並んでいる楽譜です。
(b)が(a)に対応する理論的なクロマベクトルになるんですが、実際にCメジャースケール(ハ長調)を弾いてみた音の波形である(c)からクロマグラムを作成すると(d)のようになります。
イメージとしては**音楽信号をピアノロール1オクターブ分に全て押し込んだ感じ**でしょうか。

各時刻における信号は12音の成分の強さで表せるので、**クロマベクトル**と言ったりもします。

## 結局どうやるの？

どうやってコード認識するのかという話ですが、

**"音楽信号のクロマグラムがどれくらい欲しいコードのクロマベクトルと似ているのか"**

によって各時刻のコードを決めていきます。具体的には以下のような手順になります。

**1.** 認識したいコードのクロマベクトル(テンプレートベクトル)を作る

**2.** 楽曲のクロマグラムを作る

**3.** 各時刻のクロマグラムとテンプレートベクトルの類似度を計算して、最も類似度が高いコードをその時刻のコードと推定する

## 1.テンプレートベクトルを作る

まず、認識したいコード(和音)の理想的なクロマベクトルを作っておきます。以降ではこうした理想的なクロマベクトルを**テンプレートベクトル**と呼ぶことにします。

たとえば[ダイアトニックコード](http://www.chaco2008.com/theory/1-8_diatonicchord.html)は下の図のような感じのテンプレートベクトルになりますね。ちなみにダイアトニックコードというのは、ある調の構成音を3度ずつ、つまり1つ飛ばしで堆積させた和音です(長調・短調それぞれにダイアトニックコードがありますが下の画像は長調です)。

![diatonic_chroma2.png](https://qiita-image-store.s3.amazonaws.com/0/103754/8ee5010f-cc0b-3ad4-575f-1e567511aaa5.png)

ちなみにさっきは黒かったのに上の画像のクロマベクトルの色が違っているのは、ある時刻における各ピッチクラスのパワーの和を1にそろえるためです。

上の画像の一番下に書いてある数式がCメジャーコードのテンプレートベクトルを数式で表したものになります。各次元が12種類の[ピッチクラス](https://ja.wikipedia.org/wiki/%E3%83%94%E3%83%83%E3%83%81%E3%82%AF%E3%83%A9%E3%82%B9)に対応していて、Cから順に半音ずつ上がっています。

## 2.クロマグラムを作る

先ほど紹介したクロマグラムを作ります。神ライブラリのlibrosa様に頼りましょう。ついでなのでパーカッシブ成分を分離したものを解析しています。内部で何をしているか気になる方は[神ドキュメント](http://librosa.github.io/librosa/generated/librosa.feature.chroma_cens.html?highlight=chroma_cens#librosa.feature.chroma_cens)を参照ください。

コードはGitHubに置いてありますのでよかったら試してみてください。
https://github.com/namaozi/ChromagramSample

```python:functions.py

def librosa_chroma(file_path="audios/harmony1.wav", sr=44100):
    #インポート(インストールしないと使えません)
    import librosa

    # 読み込み(sr:サンプリングレート)
    y, sr = librosa.load(file_path, sr=sr)

    # 楽音成分とパーカッシブ成分に分けます
    y_harmonic, y_percussive = librosa.effects.hpss(y)

    # クロマグラムを計算します
    C = librosa.feature.chroma_cens(y=y_harmonic, sr=sr)

    # プロットします
    plt.figure(figsize=(12,4))
    librosa.display.specshow(C, sr=sr, x_axis='time', y_axis='chroma', vmin=0, vmax=1)
    plt.title('Chromagram')
    plt.colorbar()
    plt.tight_layout()
    plt.show()

    return C
###
```
なんと数行でクロマグラムが計算できちゃいます。はぁライブラリ神…
以下が得られたクロマグラムです。
![harmony1_chroma.png](https://qiita-image-store.s3.amazonaws.com/0/103754/23dd9be2-2071-f8fe-8a24-bfa7e11bff23.png)

[このサンプル音源](http://twitsound.jp/musics/ts14XqIQU)はBPM120で

![chord_types.png](https://qiita-image-store.s3.amazonaws.com/0/103754/38c33bf5-2a4a-0235-8df5-7336f8523143.png)

というコードを順に4拍ずつ、つまり2秒ずつ伸ばしています。得られたクロマグラムを見ると、2秒ごとのコードチェンジのタイミングで色の濃淡が変わり、信号の音の成分が変わっているのが分かるかと思います。 [^1]

更に、例えば3つめ、4秒~6秒の間のGメジャーコードのクロマグラムは**D, G, B**の部分の色が濃くなっていて、これは**Gメジャーコードの構成音と一致**しています。クロマグラムを見るだけでなんとなくコードが分かりそうですよね！ [^2]

[^1]: 本当はどこからどの区間までがそのコードであるのかも計算で決めなければいけないのですが、今回のサンプル音源では2秒おきにコードが変わっていますので、その情報をプログラムに与えてしまいます。

[^2]: 興味がある方は、audiosフォルダにいくつかサンプル音源を入れていますのでそちらを聴いてみて、クロマグラムを作ってみると面白いかもしれません。harmony1.wavがオルガンのコードとドラムだけだったのに対し、harmony2.wavはベースを加え、harmony3.wavが更にメロディーを加えた音源になっています。実際にクロマグラムを見ると、メロディーやベースの音に対応する部分の色が濃くなって見えるはずです。


## 3. コード推定
あとは時間軸に沿ってコードを推定していくだけです。

```python:chordestimation.py
#coding:utf-8
from collections import OrderedDict
import numpy as np
import matplotlib.pyplot as plt
import librosa #インストールしてください!

import soundanalysis as sa

file_path = "audios/"
file_name = "harmony1.wav"

# クロマグラムを求めます
chroma = sa.librosa_chroma(file_path + file_name)

TONES = 12 # ピッチクラス,音の種類の数
sampling_rate = 44100 #音源依存

# "この設定では",こんな感じで時間軸設定を求められます
# (詳しくはドキュメントを読んで下さい)
time_unit = 512.0 / 44100 # 1フレームのクロマグラムの長さ
# stop = time_unit * (chroma.shape[1] - 1)
stop = time_unit * (chroma.shape[1])
time_ruler = np.arange(0, stop, time_unit)

###コードのテンプレートベクトルです
#メジャーとマイナーだけを考えます
#脳筋コード,時間が無かったので…
#順番を保ちたいのでOrderdDictを使います
one_third = 1.0/3
chord_dic = OrderedDict()
chord_dic["C"] = [one_third, 0,0,0, one_third, 0,0, one_third, 0,0,0,0]
chord_dic["Db"] = [0, one_third, 0,0,0, one_third, 0,0, one_third, 0,0,0]
chord_dic["D"] = [0,0, one_third, 0,0,0, one_third, 0,0, one_third, 0,0]
chord_dic["Eb"] = [0,0,0, one_third, 0,0,0, one_third, 0,0, one_third, 0]
chord_dic["E"] = [0,0,0,0, one_third, 0,0,0, one_third, 0,0, one_third]
chord_dic["F"] = [one_third, 0,0,0,0, one_third, 0,0,0, one_third, 0,0]
chord_dic["Gb"] = [0, one_third, 0,0,0,0, one_third, 0,0,0, one_third, 0]
chord_dic["G"] = [0,0, one_third, 0,0,0,0, one_third, 0,0,0, one_third]
chord_dic["Ab"] = [one_third, 0,0, one_third, 0,0,0,0, one_third, 0,0,0]
chord_dic["A"] = [0, one_third, 0,0, one_third, 0,0,0,0, one_third, 0,0]
chord_dic["Bb"] = [0,0, one_third, 0,0, one_third, 0,0,0,0, one_third, 0]
chord_dic["B"] = [0,0,0, one_third, 0,0, one_third, 0,0,0,0, one_third]
chord_dic["Cm"] = [one_third, 0,0, one_third, 0,0,0, one_third, 0,0,0,0]
chord_dic["Dbm"] = [0, one_third, 0,0, one_third, 0,0,0, one_third, 0,0,0]
chord_dic["Dm"] = [0,0, one_third, 0,0, one_third, 0,0,0, one_third, 0,0]
chord_dic["Ebm"] = [0,0,0, one_third, 0,0, one_third, 0,0,0, one_third, 0]
chord_dic["Em"] = [0,0,0,0, one_third, 0,0, one_third, 0,0,0, one_third]
chord_dic["Fm"] = [one_third, 0,0,0,0, one_third, 0,0, one_third, 0,0,0]
chord_dic["Gbm"] = [0, one_third, 0,0,0,0, one_third, 0,0, one_third, 0,0]
chord_dic["Gm"] = [0,0, one_third, 0,0,0,0, one_third, 0,0, one_third, 0]
chord_dic["Abm"] = [0,0,0, one_third, 0,0,0,0, one_third, 0,0, one_third]
chord_dic["Am"] = [one_third, 0,0,0, one_third, 0,0,0,0, one_third, 0,0]
chord_dic["Bbm"] = [0, one_third, 0,0,0, one_third, 0,0,0,0, one_third, 0]
chord_dic["Bm"] = [0,0, one_third, 0,0,0, one_third, 0,0,0,0, one_third]

prev_chord = 0
sum_chroma = np.zeros(TONES)
estimate_chords = []

result = np.zeros((TONES * 2, 8))

for time_index, time in enumerate(time_ruler):
    # 今は何番目のコードを解析しているのか
    # 2秒おきに変わるので2で割って求めます
    nth_chord = int(time) / 2

    # 次の2秒間に移る時に,前の2秒間のコードを推定します
    if nth_chord != prev_chord:
        maximum = -100000
        this_chord = ""
        # コサイン類似度が最大になるコードを調べます
        for chord_index, (name, vector) in enumerate(chord_dic.iteritems()):
            similarity = sa.cos_sim(sum_chroma, vector)
            result[chord_index][nth_chord - 1] = similarity
            if similarity > maximum:
                maximum = similarity
                this_chord = name
        # 初期化、推定したコードを格納します
        sum_chroma = np.zeros(TONES)
        estimate_chords.append(this_chord)

    else:
        # chromaのshapeに注意しながら足していきます
        for i in range(TONES):
            sum_chroma[i] += chroma[i][time_index]

    # 更新
    prev_chord = nth_chord
###

# 最終結果です
print estimate_chords


###がんばってプロットします
axis_x = np.arange(0, 16, 2)
bar_width = 0.07
colors = ["#ff9999", "#ffaf95","#fabb92","#ffd698","#fae991","#c1fc97","#97fac8","#96f9f5","#98e1fb","#9cb2ff","#b79bfe","#fa96f9", "#b36a6a", "#ab7361","#aa7d61","#ad9165","#b4a765","#8ab66b","#6ab48f","#68b0ad","#689fb3","#6979b0","#7462a3","#aa62a9"]
for i, (name, vector) in enumerate(chord_dic.iteritems()):
    plt.bar(axis_x - ((axis_x[1] - axis_x[0]) * 0.45) + bar_width * i, result[i], color=colors[i], width = bar_width, label = name, align = "center")

plt.legend()
plt.xticks(axis_x + bar_width / 24)
plt.show()

```

上のコードを実行して得られた結果！！！
![chord_result.png](https://qiita-image-store.s3.amazonaws.com/0/103754/71df7d31-fe1a-93de-043f-4bb8ead5893b.png)

というコード列が得られました！！
[入力した曲](http://twitsound.jp/musics/ts14XqIQU)のコードは
![chord_types.png](https://qiita-image-store.s3.amazonaws.com/0/103754/38c33bf5-2a4a-0235-8df5-7336f8523143.png)

だったのですが、今回はメジャーとマイナーコードしか認識できないので、2つめの**Dm7**(Dマイナーセブンス)はそもそも正解が出力できません。Dm7の構成音が[D,F,A,C]であることを考えると,DマイナーかFメジャーであればそれなりに認識出来ていることになります。出力では**Fメジャー**が得られているのでほぼ正解ですね。
7つめのコード**Am/E**も、構成音はAmと同じなので出力のAmで正解です。

こうしてみると、<font color="red">サンプル音源「harmony1.wav」では、クロマグラムとテンプレートベクトルを使って、全てのコードに対してほぼ正解を出力することが出来ました！</font>めでたしめでたし。 [^3]

[^3]: 実世界の曲は、単純なクロマグラムだけからこんなにうまく耳コピができちゃうことは残念ながらまずありません。harmony1.wavは和音+ドラムのみというごく単純な構成でしたが、実世界の曲は無限のトラック数を使った複雑な曲が多いです。それに歌声のメロディーなんかが入ると大きく影響を受けるので、和音を推定することは更に難しくなります。脚注1で述べた、「どの区間がそのコードであるのか」を推定するのも難しいので、余計に和音解析タスクが難しくなります。

各区間の全コードに対する類似度をプロットしたものが以下の図になります(上のプログラムを実行すれば得られます)、とてもきれいな可視化ではないですが、、、
![similarity_distribution.png](https://qiita-image-store.s3.amazonaws.com/0/103754/c9a93651-5741-9c89-44d5-0fc03c0b8825.png)

この図を(がんばって)見ると、例えば最初の0~2秒の区間では**C**の類似度が一番高くなっていますが、そのほかに同主単調の主和音の**Cm**, 平行調の主和音の**Am**, 属調(IIIm)の主和音**Em**などでも類似度が高くなっていることがわかると思います。

<font color="red">このように関係調の主和音の類似度が似通っていることが分かる</font>と思います。このような数理的解釈から音楽理論を再考してみるのも面白いかもしれませんね。




## まとめ
- pythonはべんり
- librosa最強
- 自分で実装する前にライブラリでできないか探しましょう
- コード認識は最強のライブラリを使っても結構難しい