本記事は、弊社AIコンサルティンググループが運営している「note」の「AIグループ@グランバレイ」にて、2019年8月に掲載した記事を元に一部修正を加えた再掲載いたしました。

はじめに

最近、音声を使って何かやってみたい（できれば音声抽出とか）と考えているのですが、知識やスキルが全然足りず苦戦しております。大学時代は音響系の研究をしてたはずなのですが、ほとんどを忘れてしまい、一から勉強中です。（そもそも大学時代に理解できていたのか怪しい）

今回、音声処理の勉強の一環として、Pythonでwavファイルを読み込み、その周波数スペクトルを描画ライブラリである「Matplotlib」を利用しモニタリングできるものを作ってみましたので、それをご紹介いたします。
シンプルに申しますと、イコライザーなどで表示されるバー？のようなものをPythonで作ってみます。

想定される方

・音声解析に興味のある人
・Pythonでwavファイルの取り扱い方に興味のある人
・Matplotlibのグラフやアニメーションの描画に興味のある人

最終形

まず、作ってみたものは以下のようなアニメーションです。

表示している周波数帯は8つです。本来は各周波数帯でフィルタリングする（バンドパスフィルターをかける）のが正しいかもしれませんが、今回は、各周波数帯における音圧の平均で可視化しております。

また、音圧が高ければ赤色、低ければ青色となるようにしています。将来的には、以下のように可視化するグラフも考えています。

実行環境

OS：Windows10
Python：3.7.3
numpy：1.16.2
Matplotlib：3.0.3

プログラムコード

プログラムコードは太字になっております。

ライブラリのインポート部

# wavファイルを扱うためのモジュール
import wave

import math
import numpy as np

import matplotlib.pyplot as plt
import matplotlib.animation as animation
%matplotlib nbagg

読み込むwavファイルの指定

file = ‘hoge.wav’

スペクトルを取得する間隔やハミング窓のサイズ（秒数）を指定します。今回は0.25秒を設定しています。

スペクトルに関するプログラムコード

# 何秒間隔でスペクトルを取るか
STEP_TIME = 0.25
# ハミング窓サイズ
WINDOW_TIME = STEP_TIME

スペクトルを取得する関数を定義します。スペクトルの取得はnumpyのrfftを使用します。これは、numpyの高速フーリエ変換を行うための関数の一つです。

# スペクトルを取得する関数
def spectrum(frames, fs):
f = np.fft.rfft(frames)
# 音圧
amp = np.abs(f)

# 周波数
freq = np.fft.rfftfreq(frames.shape[0], d=1.0/fs)

return freq, amp

次に、先ほど定義した関数（spectrum）を使い、wavファイルをスペクトルに変換していきます。Pythonでwavファイル扱うために、waveという便利なモジュールがあったので利用してみました。

wave.open()でwavファイルをオープンすると、wavファイルの様々な情報が取得できます。また音声データ自体は、f.readframes(nframes)で取得しますが、このデータはbytes型となっているため、numpyのfrombufferでnumpy配列に変換しています。今回、入力としたwavファイルは量子化ビット数が16のため、dtypeはint16を指定していますが、量子化ビット数が32の場合はint32を指定してください。（[1]の箇所）また、今回入力としたwavファイルはステレオとなっているため、左右それぞれのチェンネルにデータを分離させる必要があります。左右のデータはそれぞれ[左, 右, 左, 右…]のように交互に格納されているため、Pythonのスライスを活用して取得しています。（[2]の箇所）左右それぞれのデータに分離ができたら、ハミング窓を掛け合わせ、最後に先ほど定義したspectrum関数を使って周波数（…_X）と音圧（…_Y）のスペクトルの情報に変換します。

# 左チャンネルのスペクトルを格納
left_spectrum = [] # 右チャンネルのスペクトルを格納
right_spectrum = []

# wavファイルの読み込み
with wave.open(file, ‘rb’) as f:
# チャンネル数
nchannels = f.getnchannels()
# 1サンプルのバイト数
sampwidth = f.getsampwidth()
# サンプリング周波数
framerate = f.getframerate()
# フレーム数
nframes = f.getnframes()
# 圧縮タイプ（wavはNone）
comptype = f.getcomptype()
# 圧縮名（wavはNone）
compname = f.getcompname()
# 全フレームを読み込む
frames = f.readframes(nframes)

# Bytesをndarrayに変換（量子化ビット数16を想定）・・・[1] np_frames = np.frombuffer(frames, dtype= “int16”)

# ステップするフレーム数
step_nframes = int(framerate * STEP_TIME * 2)
# ハミング窓のフレーム数
window_nframes = int(framerate * WINDOW_TIME * 2)

# フレームのポジション
pos = 0

# ハミング窓
window = np.hamming(framerate * WINDOW_TIME)

while len(np_frames) > (pos + window_nframes):
# ステレオデータを左右に分離・・・[2] left_frames = np_frames[pos:pos+window_nframes:2] right_frames = np_frames[pos+1:pos+window_nframes:2]

# 左右それぞれのフレームにハミング窓をかける
window_left_frames = window * left_frames
window_right_frames = window * right_frames

# 左右それぞれのフレームのスペクトルを取得
left_X, left_Y = spectrum(window_left_frames, framerate)
right_X, right_Y = spectrum(window_right_frames, framerate)

# リストに格納する
left_spectrum.append((left_X, left_Y))
right_spectrum.append((right_X, right_Y))

# 次のポジション
pos += step_nframes

X軸（周波数帯）とY軸（音圧）の設定

今回は64Hz～16384Hzのうちの8帯域とします。

#プロットする周波数
Freqs = [64, 128, 256, 1024, 2048, 4096, 8192, 16384] # X軸（[0, 1, 2, 3, 4, 5, 6, 7]）
X = list(range(len(Freqs)))

Y軸（音圧）を取得する関数を定義します。0Hz～64Hzの平均、65Hz～128Hzの平均、129Hz～256Hzの平均、というように各帯域の音圧の平均をとっています。

# Y軸を作成する関数（指定された周波数帯域毎の平均をとる）
def create_Y(spectrum_X, spectrum_Y, Freqs):
Y = [] counter = 0 # Freqs用
tmp = [] for i in range(len(spectrum_X)):
# tmpに音圧を格納
tmp.append(spectrum_Y[i])
if spectrum_X[i] == Freqs[counter]: # Freqsの周波数と一致した場合
avg = sum(tmp) / len(tmp) # 平均をとる
Y.append(avg) # 平均をYに格納する
tmp = [] # tmpを空にする
counter += 1 # カウンタをインクリメント
if counter + 1 > len(Freqs): # カウンタがFreqsの長さを超えたらbreak
break

return Y

グラフにアニメーション描画

plotという関数は、アニメーションのフレーム単位で呼ばれる関数です。基本的にplot関数の中でグラフを作ると、それがパラパラ漫画のようにアニメーション化されます。Matplotlibでアニメーションを作るときの注意点ですが、cla関数というグラフをクリアする関数を最初に呼び出さないと、前に描画したグラフの上に新しいグラフが重なって描画されてしまいます。（[1]の箇所）

また、今回実装した音圧の大きさによって棒グラフの色が変わるという部分ですが、これはfor文でbar関数を1回1回呼び出して、それ毎に色（color引数）を指定することで実現しています。指定する色はカラーマップから取得しています。（[2]の箇所）なお、カラーマップは、plt.get_cmap(カラーマップ名)でMatplotlibに定義されているカラーマップ（cmap）を取得しています。このカラーマップは0.0～1.0の値を入れてやると、その値に応じた色が取り出せるのですが、今回は0.0だと青色、1.0だと赤色になる’jet’というカラーマップを使用しました。（[3]の箇所）

最後にちょっとした注意点ですが、縦軸の音圧は対数にすることです。（[4]の箇所）私が学生時代に聞いた話では、人間の感覚は対数で感じるとか。周波数の話になりますが、例えばラの音（A3）は440Hzであることが有名ですが、1オクターブ下のラの音（A2）は220Hzで、1オクターブ上のラと音（A4）は880Hzになり、2を底とした対数となります。このように、音圧に関しても、対数で表現したほうがより人間の感覚に近いイメージになります。（間違ってたらすみません。）

fig, (L, R) = plt.subplots(nrows=2, figsize=(8, 12))

# カラーマップ・・・[3] cmap = plt.get_cmap(‘jet’)

max_Y = 10**7

# グラフプロットするための関数
def plot(i):
index = i % len(left_spectrum)
##### 左チャンネルのグラフ #####
L.cla() # グラフをクリア・・・[1] L_Y = create_Y(left_spectrum[index][0], left_spectrum[index][1], Freqs)
# 棒グラフを一つ一つ設定する。・・・[2] for x, y in zip(X, L_Y):
# Y値によって色を変える
L.bar(x, y, color=cmap(math.log(y, 10)/math.log(max_Y, 10)))
L.set_title(‘Left’)
L.set_xlabel(‘Hz’)
L.set_xticks(X)
L.set_xticklabels(Freqs)
L.set_yscale(‘log’, basey=10) # 対数で表示・・・[4] L.set_ylim(0, max_Y)

##### 右チャンネルのグラフ #####
R.cla() # グラフをクリア・・・[1] R_Y = create_Y(right_spectrum[index][0], right_spectrum[index][1], Freqs)
# 棒グラフを一つ一つ設定する。・・・[2] for x, y in zip(X, R_Y):
# Y値によって色を変える
R.bar(x, y, color=cmap(math.log(y, 10)/math.log(max_Y, 10)))
R.set_title(‘Right’)
R.set_xlabel(‘Hz’)
R.set_xticks(X)
R.set_xticklabels(Freqs)
R.set_yscale(‘log’, basey=10) # 対数で表示・・・[4] R.set_ylim(0, max_Y)

# アニメーションの実行
ani = animation.FuncAnimation(fig, plot, interval=int(STEP_TIME * 1000))
# アニメーションをGIFアニメで保存する場合は以下。（ImageMagickが必要）
# ani.save(‘hoge.gif’, writer=”imagemagick”)

最後のコードをコメントアウトしてますが、これはアニメーションを保存するコードです。GIFアニメで保存する場合はImageMagickを入れる必要があるのでそこだけ注意してください。（MP4でも保存できますが、MP4の場合はffmpegが必要となります）

以上が今回のプログラムソースとなります。

まとめ

今回このプログラムを作成したことで、wavファイルの扱い方やFFT（高速フーリエ変換）の使い方、Matplotlibでのアニメーションの作り方のノウハウを得ることができました。

難しく感じますが、Pythonを使えば容易に可視化できることがわかると思います。

最後にお知らせになりますが、このようなPythonのプログラムを覚えたい方向けに、グランバレイでは、「人工知能・機械学習講座（初級編-Python機械学習入門コース）」を開催中です。初歩からお教えしますので、是非ともご検討ください。

---

本記事は、弊社AIコンサルティンググループが運営している「note」内の「AIグループ@グランバレイ」の記事を一部修正を加え転載しております。
https://note.com/gvaiblog/n/n2e426f233003

※ その他の会社名、製品名は各社の登録商標または商標です。
※ 記事の内容は記事公開時点での情報です。閲覧頂いた時点では異なる可能性がございます。