この記事はプリッカソン Advent Calendar 2018 21日目の記事です。
飛び入り参加失礼します。

はじめに

始めまして、四国の山の中で~~冴えない男子工学部生~~女児をしているあいうえと申します。
普段は地球物理っぽい分野の勉強を行っています。
最近は趣味でねっとわーくの勉強を行ったりちょろっとGolangを書いたりしています。

この記事は全くPythonに触れたことのない筆者がPythonを使って音楽データのスペクトログラムを出力してみるという記事です。
Pythonバージョンは3.7.1です。

ぷり関係ないって？
私は「プリパラが好きぷり!」
大丈夫って事にして頂きたい…

やること(やりたかったこと)

適当な音楽データのスペクトログラムを見て眺めるだけでは「あぁ～この波形かわいいなぁ！！！」という事しか出来ないのでWITHの楽曲「Giraギャラティック・タイトロープ」を題材にしてMY☆DREAMバージョンとの比較、相互相関を取ってみれば面白いかなあと思いました。
ですが今回は相互相関を取るところまで実装出来ませんでした。
ただ波形を見て「いいぜ！」「かわいいなあ！」というだけの内容になります。

フーリエ変換

フーリエ変換は以下の式で表されます。

$F(ω) = \mathfrak{F}[f(t)]= \int^{\infty}_{-\infty}f(t)\exp(-j2\pi f)dt$

信号解析を睡眠学習で履修した私にはこの式を見てもちんぷんかんぷんです。

この式が表すことは時間の関数を角周波数の関数に変換出来るということです。
変換すると、ある角周波数の時の振幅は○○だ！！と分解することが可能という事です。
角周波数成分に変換してそれぞれがどれくらいの強さであるか分かると何が嬉しいのかというと、
我々が普段聞いている音はsin波の足し合わせで出来ており、それを数値で表現できるという事です。

スペクトログラム

スペクトログラムは音を可視化したグラフです。
時間、周波数、信号の強さの3次元のグラフで、表されます。時間をX軸に、周波数をY軸に持ち、信号の強さを色の違いで表します。
刑事ドラマなんかでよく見る声紋鑑定はスペクトログラムを用いた鑑定のことです。
今回はmatplotlibの機能を用いてグラフ出力してみました。
モジュールの豊富、Python羨ましいです。

具体的に中で何をしているのか、実装はおそらく短時間フーリエ変換を行っているのだと思います。(モジュールを使うときは中身も確認しましょう。します…)
短時間フーリエ変換はある信号に対してフーリエ変換をする幅を決め窓関数をかける、幅の位置を少しずらして窓関数をかける、という行為を繰り返していくものです。

実装

今回作成したプログラムを以下に示します。

# モジュール読み込み
import wave
import numpy as np
import matplotlib.pyplot as plt

# waveモジュールを用いてwave各要素を読み込み
file_path = "15.wav"
wr = wave.open(file_path,"rb")
ch = wr.getnchannels()
width = wr.getsampwidth()
fr = wr.getframerate()
fn = wr.getnframes()
buf = wr.readframes(wr.getnframes())
wr.close()
length = 1.0 * fn / fr

# データを正規化
data = np.frombuffer(buf, dtype = 'int16')/32768.0

# 左右チャネルに分離
lc = data[::2]
rc = data[1::2]

# 確認用
print('チャンネル', ch)
print('総フレーム数', fn)
print('フレームレート', fr)
print('フレーム長', width)
print('再生時間', length)
print(lc)
print(rc)

# サンプル数と窓関数定義
N = 1024
h_window = np.hamming(N)

# 左チャネルについて
plt.subplot(2, 1, 1)
plt.title("left_channel")
pxx, freqs, bins, im = plt.specgram(lc, NFFT=N, Fs=wr.getframerate(), noverlap=0, window=h_window)
plt.axis([0, length, 0, wr.getframerate() / 2])
plt.xlabel("time [second]")
plt.ylabel("frequency [Hz]")

# 右チャネルについて
plt.subplot(2, 1, 2)
plt.title("right_channel")
pxx, freqs, bins, im = plt.specgram(rc, NFFT=N, Fs=wr.getframerate(), noverlap=0, window=h_window)
plt.axis([0, length, 0, wr.getframerate() / 2])
plt.xlabel("time [second]")
plt.ylabel("frequency [Hz]")

plt.tight_layout()
plt.savefig('figure.png')

やってること

waveモジュールを用いてwave形式な音楽データの各要素の読み込み
データの正規化で取り得る範囲を-1~1に
読み込むデータだ2ch音声と仮定しチャネル分離
サンプル数(幅)と窓関数(今回はハミング窓)の決定
それぞれチャネルについてspecgram関数を用いてプロット
画像を保存

実行した結果が下です。
一回目の実行がWITH、二回目の実行がマイドリです。

[user@PC] % python fft_01.py
チャンネル 2
総フレーム数 11791164
フレームレート 44100
フレーム長 2
再生時間 267.37333333333333
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
/Users/user/.pyenv/versions/3.7.1/lib/python3.7/site-packages/matplotlib/axes/_axes.py:7609: RuntimeWarning: divide by zero encountered in log10
  Z = 10. * np.log10(spec)
[user@PC] % python fft_03.py
チャンネル 2
総フレーム数 11747064
フレームレート 44100
フレーム長 2
再生時間 266.37333333333333
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
/Users/user/.pyenv/versions/3.7.1/lib/python3.7/site-packages/matplotlib/axes/_axes.py:7609: RuntimeWarning: divide by zero encountered in log10
  Z = 10. * np.log10(spec)

なにやら怒られています。
これは0(無音時間)は対数取ると無限大になるよ！と言われていますね。
これでも私の環境ではグラフをプロットすることは可能でした。

結果

出力された画像を以下に示します。
上の画像がWITH、下の画像がマイドリです。

f:id:ittyo3103:20181217203645p:plain f:id:ittyo3103:20181217203657p:plain

WITHの方が全体的に濃く色づいており「いいぜ！」って感じですね。
「できるできる Galaxy」、かわいい。

おわりに

今回の記事ではスペクトログラムを表示させただけとなっています。

Pythonいいですね。文法全く理解して無くてもドキュメント見るだけでなんとなく動きました。
こう書いた方が良いって意見あればコメントください。

スペクトログラムを出すだけであればプログラムを書くまでもなくSoXというソフトを用いれば簡単に出すことが可能です。
より高度？なことをする為にPythonに触れて見ましたが時間がありませんでした。

相互相関を取るとキーの違い(♭5)である点や歌声の違いで面白い結果が出るんだろうなとか思います。

今後の発展としてはボイスのみを抽出し音声の特徴量(メル周波数ケプストラム係数)抽出なんかをするとパート分けを聞き分けたり出来そうですね(やってる人がおられます)
めるめるだったらチョチョイのチョイですね

追伸
ウィンターライブ最高でした！

Pripara friendship tour 2019 大阪昼夜当選してました :iwai:
シルバークラスな若輩者ですがフォロチケ交換して下さい><

blog

Giraギャラティック・タイトロープ可視化やってみた！

はじめに

やること(やりたかったこと)

フーリエ変換

スペクトログラム

実装

結果

おわりに