python

pythonのnumpyで自作のSTFT関数とISTFT関数を実装!

今年の4月に某大学のデジタル信号を扱う研究室に配属されました!

希望が通り一安心です。

センサデータに強いデータサイエンティストになりたい!

そんな思いで研究をしています。

このブログでは学習した内容をまとめていきたいと思います!

 

今回は、信号データの基本である短時間フーリエ変換(STFT), 短時間逆フーリエ変換(ISTFT)を実装します。

今回は、np.fft.fftは使用します。

そう遠くない未来に、C言語でft, fftは実装する予定です。

信号をnumpy形式で読み込む

windows10 homeの標準ソフトボイスレコーダーを用いて音声を自分のノートPCで録音

  •   サンプリング周波数 48000 [Hz]
  •   16bit
  •   2チャネル
from scipy import signal
import librosa
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

name = "signal.wav"
wav, sr = librosa.load(fname, sr=48000)
こんにちは。初めまして

と言っています。

窓関数、窓幅、シフト幅を決定

window = 1024
step = 512
i = 0
win_fc = signal.hamming(window)
tmp = wav[i*step:i*step+window]
tmp = tmp*win_fc

 

どんなことをしているのかを画像で確認をすると、$i = 0$のときでは

左から、窓関数、オリジナルの波形、窓関数を適応した波形です。

これをシフト幅ごとに連続で繰り返してあげると短時間フーリエ変換が完了します!

それでは、関数化します。

def STFT(x, window=1024, step=512, win_fc=1.):
    """
    Compute the Short Time Fourier Transform (STFT).
    Parameter
    -----------
    x : array-like
        Time series of measurement values
    window : int, default 1024
        Length of window size
    step : int, default 512
        Length of step size
    win_fc : array-like
        window function
    Return
    ------------
    Z : ndarray
        STFT of x. By default, the last axis of Zxx corresponds to the segment times.
    """
    Z = []
    for i in range((x.shape[0] - window) // step):
        tmp = x[i*step : i*step + window]
        tmp = tmp * win_fc
        tmp = np.fft.fft(tmp)  # get fft data [i*step:i*step+window]
        Z.append(tmp)
    Z = np.array(Z)
    return Z

 

同様に短時間逆フーリエ変換も定義します。

def ISTFT(Z, window=1024, step=512):
    """
    Compute the Invert Short Time Fourier Transform (ISTFT).
    Parameter
    -----------
    Z : ndarray
        STFT of x. By default, the last axis of Zxx corresponds to the segment times.
    window : int, default 1024
        Length of window size
    step : int, default 512
        Length of step size
    Return
    ------------
    x : array-like
        Time series of measurement values
    """
    ifftList = [np.fft.ifft(t) for t in Z]
    x = np.zeros(Z.shape[0] * step + window)
    for i, ifft in enumerate(ifftList):
        x[i*step : i*step + window] += ifft.real
    return x

 

使用例

実装

def main():
    # load file
    fname = "signal.wav"
    wav, sr = librosa.load(fname, sr=48000)
    # calculate STFT
    window = 1024
    step = 512
    win_fc = signal.hamming(window)
    ampList = STFT(wav, window=window, step=step, win_fc=win_fc)
    # freq range
    freq = np.fft.fftfreq(wav.shape[0], 1.0 / sr)
    freq = freq[:int(freq.shape[0] / 2)]
    cols = [freq[j] for j in range(0, len(freq), len(freq) // step)]
    cols = [int(c) for c in cols]
    # time course
    time = np.arange(0, ampList.shape[0], 1) * step / sr
    time = [f'{t:.2f}' for t in time]
    # plot spectrogram
    df_amp = pd.DataFrame(
        data=abs(ampList[:, :ampList.shape // 2]), index=time, columns=cols)
    plt.figure(figsize=(12, 4))
    plt.title('Spectrogram')
    grid = sns.heatmap(data=20 * np.log10(df_amp.T), xticklabels=50,
                       yticklabels=100, cmap=plt.cm.gist_rainbow_r)
    grid.set(xlabel='Time [sec]', ylabel='Frequency [Hz]')
    plt.savefig('spectrogram.png')
    # calculate STFT
    wav_hat = ISTFT(ampList, window=window, step=step)
    # compair original signal and re-synthesized signal
    plt.figure(figsize=(15, 10))
    plt.subplot(2, 1, 1)
    plt.plot(wav)
    plt.title('Original signal')
    plt.xlabel('Time [sec]')
    plt.ylabel('Magnitude')

    plt.subplot(2, 1, 2)
    plt.plot(wav)
    plt.title('Re-synthesized signal')
    plt.xlabel('Time [sec]')
    plt.ylabel('Magnitude')
    plt.savefig('compair_signal.png')
    # check the correlation
    print("corr:{}".format(np.corrcoef(x=wav[:len(wav_hat)], y=wav_hat)))

if __name__ == "__main__":
    main()

 

上記では、

  1. 波形からスペクトログラムを作成してプロット
  2. 短時間フーリエ変換した結果を逆変換して元の波形との相関係数を調べて正しく動作しているかどうかを調査

しています。

短時間フーリエ変換の結果は、窓幅に対して対称な出力となるので、スペクトログラムを作成する際は、半分だけを使用する点に注意が必要です。

それぞれの結果を見てみます。

波形からスペクトログラムを作成してプロット

seabornは出力がきれいなので好きです。ただ、画像を出力するまでにかなり時間がかかるので、matplotlib.pyplot.imshowを使った方が良いかもしれません。

もし、スペクトログラムを関数化するならそうします。

短時間フーリエ変換した結果を逆変換して元の波形と比較

相関係数は0.9999996896225544となりました。

ほぼ完璧に復元できていることが分かります。

最後に

短時間フーリエ変換STFT、短時間逆フーリエ変換ISTFTを実装しました。

意外と実装簡単だなというのが感想です。

 

最近だと、Parallel WaveGANの損失関数にも使われていました。

最新の知識を取り入れるためにも基礎的な知識の積み重ねが大切ですね。

pythonで自作短時間フーリエ変換を用いたスペクトログラム関数作成!こちらの記事で作成をした短時間フーリエ変換の関数を用いて、スペクトログラム関数を自作していきます! https://ibuki-s...

こちらの記事でスペクトログラム関数を作成しました。

もし、興味があれば見てみてください!

オススメのプログラミングスクールをご紹介

タイピングもままならない完全にプログラミング初心者から

アホいぶきんぐ
アホいぶきんぐ
プログラミングってどこの国の言語なの~?

たった二ヶ月で

いぶきんぐ
いぶきんぐ
え!?人工知能めっちゃ簡単にできるじゃん!

応用も簡単にできる…!!

という状態になるまで、一気に成長させてくれたオススメのプログラミングスクールをご紹介します!

テックアカデミーのPython+AIコースを受講した僕が本音のレビュー・割引あり! というプログラミング完全初心者だった僕が Tech Academy(テックアカデミー)のPython×AIコース を二ヶ月間...

COMMENT

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です