いろんなものはつながっている

pythonでスペクトログラムを表示してみる

scipy.signal.spectrogramでスペクトログラムデータを作成してpcolormeshで表示する。

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


# 指定桁で切り上げする
def ciel_quantize(src, digit):
    tmp = src/digit;
    tmp = np.math.trunc(tmp)
    if tmp >= 0: 
        tmp = (int)((tmp + 1) * 10)
    else:
        tmp = (int)((tmp - 1) * 10)

    return tmp

    
# 信号にノイズを加える
def add_noise(sig, sn):
    # 信号の電力(dB)
    norm_sig=np.linalg.norm(sig)
    pw_sig_dB=20*np.math.log10(norm_sig)
    # ノイズ生成
    sig_length=len(sig)
    noise=np.random.randn(sig_length)
    pw_noise_dB=20*np.math.log10(np.linalg.norm(noise))
    # ノイズの電力を指定のSNになるように調整する
    pw_noise_dB_delta=(pw_sig_dB-sn)-pw_noise_dB
    noise=noise*10**(pw_noise_dB_delta/20)
    sig_with_noise=sig+noise

    return sig_with_noise


if __name__ == '__main__':

    #サンプル信号生成
    sa=16000
    fft_bin=2048
    sn=10
    sig=np.exp(2*np.pi*1j*(100/fft_bin)*np.arange(0,sa*10))
    norm_sig=np.linalg.norm(sig)
    #sig=sig/norm_sig
    sig = add_noise(sig, sn)
   
    #スペクトログラム生成
    f_list, t_list, sxx = signal.spectrogram(sig, sa, nperseg=fft_bin,\
                          noverlap=28, return_onesided=False)

    #データの並びは0~sa/2, -1*sa/2~0 になっているので
    #-1*sa/2~sa/2に並びかえる
    #また、転置をして(周波数,時間)にする
    sxx = sxx.T
    sxx = np.fft.fftshift(sxx)
    f_list = np.fft.fftshift(f_list)

    #対数化
    sxx_log10=20*np.log10(np.abs(sxx))

    # 表示のためにノイズレベル、最大値を求める
    # 最頻値をノイズレベルと仮定する
    sudo_noise_level_dB = np.median(sxx_log10.flatten()) 
    sudo_noise_level_dB =ciel_quantize(sudo_noise_level_dB, 10)
    # 電力の最大値
    max_pw_dB = np.max(sxx_log10.flatten())
    max_pw_dB = ciel_quantize(max_pw_dB, 10)

    #スペクトログラムを表示する
    plt.pcolormesh(f_list_, t_list, sxx_log10, \ 
vmin=sudo_noise_level_dB, vmax = max_pw_dB)
    plt.colorbar()
    plt.show()

関連記事

コメント

  1. この記事へのコメントはありません。

  1. この記事へのトラックバックはありません。

スポンサード リンク

カテゴリー

スポンサード リンク