【python】音声データをFFT解析してスぺクトログラムを作成する方法

はじめに

以前は自分で作った周波数不変のサイン波でFFT解析を行いました。

今回は音声データを用いてFFT解析を行っていきたいと思います。
音声データでは遂次周波数成分が変わっていくので短時間で区切ってFFT解析を行います。
それを時系列でプロットしたものがスペクトログラムとよばれ、今回はこれを作っていきます。
やっていきましょう

はじめに
音声データの作成
音声データのFFT解析
- 音声データの読み込み
- 短時間フーリエ変換
おわりに

音声データの作成

今回はちゃんとスペクトログラムができているかの確認のため、音楽とかじゃなくて自分で作った音声信号データを用います。
フリーソフトのaudacityを使って下記の音声データを作成しました。
tone.wavというファイル名で書きだしておきます。

サンプルレート44100Hzのモノラル信号
440Hz, 340Hz, 540Hz, 640Hzのサイン波を各1秒間続けて鳴らす

audacityは以下リンクからダウンロードできます

音声データのFFT解析

まず、音声データの読み込みから行います。

音声データの読み込み

wavファイルをnumpyで扱えるような形式で読み込む方法はいくつかありますが、Pysoundfileというライブラリを使うのが簡単なようなのでそれを使います。

import soundfile as sf
import numpy as np
import matplotlib.pyplot as plt

filepath = "tone.wav"
data, fs = sf.read(filepath)

これだけでdataに音声の波形データがnumpy array形式で、fsにサンプリング周波数を取得できます(モノラル音声の場合)
ステレオの場合は、dataが2次元行列になってます。
試しに波形データをプロットしてみましょう。
時系列データはdataのデータ数をサンプリング周期から作成します。

t = np.arange(0, len(data)/fs, 1/fs)
fig, ax = plt.subplots(2, 2)
ax[0,0].plot(t, data)
ax[0,0].set_xlim([0, 0.1])
ax[0,0].set_title('440Hz')
ax[0,0].set_xlabel('time [s]')
ax[0,1].plot(t, data)
ax[0,1].set_xlim([1.1, 1.2])
ax[0,1].set_title('340Hz')
ax[0,1].set_xlabel('time [s]')
ax[1,0].plot(t, data)
ax[1,0].set_xlim([2.4, 2.5])
ax[1,0].set_title('540Hz')
ax[1,0].set_xlabel('time [s]')
ax[1,1].plot(t, data)
ax[1,1].set_xlim([3.5, 3.6])
ax[1,1].set_title('640Hz')
ax[1,1].set_xlabel('time [s]')
plt.tight_layout()
plt.show()

短時間フーリエ変換

取得した音声データを短時間ずつ区切ってFFT解析を行っていきます。
高速フーリエ変換を行いたいのでデータ数は2の累乗の値を使用します。
今回はN=1024とします。

音声データを時系列順に少しオーバーラップさせながらN個抽出し、
得られたパワースペクトルを行列に追加していきます。
オーバーラップ量は適当ですが、Nの1/4程度にしておきます。
コードは下記です。

from scipy import signal
# fft
N =1024            # サンプル数
window = signal.hann(N)  # ハニング窓関数

# 高速フーリエ変換
st = 0
shift = N // 4

spec = np.zeros([len(data) // shift - (N // shift - 1), N // 2])

for i in range(len(data) // shift - (N // shift - 1)):
    st = i * shift
    x = data[st:st + N] * window
    F = np.fft.fft(x)
    freq = np.fft.fftfreq(N, d=1 / fs) # 周波数スケール
    # フーリエ変換の結果を正規化
    F = F / (N / 2)
    F = F * (N / sum(window))
    Amp = np.abs(F)
    # パワースペクトルの計算（振幅スペクトルの二乗）
    Pow = Amp ** 2
    # デシベル変換
    Pow = 20 * np.log10(Pow / (20 * 10 ** (-6)))
    spec[i,:] = Pow[:N // 2]

y_max = 8000
fig, ax = plt.subplots()
im = ax.imshow(spec.T, origin='lower',
               extent=(0, (st + N) / fs, 0, freq[N // 2 - 1]),
               aspect = ((st + N) / fs) / y_max,
               vmin=0, vmax=100)
ax.set_ylim([0, y_max])
ax.set_title('Power Spectrogram')
ax.set_xlabel('time [s]')
ax.set_ylabel('Frequency [Hz]')
plt.tight_layout()
plt.colorbar(im)
plt.show()