pythonで簡単な音声認識をやってみたいぞ。
そもそも何から始めればいいのかしら。
今回は，基本的な音声認識をpythonで行う方法をお伝えしていこうと思います。本記事はpython実践講座シリーズの内容になります。その他の記事は，こちらの「Python入門講座/実践講座まとめ」をご覧ください。
お題
STFTをpythonで実装してみよう。
流れ
今回対象とするデータは「あいうえお」ではなく，データ長が2をべき乗に制限した「ドレミファソラシド」にしようと思います。以下のファイルを対象とします。
短時間フーリエ変換の概要
STFTの簡単な概要は以下のページで説明しています。ここでは，実装方法をメインにお伝えしていこうと思います。
実装
前提として，窓関数は既成のライブラリを使用します。また，FFTの実装はこちらのページで説明しているものを利用します。また，wavファイルの読み込みと書き出しにはlibrosaを利用します。
必要なライブラリのインポート
import numpy as np
import matplotlib.pyplot as plt
import librosa
from scipy import hammingFFTの定義
import numpy as np
def FFT(x):
N = x.shape[0]
# 再起動作の一番最後
if N==1:
return x[0]
x_even = x[0:N:2]
x_odd = x[1:N:2]
# ここで再帰動作をする
X_even = FFT(x_even)
X_odd = FFT(x_odd)
# DFTと同じようにWを求める
W = []
for t in range(N//2):
W.append(np.exp(-1j * ((2*np.pi*t) / N)))
W = np.array(W)
# ここで型をcomplexに指定しないとエラーを吐くので注意
X = np.zeros(N, dtype="complex")
X[0:N//2] = X_even + W*X_odd
X[N//2:N] = X_even - W*X_odd
return XSTFTの定義
librosaの公式ドキュメントにもある通り，librosaのSTFTでは窓関数を端っこに適用するときは反転パディングを行っています。今回もそれに習いましょう。また，実信号に対するFFTの結果が対象であることに注意して，後ろ半分の情報を捨ててしまうように実装しています。
def STFT(x, win_length, hop=0.5):
'''
今回はwin=Mとして実装。
データサイズも簡単のため2のべき乗に制限。
hopはデフォルトで窓幅の半分とする。
'''
hop_length = int(win_length * hop)
# 窓関数をかける時に端点が問題になります。
# 今回はlibrosaのデフォルトに習って反転パディングをしてみます。
pad_first = x[:hop_length]
pad_last = x[-hop_length:][::-1]
x_pad = np.concatenate([pad_first, x, pad_last])
# データのサンプル数
N = x_pad.shape[0]
# FFTの結果を半分にした長さ
M = int(win_length//2) + 1
# 窓関数を適用する回数
T = int((N - hop_length)/hop_length)
# 今回はハミング窓（ハン窓）を利用します
han = hamming(win_length)
# 結果を格納する箱です
spec = np.zeros((M,T), dtype="complex")
for t in range(T):
# まず窓関数を適用します
windowed_x = x_pad[t*hop_length:t*hop_length+win_length] * han
# 次にFFTを実行します。
spec[:,t] = FFT(windowed_x)[:int(win_length//2)+1]
return specデータの読み込み
# [path_to_wavfile]は皆さんの環境に合わせて変更してください。
y, sr = librosa.load("[path_to_wavfile]")STFTの実行
スペクトログラムの描画にはlibrosaを使います。また，こちらのページでdb（デシベル）への変換は説明していますので，今回は振幅からdbへの変換はlibrosaに頼ってしまいます。
spec = STFT(y, 1024, 0.5)
spec_db = librosa.amplitude_to_db(np.abs(spec))
librosa.display.specshow(spec_db, y_axis="log")
逆変換
短時間フーリエ変換の逆変換を実装するのも，よい練習になります。今回私も実装してみたのですが，復元された「ドレミファソラシド」の音質が非常に悪くなってしまいました。あとでlibrosaの逆変換と比べてみます。まずは，実装からみてみます。
実装
高速フーリエ変換の逆変換はこちらのページで説明しています。そのままのコードを記しておきます。
def IFFT(x):
x = x.conjugate()
X = FFT(x)
return (1/N) * X.conjugate()こちらのIFFTを利用して，本当に文字通りSTFTとは反対の操作を施していきます。
def ISTFT(X, win_length, hop=0.5):
# スペクトログラムの横（時間方向）の長さ
T = X.shape[1]
# 窓関数のホップ長
hop_length = int(win_length * hop)
# 今回もハン窓を利用します
han = hamming(win_length)
# STFTでFFTの結果の半分を捨てていましたので対称性に気をつけながら復元します。
X_recover = np.concatenate([X[:-1], X[::-1, :][1:]])
# 改めて復元するデータの長さを計算します
# STFTでTを計算した式の逆を求めているだけです
N = int((T * hop_length) + hop_length) - win_length
# 改めて窓関数を適用する回数を計算します
nm_hop = int((N - hop_length)/hop_length)
# データを格納する箱です
x = np.zeros(N, dtype="complex")
# 端点だけ場合分けしてIFFTを施していきます
for i in range(nm_hop):
if i==0:
x[i*hop_length:i*hop_length + int(win_length//2)] += IFFT(X_recover[:,i])*han[int(win_length//2):]
if i==X_recover.shape[0]-1:
x[i*hop_length:i*hop_length + int(win_length//2)] += IFFT(X_recover[:,i])*han[:int(win_length//2)+1]
else:
x[i*hop_length:i*hop_length + win_length] += IFFT(X_recover[:,i])*han
return x結果
まず，librosaで復元した音源を聞いてみましょう。
spec = librosa.stft(y)
y_recover_librosa = librosa.istft(spec)
librosa.output.write_wav("test_librosa.wav", y_recover_librosa, sr)次に，自作関数による復元結果を聞いてみましょう。最後に「np.asfortranarray」でfortranarrayに変換しておかないとエラーを吐くので注意です。
x_recover = np.real(ISTFT(spec, 1024, 0.5))
librosa.output.write_wav('test.wav', np.asfortranarray(x_recover), sr)同じディジタル信号処理のはずなのに，なぜここまで差が出るのかについては要検討です。
