FFTは時系列データの周波数分析に重宝しますが、時間の変化はわかりません。スペクトログラムは周波数と時間変化の両方を分析可能な手法です。しかし計算方法は少々手順を踏みます。ここではPythonによるコーディングを1つずつ解説して行きます。
こんにちは。wat(@watlablog)です。
スペクトログラムは特に音声データ分析で絶大な効果を発揮する手法です。
ここではPythonによるスペクトログラムコードを自作したので紹介します！
ここでは周波数、時間、振幅の変化がひと目でわかるスペクトログラムを学びます。
しかし、Pythonの基本的なプログラミング手法を細かく説明しているわけではありません。もしPythonの基本文法等の学習をしたい場合はPython特化型学習サービスの「PyQ（パイキュー）」が大変オススメです！
目次（項目クリックでジャンプできます）
スペクトログラムは時系列データから多彩な情報を得ることができる！
スペクトログラムとは？
スペクトログラムとは、音声をはじめとする様々な時系列データの分析に使われる可視化手法のことです。
このスペクトログラムは時間×周波数×振幅レベルという3つの要素で構成された3次元グラフです。
このページの後半ではPythonでスペクトログラムを計算するコードを載せますが、以下にそのプログラムで得られるグラフの例を先に紹介します。
スペクトログラム表示する前の入力データはこちら。
この波形はグランドピアノの「A（ラ）」の音を最低音から最高音までオクターブ毎に順番に打鍵していったものです。
wavファイルの音声をお聞きしたい方は、以下のメディアプレーヤーから再生してみて下さい。
※生のグランドピアノの音が再生されます。
そして以下が上の時系列データをスペクトログラム表示した結果です。
このように、スペクトログラム表示をすることで時系列データだけでは時間と振幅の関係しかわからなかったものに、周波数軸という奥行きを持たせることができます。
キレイだ！
この時系列データの処理に今回はSTFT（短時間フーリエ変換）という手法を用いました。
STFT(短時間フーリエ変換)とは？
STFT（Short-Time Fourier Transform : 短時間フーリエ変換）とは、時系列データを細かく窓（フレーム）で区切ってFFTし、その窓を順次横にずらしながらFFTを繰り返していく解析方法です。
以下は時系列データからフレームでデータを切り出し、順次FFTしていくSTFTからスペクトログラムまでの過程を示す図です。
このように時系列データを抜き出す際にフレームをオーバーラップさせて計算するのが一般的です。
フレームで抜き出してきた短時間のフーリエ変換の結果は、時系列順に並べることで3次元の配列データになりますが、FFTの振幅レベルを色に変換することで横軸時間、縦軸周波数の2次元表示が可能になります。
スペクトログラムから得られる情報
スペクトログラムは世間で広く時系列データの分析に使われていますが、一般的にどんな情報を読み取っているのかを紹介します。
通常のFFTが周波数と振幅レベル（＋位相）という情報を持っているのに対し、スペクトログラムは時間の情報がプラスされています。
そのためスペクトログラムで主に分析する観点は「時間変化」になります。
以下は振幅レベルと周波数の時間変化の例です。
左の図は各ピーク（振幅レベルが卓越している部分）の振幅レベルが小さくなっていく（色が青くなっていく）時間がそれぞれ異なることを示しています。
右の図は時間が経つにつれピークが周波数方向に変化していくことを示しています。
このようにスペクトログラムではFFTで得られる情報に加えて、時間変化の情報を得ることができます。
Pythonでスペクトログラムを計算する方法
それではいよいよPythonによるSTFT計算とスペクトログラム表示までのプログラムを紹介します。
…とは言っても、実はWATLABブログでこれまで習得してきた内容の集大成となっていますので、それぞれ関連の記事を紹介していきます！
例により、ここに記載のプログラムは僕個人が学習して作った自作品です。間違いがある可能性もありますので使用には注意して下さい。
importするパッケージ
今回のSTFTとスペクトログラム表示のために、Pythonの外部ライブラリパッケージを宣言します。
インストールをしていない場合はpipで簡単にインストール可能なので、「Pythonのパッケージ管理ツール pipの使い方とコマンド集」の記事を確認して下さい。
|
1 2 3 4 |
import numpy as np from scipy import signal from scipy import fftpack import soundfile as sf |
配列処理のためのNumPyと科学技術計算のためのSciPy、wavファイルを使うためのPySoundFileを用意します。
時間波形データの用意(wavファイル読み込み)
まずは分析する時系列データを準備します。
具体的な方法は「PythonのPySoundFileでwavファイル波形表示」の記事で紹介した通りですが、コードを以下に示します。
|
1 2 3 |
def wavload(path): data, samplerate = sf.read(path) return data, samplerate |
オーバーラップ切り出し処理
時間波形を切り出す時にはオーバーラップ処理をできるようにしておきます。
オーバーラップ切り出し処理については「PythonでFFTをする前にオーバーラップ処理をしよう！」の記事に詳細を書いていますので参照下さい。
但し、上記記事とは戻り値に「最終時刻：final_time」を設定している部分に違いがあります。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
def ov(data, samplerate, Fs, overlap): Ts = len(data) / samplerate # 全データ長 Fc = Fs / samplerate # フレーム周期 x_ol = Fs * (1 - (overlap / 100)) # オーバーラップ時のフレームずらし幅 N_ave = int((Ts - (Fc * (overlap / 100))) /\ (Fc * (1 - (overlap / 100)))) # 抽出するフレーム数（平均化に使うデータ個数） array = [] # 抽出したデータを入れる空配列の定義 # forループでデータを抽出 for i in range(N_ave): ps = int(x_ol * i) # 切り出し位置をループ毎に更新 array.append(data[ps:ps + Fs:1]) # 切り出し位置psからフレームサイズ分抽出して配列に追加 final_time = (ps + Fs)/samplerate #切り出したデータの最終時刻 return array, N_ave, final_time # オーバーラップ抽出されたデータ配列とデータ個数、最終時間を戻り値にする |
最終時刻とは、以下の図に示すように、フレームをオーバーラップして切り出してきた時の最終フレームの終端における時刻です。
この時刻はスペクトログラムの時間軸を作成する時に使います。
フレームサイズを決めてデータを切り出していくと、データ長とフレームサイズ、オーバーラップ率によってはフレーム切り出しできない余った部分がでます。
スペクトログラム表示の正確な時間軸を作るためには、この余った部分を除外する必要があるのでこのような計算をしています。
STFT計算
上記までで時間波形切り出しの関数は準備できたので、続いて周波数分析としてFFTを行うための関数を用意します。
窓関数と振幅補正係数
FFTをかけるために窓関数としてハニング窓の関数を用意します。詳細は「PythonでFFT！SciPyで窓関数をかける」をご確認下さい。
そして窓関数を使って振幅レベルを適正に評価するためには、振幅補正係数acf（Amplitude Correction Factor）が必要です。acfの詳しい内容は「窓関数使用時の補正！FFTの時に忘れがちな計算とは？」に書きました。
|
1 2 3 4 5 6 7 8 9 |
def hanning(data_array, Fs, N_ave): han = signal.hann(Fs) # ハニング窓作成 acf = 1 / (sum(han) / Fs) # 振幅補正係数(Amplitude Correction Factor) # オーバーラップされた複数時間波形全てに窓関数をかける for i in range(N_ave): data_array[i] = data_array[i] * han # 窓関数をかける return data_array, acf |
デシベル(dB)変換
スペクトログラム表示をする場合は色成分をログスケールやdBスケールで表示させないと、音のような広い周波数帯域のデータを包括的に分析するには不向きになります（ピークの振幅レベルが10倍、100倍といった度合いで変化するのが普通なので）。
dBについての詳細は「Pythonで音圧のデシベル(dB)変換と逆変換式！」に記載しました。
|
1 2 3 |
def db(x, dBref): y = 20 * np.log10(x / dBref) # 変換式 return y # dB値を返す |
聴感補正（A特性補正）
聴感補正は音のデータを分析する時に通常使用します。人間の耳の特性に近づけた分析をすることで、実際の活用に役立てます。
聴感補正に関しては、「Pythonで聴感補正(A特性)の曲線を作る！」で式と実装を解説しました。
|
1 2 3 4 5 6 7 8 9 10 11 12 |
def aweightings(f): if f[0] == 0: f[0] = 1 else: pass ra = (np.power(12194, 2) * np.power(f, 4)) / \ ((np.power(f, 2) + np.power(20.6, 2)) * \ np.sqrt((np.power(f, 2) + np.power(107.7, 2)) * \ (np.power(f, 2) + np.power(737.9, 2))) * \ (np.power(f, 2) + np.power(12194, 2))) a = 20 * np.log10(ra) + 2.00 return a |
FFT処理
FFTの関数は「PythonでFFT！SciPyのFFTまとめ」に記載した内容とほぼ同じですが、今回は聴感補正を考慮するので、最初からdB形式で計算するため、以下の順番で演算します。
①振幅補正係数acfをかける。
②dB変換する（聴感補正がdB値なので、dB同士のデータにする）。
③聴感補正曲線を足す。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
def fft_ave(data_array, samplerate, Fs, N_ave, acf): fft_array = [] fft_axis = np.linspace(0, samplerate, Fs) # 周波数軸を作成 a_scale = aweightings(fft_axis) # 聴感補正曲線を計算 # FFTをして配列にdBで追加、窓関数補正値をかけ、(Fs/2)の正規化を実施。 for i in range(N_ave): fft_array.append(db\ (acf * np.abs(fftpack.fft(data_array[i]) / (Fs / 2))\ , 2e-5)) fft_array = np.array(fft_array) + a_scale # 型をndarrayに変換しA特性をかける fft_mean = np.mean(fft_array, axis=0) # 全てのFFT波形の平均を計算 return fft_array, fft_mean, fft_axis |
ndarray形式であればforループを使わなくても配列演算ができるので、効率的ですね。
メイン実行文とグラフ表示
以上までのdef文で作成したコードを「function.py」というファイルで用意し、以下に示すメインコード（main_spectrogram.py）を実行することで冒頭に紹介したスペクトログラムを得ることができます（import文で関数をまとめたファイルを宣言する必要はあります）。
関数の使い方は「Pythonの関数 def文の使い方！引数や別ファイル式も解説」で説明しています。
wavファイルが無くても、自分で波形を作成（samplerate, x, dataを決めるだけ）して使っても動作しますので、是非お試し下さい。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 |
import function import numpy as np from matplotlib import pyplot as plt path = 'ff_A.wav' #ファイルパスを指定 data, samplerate = function.wavload(path) #wavファイルを読み込む x = np.arange(0, len(data)) / samplerate #波形生成のための時間軸の作成 # Fsとoverlapでスペクトログラムの分解能を調整する。 Fs = 4096 # フレームサイズ overlap = 90 # オーバーラップ率 # オーバーラップ抽出された時間波形配列 time_array, N_ave, final_time = function.ov(data, samplerate, Fs, overlap) # ハニング窓関数をかける time_array, acf = function.hanning(time_array, Fs, N_ave) # FFTをかける fft_array, fft_mean, fft_axis = function.fft_ave(time_array, samplerate, Fs, N_ave, acf) # スペクトログラムで縦軸周波数、横軸時間にするためにデータを転置 fft_array = fft_array.T # ここからグラフ描画 # グラフをオブジェクト指向で作成する。 fig = plt.figure() ax1 = fig.add_subplot(111) # データをプロットする。 im = ax1.imshow(fft_array, \ vmin = 0, vmax = 60, extent = [0, final_time, 0, samplerate], \ aspect = 'auto',\ cmap = 'jet') # カラーバーを設定する。 cbar = fig.colorbar(im) cbar.set_label('SPL [dBA]') # 軸設定する。 ax1.set_xlabel('Time [s]') ax1.set_ylabel('Frequency [Hz]') # スケールの設定をする。 ax1.set_xticks(np.arange(0, 50, 5)) ax1.set_yticks(np.arange(0, 20000, 2000)) ax1.set_xlim(0, 30) ax1.set_ylim(0, 10000) # グラフを表示する。 plt.show() plt.close() |
追記：wavファイルが無いけどすぐ試したい方は自分で時間波形を作ろう！
メインのコードを変更
スペクトログラムはwavファイルでなくても描くことができます。上記のコードはモノラルのwavファイルを使うことが前提の内容でしたが、wavファイルって一回録音して、ファイル移動して…ってひと手間が面倒ですよね。
手っ取り早くお手元の環境でSTFT解析を行い、スペクトログラム表示を試したい方のために、以下にwavファイルが無くても動作するメインの.pyファイルを示します。
コードの冒頭でscipyをimportし、chirp信号という時間で周波数が変化する信号を生成しています。
チャープ信号については「Pythonでチャープ信号！周波数スイープ正弦波の作り方」に詳しく書いてありますので、是非併せてご覧下さい！
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 |
import function import numpy as np from matplotlib import pyplot as plt from scipy.signal import chirp # チャープ信号を作成 t = np.linspace(0, 5, 128000) data = chirp(t, f0=1, f1=2000, t1=5, method='linear') samplerate = 25600 x = np.arange(0, len(data)) / samplerate #波形生成のための時間軸の作成 # Fsとoverlapでスペクトログラムの分解能を調整する。 Fs = 4096 # フレームサイズ overlap = 75 # オーバーラップ率 # オーバーラップ抽出された時間波形配列 time_array, N_ave, final_time = function.ov(data, samplerate, Fs, overlap) # ハニング窓関数をかける time_array, acf = function.hanning(time_array, Fs, N_ave) # FFTをかける fft_array, fft_mean, fft_axis = function.fft_ave(time_array, samplerate, Fs, N_ave, acf) # スペクトログラムで縦軸周波数、横軸時間にするためにデータを転置 fft_array = fft_array.T # ここからグラフ描画 # グラフをオブジェクト指向で作成する。 fig = plt.figure() ax1 = fig.add_subplot(111) # データをプロットする。 im = ax1.imshow(fft_array, \ vmin = -10, vmax = 60, extent = [0, final_time, 0, samplerate], \ aspect = 'auto',\ cmap = 'jet') # カラーバーを設定する。 cbar = fig.colorbar(im) cbar.set_label('SPL [dBA]') # 軸設定する。 ax1.set_xlabel('Time [s]') ax1.set_ylabel('Frequency [Hz]') # スケールの設定をする。 ax1.set_xticks(np.arange(0, 120, 1)) ax1.set_yticks(np.arange(0, 20000, 500)) ax1.set_xlim(0, 5) ax1.set_ylim(0, 2000) # グラフを表示する。 plt.show() plt.close() |
実行すると、以下のグラフが表示されると思います。5[s]間で0[Hz]から2000[Hz]まで周波数が増加している様子がはっきりと確認できますね！
補足：カラーマップの種類について
スペクトログラムは色を使ってデータを表現します。
この色のパターンは非常に多くの種類があり、imshow関数内で引数としているcmapの指定で変更することが可能です。
以下の図はcmap違いのスペクトログラムです。「参考外部サイト様：color exsample code」に良いリストがありますので、お好みの設定を見つけて下さい。
まとめ
少し長くなってしまいましたが、本記事では音声データのスペクトログラムを計算する一連の方法を書きました。
しかし新しいことはほとんど無く、これまでWATLABブログで扱ってきた内容（wavファイルの取り扱い、波形切り出し、窓関数の取り扱い、dB変換、聴感補正、多段のFFT）が多くありました。
カラーマップを変更することでスペクトログラムの見え方を大きく変えることができました。
今回は音声データの強力な分析手法であるSTFTやスペクトログラム表示を習得しました！これで色々なデータから様々な情報を得ることができるようになるね！
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！
質問させてください！
「PythonでFFTをする前にオーバーラップ処理をしよう！」でも記述してあったのですが、
FFTをした後に平均化処理をする理由がよく理解できません。お教えしていただけたら幸い
です。よろしくお願いします。
watです。記事を読んで頂きありがとうございます。
平均化ですが、ここに書いてある内容は主に実験データを想定しています。
理想波形であれば1フレームのFFTも複数フレームを平均化したFFTも変わりません。
しかし、実験データは時間によって主要次数の振幅が安定していないことも多々あります。
平均化処理は、不安定なデータを綺麗にするために行います。
おそらく大抵の商用FFTアナライザー（100万円くらいするやつ）にはこのような機能が付いていると思います。
回答が的を射てなかったらすみません…。その場合はまた質問して下さい！
なるほど！！
理解することができました！ありがとうございます！
ご連絡ありがとうございます！
また何かわからないことがありましたらコメントよろしくお願いします！
突然失礼します。ここにあるコードでfftを行ったのですが、 data_array[i]= data_array[i] * han # 窓関数をかける
ValueError: operands could not be broadcast together with shapes (4096,2) (4096,) というエラーが出て困っています。このエラーの解決方法について助言をいただけたら幸いです。
コメントありがとうございます。
過去に似たような質問を頂いたことがあるのですが、ステレオ録音して作成したwavファイルをご使用していますでしょうか？
ここに記載のコードはモノラルのwavにのみ対応していますので、まずはそこをご確認下さい。
dataを定義した直後に、print(data.ndim)とすることでdataの次元を調査することも可能です。その次元が2とかですとステレオのwavである可能性が高いです。
また、もしwavファイルが原因で発生しているエラーでしたら、記事内の「追記：wavファイルが無いけどすぐ試したい方は自分で時間波形を作ろう！」の項目に記載してあるコードをお試し下さい。
このコードはwavファイルを使わないでSTFT計算をするバージョンです。
もしこのコードがエラーなく実行できれば、他は問題ないということになります。
もし僕のこの回答が検討はずれである場合は再度コメント頂ければと思います！よろしくお願いします！
本当にありがとうございます！
ご指摘の通り、ステレオ録音で作成したwavファイルだったので、
data, samplerate = function.wavload(path) #wavファイルを読み込む
と
x = np.arange(0, len(data)) / samplerate #波形生成のための時間軸の作成
の間に
♯ステレオ2chの場合、LchとRchに分割
data_l = data[:, 0]
data_r = data[:, 1]
♯ 入力をモノラル化
data = (0.5 * data_l) + (0.5 * data_r)
を入れて、ファイルをモノラルに変えたら無事動きました。
ご連絡ありがとうございます。
解決したようでよかったです！また何かご不明点ありましたらコメント下さい！