今回はSTFT(Short-Time Fourier Transform: 短時間フーリエ変換)を使って読み込んだ音楽ファイルのスペクトログラムを見てみましょう。
音楽ファイルの読み込みについては前回の記事を参照してください。
では早速本題です。
目次
STFT(短時間フーリエ変換)とは
時間領域の信号に対してフーリエ変換をすることで，その信号を純音(sinやcos)の足し合わせで表現できます。
ここで一度フーリエ変換の式を復習しておきましょう。
$$F(\omega) = \int_{-\infty}^{\infty} f(t) e^{-i\omega t} dt$$
これは時間領域の信号を周波数分析出来ることを表しています。
更に周波数領域での分析は人間の聴覚の仕組み(基底膜の振動の原理)と非常に親和性が良く，そのためイコライザーなんかは良く使われます。
しかしこのフーリエ変換は対象の信号全体を周波数分析することに相当しており，局所時間での周波数情報は落ちてしまっています。
これは例えば1～5秒の区間では440Hzの音が鳴っていて5～10秒の区間では800Hzの音が鳴っているという状況をフーリエ変換では見ることができないということになります。
そこでSTFTの登場です！
STFTは短時間についてフーリエ変換を行うことで，局所時間の周波数分析を行うことが出来ます。
これだけ聞くとメリットしかないように感じるかもしれませんが，実はデメリットもあります。
それは不確定性原理と言われるもので，簡単に言うと短時間で周波数解析するために時間分解能は良くなるが周波数分解能が悪くなるという問題です。
これは時間周波数解析手法で必ず問題になるもので，ウィグナー分布やウェーブレット変換など様々な手法が提案されています。
ではSTFTの式を見ていきましょう。
今回は計算機で実装することを考えているので離散時間のSTFTのみご紹介します。
※離散時間という事を意識するために変数を[]で表記しています
$$F[m,k] = \sum_{n=0}^{N-1} f_{m}[n] e^{-i2\pi kn/N} $$
ここで$m$は時間フレームのインデックス，$k$は離散周波数を表すインデックス，$N$はフレーム長を表します。
更にフレームごとの関数$f_{m}[n]$は
$$f_{m}[n] := w[n] f[n + mS]$$
と窓関数$w[n]$(${\rm supp} (w[n]) = [0,N-1]$)と信号の掛け算の形で定義しています。
ここで$S$はフレームシフト量です。
離散フーリエ変換と見比べてみると変数が時間インデックス$m$と周波数インデックス$k$の2つになっているところが大きく変わってますね。
スペクトログラムでわかり易く表現
スペクトログラムとは信号を時間周波数分析をすることで得た情報(時間，周波数，信号の強さ)を図として描画したものです。
3次元情報を2次元描画するために，基本的には横軸に時間，縦軸に周波数，色によって信号の強さを表現します。
ここで信号の強さとは$F[m,k]$の絶対値の2乗振幅($|F[m,k]|^{2}$)を表すことが多いです。※これをパワースペクトル密度と呼んだりもします
では実際に信号のスペクトログラムを見てみましょう。
以下はmacに標準搭載の音声読み上げ機能に「もろみ先輩」と読み上げてもらった信号のスペクトログラムです。
ちなみに以下のコマンドを使いました。
$ say -o moromi_senpy.wav --data-format=LEI24@48000 "もろみ先輩"
これは「もろみ先輩」と読み上げた音声をカレントディレクトリに「moromi_senpy.wav」というファイル名で，データフォーマットを量子化ビット数24bitのサンプリング周波数48000Hzで保存して下さいという意味の命令です。
音声はこちらです。
音声研究をされてる人が見たらスペクトログラムから何を言っているか分かるのかもしれませんが，僕には全くわかりません笑
pythonによるSTFTの実装
それでは実際にpythonで実装してみましょう。
今回は細かいところは解説しませんが，もし何か質問があればコメントでお尋ねください。
# STFT (s: signal(1D-array), Lf: length of frame(window), noverlap: number of overlap)
def STFT(s, Lf, noverlap=None):
if noverlap==None:
noverlap = Lf//2
l = s.shape[0]
win = np.hanning(Lf)
Mf = Lf//2 + 1
Nf = int(np.ceil((l-noverlap)/(Lf-noverlap)))-1
S = np.empty([Mf, Nf], dtype=np.complex128)
for n in tqdm(range(Nf)):
S[:,n] = np.fft.rfft(s[(Lf-noverlap)*n:(Lf-noverlap)*n+Lf] * win, n=Lf, axis=0)
return S
これが信号sについてSTFTを行った結果を返す関数です。
Lfはフレーム長(上式のN)，noverlapはフレームシフト量(上式のS)です。
窓関数$w[n]$にはハニング窓を使用しております。
また，最後のフレームについては今回切り捨てる仕様にしておりますが，必要であれば簡単に変更できるかと思います。
※tqdmという関数は進捗状況を見るための関数ですので消してしまっても構いません
スペクトログラム描画部分
では次はこのSからスペクトログラムを描画するコードです。
# plot spectrogram (fs: sampling frequency, s: signal(1D-array), Lf: length of frame(window), noverlap: number of overlap)
def plot_spectrogram(fs, s, Lf, noverlap=None):
S = STFT(s, Lf, noverlap)
S = S + 1e-18
P = 20 * np.log10(np.abs(S))
P = P - np.max(P) # normalization
vmin = -150
if np.min(P) > vmin:
vmin = np.min(P)
m = np.linspace(0, s.shape[0]/fs, num=P.shape[1])
k = np.linspace(0, fs/2, num=P.shape[0])
plt.figure()
plt.pcolormesh(m, k, P, cmap = 'jet', vmin=-150, vmax=0)
plt.title("Spectrogram of Sound")
plt.xlabel("time[s]")
plt.ylabel("frequency[Hz]")
plt.colorbar()
plt.tight_layout()
plt.show()
ちょっとした工夫としては対数を取るためにSの要素が0にならないよう微小な数を足しているところ，最大が0[dB]になるよう正規化しているところくらいでしょうか。
また，-150[dB]以下は見えないようにしているので使う場面によってここは変更しても良いかもしれません。
関数の使用例
最後にこれら関数の使用例を見ていきましょう。
当然ファイルを読み込む必要がありますのでご注意ください。
以下サンプルコードです。
file_name = "./sample.wav" start = 0 # [s] end = 0 # [s] (if end==0 => full size) Lf = 2**8 # length of frame(window) noverlap = None s = wf.readWave(file_name, start, end) params = wf.getParams(file_name) # params = (nchannels, sampwidth, framerate, nframes, comptype, compname) plot_spectrogram(params[2], s, Lf, noverlap)
ここではフレーム長を2の8乗(256)に，フレームシフト量についてはこれの半分(128)にしています(STFT関数をご参照ください)。
ここの設定については時間分解能や周波数分解能に関わる部分になりますので状況によって変更が必要です。
蛇足になるかもしれませんが，最後に必要なパッケージのインポート部分をお示ししておきます。
import numpy as np from tqdm import tqdm import matplotlib.pyplot as plt import wave_func as wf
tqdmが入っていない場合にはpip等でインストールするか，若しくは使用しないというのもありです。
まとめ
今回はpythonによるSTFT及びスペクトログラム描画の実装をご説明しました。
少し高度な内容ですが，まずはプログラムなどの実践を通して慣れていくのも工学に重要かと思います。
この記事が誰かの趣味や研究の参考になれば幸いです。
ではお疲れ様でした。
最後に，STFTについて学ぶのであれば以下の本がわかりやすいです。
※時間周波数分析全般を詳しく書いている和書はそれ程ありません
