今回も基本シリーズ。
昨夜同様に、STFTを正しく表示させることを主眼にまとめる。
以下で解説するコードを使うと、’ひらけごま’の音声wavfileに対して次の図が得られる。

コード解説

STFTのために第一行目のfrom scipy import signalを利用する以外は、基本的に利用するものは前回のFFTと同様です。

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
import wave
import matplotlib
from scipy.fftpack import fft, ifft

wavfileの読み込みもFFTの場合と同じです。

wavfile = 'hirakegoma'

def wave_input(wavfile):
    wr = wave.open(wavfile+'.wav', "rb")
    ch = wr.getnchannels()
    width = wr.getsampwidth()
    fr = wr.getframerate() #sampling freq ; RATE
    fn = wr.getnframes()  #sampling No. of frames; CHUNK
    fs = fn / fr  #sampling time
    origin = wr.readframes(fn)
    data = origin  #[:fn]
    print("fs",fs)
    print("fr",fr)
    # ステレオ前提 > monoral
    sig = np.frombuffer(data, dtype="int16")  /32768.0
    t1 = np.linspace(0,fs, fn) #, endpoint=False)
    return t1,sig,fn,fr 
t1,sig,fn,fr=wave_input(wavfile)

描画するために座標等を定義します。

amp=np.max(sig)
matplotlib.rcParams.update({'font.size': 18, 'font.family': 'sans', 'text.usetex': False}) 
fig = plt.figure(figsize=(12,12))  #(width,height)

axes1 = fig.add_axes([0.1, 0.6, 0.7, 0.4]) #sig vs t;normal size
axes2 = fig.add_axes([0.15, 0.75, 0.2, 0.2]) #sig vs t;insert 
axes3 = fig.add_axes([0.1, 0.1, 0.7, 0.4]) #stft vs t;main
axes4 = fig.add_axes([0.81, 0.1, 0.03, 0.4]) #colorbar(main)
axes5 = fig.add_axes([0.15, 0.25, 0.2, 0.2]) #stft vs t;insert
axes6 = fig.add_axes([0.91, 0.1, 0.03, 0.4]) #colorbar(insert)
axes7 = fig.add_axes([0.55, 0.25, 0.2, 0.2]) #FFT vs freq;insert

信号を描画します。描画方法は前回までと同様です。
※sig vs tの挿入図の範囲は[1,1.1]としており、以下のFFTの挿入図の範囲は[2.5,2.8]であり異なっていますが、このまま掲載します

axes1.plot(t1, sig)
axes1.grid(True)
axes2.set_xlim([1, 1.1])
axes2.set_ylim(-amp*1.1,amp*1.1)
axes1.set_xlim(0,5)
axes2.plot(t1,sig)

以下がFFT関数です。

def FFT(sig,fn,fr):
    freq =fft(sig,fn)
    Pyy = np.sqrt(freq*freq.conj())/fn
    f = np.arange(0,fr,fr/fn)

以下でピーク振動数と振幅を求めます。ここも前回と同様です。
order=10は（両サイドの）連続する10個のデータ中相対値が一番大きな値を求めます；10回で勝ち抜け

    ld = signal.argrelmax(Pyy, order=10) #相対強度の最大な番号をorder=10で求める
    ssk=0
    fsk=[]
    Psk=[]
    maxPyy=max(np.abs(Pyy))
    for i in range(len(ld[0])):  #ピークの中で以下の条件に合うピークの周波数fと強度Pyyを求める
        if np.abs(Pyy[ld[0][i]])>0.25*maxPyy and f[ld[0][i]]<20000: # and f[ld[0][i]]>20:
            fssk=f[ld[0][i]]
            Pssk=np.abs(Pyy[ld[0][i]])
            fsk.append(fssk)
            Psk.append(Pssk)
            ssk += 1
    print('{}'.format(np.round(fsk[:len(fsk)],decimals=3))) #標準出力にピーク周波数fskを小数点以下二桁まで出力する
    print('{}'.format(np.round(Psk[:len(fsk)],decimals=4))) #標準出力にピーク強度Pskを小数点以下6桁まで出力する
    return  freq,Pyy,fsk,Psk,f

以下FFTの描画範囲を狭めて、挿入図として描画するために入力の信号と周波数の範囲を指定しています。上で書いていますが、FFTの時間範囲は[2.5,2.8]を取っています。
※この周波数分割を適当に実施して順に計算してFFTを実施し全体をFFT強度vs freqで横軸を時間で強度分布を描画したものがSTFTとなります。
今回はSTFTについては、signal.stft関数を使うこととします。

sig1=sig[int(2.5*fr):int(2.8*fr)]
fn1=int(fn*0.3/5)
freq,Pyy,fsk,Psk,f=FFT(sig1,fn1,fr)

以下でFFTの描画を実施しています。基本的に前回のFFTの図と同様に実施しています。
描画範囲は、上記の関数で得られたピーク振動数とピーク値を利用して
[min(fsk)*0.9, max(fsk)*1.1, 0,max(Pyy_abs)*1.5]
で描画しています。

Pyy_abs=np.abs(Pyy)
axes7.plot(f,Pyy_abs)
axes7.axis([min(fsk)*0.9, max(fsk)*1.1, 0,max(Pyy_abs)*1.5])  #0.5, 2
axes7.grid(True)
axes7.set_xscale('log')
axes7.set_ylim(0,max(Pyy_abs)*1.5)
axes7.set_title('{}'.format(np.round(fsk[:len(fsk)],decimals=1))+'\n'+'{}'.format(np.round(Psk[:len(fsk)],decimals=4)),size=10)  #グラフのタイトルにピーク周波数とピーク強度を出力する
axes7.plot(fsk[:len(fsk)],Psk[:len(fsk)],'ro')  #ピーク周波数、ピーク強度の位置に〇をつける
# グラフにピークの周波数をテキストで表示
for i in range(len(fsk)):
    axes7.annotate('{0:.1f}'.format(fsk[i]),
                 xy=(fsk[i], Psk[i]),
                 xytext=(10, 20),
                 textcoords='offset points',
                 arrowprops=dict(arrowstyle="->",connectionstyle="arc3,rad=.2")
                )

以下でSTFTを求めます。
npersegの値が以前不確定性原理として記述した周波数と時間の不確定性を決める変数です。
前回は周波数の不確定性を許していましたが、今回はどちらかというと周波数の精度を上げるために10000としています。
※この変数を変えた時の見え方の差をあとで示します

nperseg=10000
f, t, Zxx = signal.stft(sig, fr, nperseg=nperseg)

得られるZxxは虚数なので、以下のように絶対値を求めます。
また、最大値を計算しておきます。

rr=np.abs(Zxx)
max_rr=np.max(rr)

以下で強度を色階調で示す二次元分布で描画できます。
また、colorbarをaxes4の位置に描画しています。

img=axes3.pcolormesh(t, f, rr, vmin=0, vmax=max_rr,cmap="hsv")
axes3.set_yscale('log')
axes3.set_ylim([20, 20000])
axes3.set_xlim([0, 5])
axes3.set_title('STFT Magnitude')
axes3.set_ylabel('Frequency [Hz]')
axes3.set_xlabel('Time [sec]')
fig.colorbar(img, cax=axes4)

同様に、同じく強度を色階調の二次元分布を描画します。
ただし、今回は上記の信号vs時間の挿入図に合わせて描画範囲を[1,1.1]をしています。
同じく、colorbarをaxes6の位置に描画しています。
ここでは色コードを上記ではcmap="hsv"、以下ではcmap="jet"を使ってみました。
※ここは変える必要はありません

img1 = axes5.pcolormesh(t, f, rr, vmin=0, vmax=max_rr,cmap="jet")
axes5.set_yscale('log')
axes5.set_ylim([20, 20000])
axes5.set_xlim([1, 1.1])
fig.colorbar(img1, cax=axes6)

あとは、描かれた画像を保存して終了です。

plt.pause(1)
plt.savefig("./stft/stft_fig_"+wavfile+"_"+str(nperseg)+".jpg", dpi=200)
plt.close()

結果

上記のコードで上にしめしたような図が得られます。
さらに、前回と同様に基準振動400Hzに以下の周波数を重畳すると以下のような図が得られました。前回と比較してよく見ると、FFTの分解能が悪いのに気が付くと思います。これはこの時間範囲でFFTを実施しているため、不確定性原理が表れ始めていると言えます。
list_=[0.01,0.1,0.5,1,2,5,10,50,100,200,300,350,400.1,400.2,400.5,401,402,405,410,420,450,500,600,800,1000,1200,1400,1600]

また、不確定原理で示したようにnpersegを変更すると以下のような図が得られました。
こちらは、入力信号は変更せず、単にSTFTの分析する時間範囲のみを変更しています。
これは不確定性原理が陽に見えていると言えます。
list_=[100,200,500,1000,2000,5000,10000,20000,50000,100000,200000,500000]

この図から、だいたいこの周波数領域で音声を分析する場合は、
nperseg=2000-10000
辺りが適切な範囲であることが分かります。
これは、時間にすると、fr=44100とすると
nperseg/fr=0.045-0.23sec
程度がいいことが分かります。

コードは以下におきました

・AudioAutoencoder/stft.py

まとめ

・比較的簡単なコードで正しくSTFTできるようになった
・近接した周波数のSTFTもできることが分かった
・時間的に変化する音声を分解して分析できることが分かる
・音声などを分析するには、fr=44100/secのサンプリング周波数でnperseg=2000-10000程度で分析するのがよさそうである

・STFTをFFTから作成する
・フォルマント合成を本格的に実施する
・waveletで分析する