音声ファイル特徴量変換（その1）STFT

音声データを使う際は特徴量に変換する。
よく使われているのは「MFCC」だが、深層学習では「log-mel spectrogram」を使う実装例も出ている。

特徴量	実装例
STFT	Looking to Listen at the Cocktail Party
log-mel spectrrogram	Deep Learning for Audio Signal Processing
MFCC	HTK, Julius, Kaldi

今回は「STFT」（短時間フーリエ変換）について見てみる。
「TensorFlow API」を使う例はこちら
「log-mel spectrrogram」については、「libROSA」を使う例または「TesorFlow API」を使う例
「MFCC」については、「libROSA」を使う例または「TesorFlow API」を使う例

音声データは前回と同じく、「yes」という一秒間の発話データ。

フーリエ変換（Fourier Transform）

フーリエ変換（正規化無し）

import numpy as np
import librosa

# Audio Data
audio_path = 'speech_dataset/yes/0a7c2a8d_nohash_0.wav'
data, sr = librosa.load(
    audio_path,
    sr=16000)

print(data.shape)  # --> (16000,)
print(sr)  # --> 16000

# 高速フーリエ変換
F = np.fft.fft(data)

データの中身を確認。

print(F[:3])

[ 0.04498291+0.j  -0.02268502+0.00299779j  0.02986105+0.00203091j ]

振幅に変換する。

# 振幅に変換
amp = np.abs(F)

データの中身を確認。

print(amp[:3])

[0.04498291 0.02288224 0.02993003]

3番めのデータ（インデックス「2」）で検算。
「実数の二乗」と「虚数の二乗」を足して平方根をとる。

>>> math.sqrt(( 0.02986105 ** 2 ) + ( 0.00203091 ** 2))
0.02993003345354963

平方根をとる前の状態を「power spectrum」、平方根をとった後の状態を「amplitude spectrum」として区別する。

用語	類似表現	意味	備考
amplitude spectrum	magnitude	振幅
power spectrum	Energy	パワー	振幅の二乗

「スペクトラム」と「スペクトル」は同じ意味。

言語	単語	品詞
英語	spectrum	名詞
	spectral	形容詞
フランス語	spectre	名詞

プロットしてみる。

# 周波数軸のデータ作成
freq = np.linspace(0, sr, sr)

# print(freq.shape)  # --> (16000,)

# plot
import matplotlib.pyplot as plt

plt.title('yes/0a7c2a8d_nohash_0.wav')
plt.xlabel('frequency(Hz)')
plt.ylabel('amplitude')
plt.plot(freq[:int(sr/2)+1], amp[:int(sr/2)+1]) # ナイキスト定数まで表示
plt.grid()
plt.tight_layout()
plt.show()

振幅が「1」を超えてしまっているので、正規化する。

正規化あり

# 振幅に変換
amp = np.abs(F)

# 正規化
N = len(data)  # データ数
amp_normal = amp / (N / 2)
amp_normal[0] /= 2 

短時間フーリエ変換(Short-Time Fourier Transform)

確認のため、音声データ全体を１フレームとみなして適用してみる。

# 短時間フーリエ変換
S_F = librosa.stft(data,
                  n_fft=16001,
                  win_length=16000,
                  hop_length=16000)

「n_fft」パラメータを「16000」にすると、2フレームになってしまうため、「16001」を指定する。

# 周波数binの確認
bins = librosa.fft_frequencies(sr=16000, n_fft=16001) 
print('bins: ', bins.shape)  # --> (8001,)

真ん中のDFT indexを境にして前半部分だけを使うため、周波数binの数としては半分になる。
work-in-progress.hatenablog.com

振幅に変換し、正規化してからプロットしてみる。

# 振幅に変換
amp = np.abs(S_F)
print(amp.shape)  # --> (8001, 1)

# 正規化
N = len(data)  # データ数
amp_normal = amp / (N / 2)
amp_normal[0] /= 2 

# plot
import matplotlib.pyplot as plt
import librosa.display

librosa.display.specshow(amp_normal,
                         sr=16000,
                         hop_length=16000, 
                         y_axis='linear')

plt.title('yes/0a7c2a8d_nohash_0.wav')
plt.colorbar(format='%+0.05f')
plt.ylim(0, 8000)
plt.tight_layout()
plt.show()

1000Hzまでの周波数帯を表示。

先ほどのFFTの結果と並べてみる。

相対値（デシベル）で表示してみる。
基準となる振幅は「0」とする。

# 振幅からデシベルへ変換
amp_normal_db = librosa.amplitude_to_db(amp_normal, ref=0)

# plot
import matplotlib.pyplot as plt
import librosa.display

librosa.display.specshow(amp_normal_db,
                         sr=16000,
                         hop_length=16000, 
                         y_axis='linear')

plt.title('yes/0a7c2a8d_nohash_0.wav')
plt.colorbar(format='%+2.0f').set_label('[dB]')
plt.ylim(0, 8000)
plt.tight_layout()
plt.show()

「1000」Hzまでを表示。
左が相対表示（デシベル）、右が絶対表示（振幅）

相対表示の基準を「振幅の最大値」に変えてみる。
（「librosa.amplitude_to_db」関数の第二引数「ref」を変更）

# 振幅からデシベルへ変換
amp_normal_db = librosa.amplitude_to_db(amp_normal, ref=np.max)

左が「最大値を基準」、右が「0を基準」

正規化をやめてみる。
（以下のコードをコメントアウトする）

# 正規化
N = len(data)  # データ数
amp_normal = amp / (N / 2)
amp_normal[0] /= 2 

左が「正規化なし」、右が「正規化あり」

基本的な動作が確認できたので、実践的な使い方をしてみる。

• DFTのブロック数： 「512」個（DFTのindexの数）
• フレームサイズ： 「480」サンプル（切り出す量、サンプリング周波数が「16000」の場合、時間換算で「30ms」）
• ホップサイズ： 「160」サンプル（ズラす量、サンプリング周波数が「16000」の場合、時間換算で「10ms」）

# 音声ファイルをロードする部分はこれまでと同様

# 周波数binの確認
bins = librosa.fft_frequencies(sr=16000, n_fft=512) 
print('bins: ', bins.shape)  # --> (257,)

# 短時間フーリエ変換
S_F = librosa.stft(data,
                  n_fft=512,
                  win_length=480,
                  hop_length=160)

amp = np.abs(S_F)
print(amp.shape)  # --> (257, 101)

# 振幅をデシベルに変換（基準値は「振幅の最大値」）
amp_db = librosa.amplitude_to_db(amp, ref=np.max)

# plot
import matplotlib.pyplot as plt
import librosa.display

librosa.display.specshow(amp_db,
                         sr=16000,
                         hop_length=160, 
                         y_axis='linear',
                         x_axis='time')

plt.title('yes/0a7c2a8d_nohash_0.wav')
plt.colorbar(format='%+2.0f').set_label('[dB]')
plt.ylim(0, 8000)
plt.tight_layout()
plt.show()

これが、「スペクトログラム」（Spectrogram）にあたる。
３次元のデータで、横軸に時間、縦軸に周波数、信号成分の"強さ"を色で表現する。

STFTをインプットとして使う実装例としては以下をご参照。
work-in-progress.hatenablog.com

matplotlibライブラリにある「specgram」関数でも「スペクトログラム」を表示できる。
「noverlap」は、隣合うフレームが重なる部分のサンプル数。
「DFTのブロック数」（「NFFT」パラメータ）と「フレームサイズ（＝フレームのサンプル数）」は同じになり、別々に指定できない模様。

import soundfile as sf

# Audio Data
audio_path = 'speech_dataset/yes/0a7c2a8d_nohash_0.wav'
data, sr = sf.read(audio_path)

import matplotlib.pyplot as plt
fig = plt.figure()

n_fft = 512      # frame sizeはこの値と同じになる 
hop_len = 160

spectrum, freqs, t, im = plt.specgram(data,
                                      NFFT=n_fft,
                                      Fs=sr,
                                      noverlap=n_fft-hop_len,
                                      mode='magnitude')

print('spectrum shape: ', spectrum.shape)  # --> (257, 97)
print('spectrum type: ', type(spectrum))   # -->  <class 'numpy.ndarray'>

# plot
plt.title('yes/0a7c2a8d_nohash_0.wav')
plt.xlabel('Time[sec]')
plt.ylabel('Frequency[Hz]')
fig.colorbar(im).set_label('[dB]')
plt.tight_layout()
plt.show()

プロット結果。
右は先ほどの「LibROSA」パッケージを使ったプロットの再掲。

以下はmodeを「PSD」（デフォルト）に変えたもの。

「"Power Spectral(Spectrum) Density"」については以下をご参照。
work-in-progress.hatenablog.com

コメントをいただいたので追記。

最大出力（振幅または音圧）となる値とその周波数を求めるには？

短時間フーリエ変換を行って、振幅を求めた2次元配列ampに対して、「最大値となるindexを探す」に読み替えられる。

amp = np.abs(S_F)
# --> shape: (DFT_index, frame_index)

最大値

np.max(amp)

index

(np.argmax(np.max(amp, axis=1)), np.argmax(np.max(amp, axis=0)))

「DFT_index」の刻み幅

(sampling_rate / DFT_size ) = (16000 / 512) =  31.25 [Hz]

「frame_index」の刻み幅

hop_length / sampling_rate = 160 / 16000  = 0.01 [sec]

音声「yes」でindexが(15, 34)と求まったとすると、

31.25  * 15 = 468.75 [Hz]  
0.01 * 34 = 0.34 [sec]

において、最大出力となる。