pythonで自作短時間フーリエ変換を用いたスペクトログラム関数作成！

こちらの記事で作成をした短時間フーリエ変換の関数を用いて、スペクトログラム関数を自作していきます！

pythonのnumpyで自作のSTFT関数とISTFT関数を実装！今年の4月に某大学のデジタル信号を扱う研究室に配属されました！希望が通り一安心です。そんな思いで研究をしています...

Contents

自作STFT関数
自作スペクトログラム関数
使用例
- 実装
- 結果
最後に

自作STFT関数

自作のSTFTです。

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

これを用いて、関数を作成します。

自作スペクトログラム関数

色ではデシベル表記になるように数値を計算し直しています。

import numpy as np
import matplotlib.pyplot as plt
from scipy import hamming

def my_spectrogram(x, window=1024, step=512, sr=48000, title="Spectrogram"):
    """
    plot spectrogram
    Parameters
    ----------
    x : array_like
        First one-dimensional input array.
    window : int, default 1024
        Length of window size
    step : int, default 512
        Length of step size
    sr : int, default 48000
        sampling rate of x
    title : str, default "Spectrogram
        title of the picture
    Returns
    -----------
    """
    win_fc = hamming(window)
    spectrogram = STFT(x, window=window, step=step, win_fc=win_fc)
    spectrogram = 20 * np.log10(np.abs(spectrogram.T))[window // 2:]

    # display range
    # left : 0[s], right : max_time[s]
    # bottom : 0[Hz], top : max_freq [Hz]
    max_time = x.shape[0] // sr
    max_freq = sr // 2 
    extent = [0, max_time, 0, max_freq]

    # display spectrogram
    plt.figure(figsize=[12, 6])
    ax = plt.subplot()
    im = plt.imshow(spectrogram, aspect="auto", extent=extent)
    ax.set_title(title, size=16)
    ax.set_xlabel('Time [s]', size=14)
    ax.set_ylabel('Frequency [Hz]', size=14)
    cbar = plt.colorbar(im, format='%+2.0f dB', aspect=100)
    cbar.set_label("Magnitude [dB]", size=14)
    plt.show()

こだわりポイントとしては、colorbar(カラーバー)を良い感じの位置に配置することにこだわりました。笑

使用例

実装

def main():
    fname = "signal.wav"
    x, sr = librosa.load(fname, sr=48000)  # [sec], [Hz]
    my_spectrogram(x, window=1024, step=512, sr=48000, title="Original signal spectrogram")
    plt.savefig("Original_signal_spectrogram.png")
    plt.close()

if __name__ == "__main__":
    main()