【Scipy】FFT、STFTとwavelet変換で遊んでみた♬～音声分離アプリ

いよいよ、STFTを利用したアプリを作ってみよう。ということで、まずは音声の分離アプリを作ってみました。
###とりあえずの目標
以下のとおり、
※リンクされているものは記事にしたものであり、その記事を前提知識として書いて行くので、参照すると理解し易いと思います。
①Scipy環境を作る
　・環境の確認
②不確定原理について
・FFT変換・逆変換してみる；時間軸が消える
・STFT変換・逆変換してみる；窓関数について
・wavelet変換・逆変換してみる；スペクトログラムの解釈と不確定原理
③音声や地震データや株価や、。。。とにかく一次元の実時系列データに応用する
音声データ入力編
④FFTからwavelet変換まで簡単にたどってみる（上記以外のちょっと理論）
⑤二次元データに応用してみる
⑥天体観測データに応用してみる
⑦リアルタイムにスペクトログラムしてみる
高速化（バグあり）
高速化完成版
###コードは以下に置きました
・Scipy-Swan/out_wav.py
###やったこと
・録音する
・音声の分離
###録音する
③音声データ入力編で利用したコードを使って、音声を入力します。
今までの流れで「おはよう」５秒と「開けゴマ」１０秒録音しました。
今回は、リアルタイムにはこだわらないこととします。

ということで、録音された音声は以下のとおりです。
「おはよう」

「開けゴマ」

###・音声の分離
上記の音声は何度か「おはよう」と「開けゴマ」を繰り返しているので、これを分離します。
コードは以下のとおりです。
流れとしては、
①録音の全体を読み
②生データ、STFT、そしてFFTを表示する
③切り出したい時間軸を選ぶ
を繰り返す。
使うのは以下のとおり、見慣れたものです。

#coding: utf-8
import wave
import pyaudio
import matplotlib.pyplot as plt
import numpy as np
import time
import wave
from scipy.fftpack import fft, ifft
from scipy import signal

録音したwavファイルを読み込み、諸量を出力します。

def printWaveInfo(wf):
    """WAVEファイルの情報を取得"""
    print( "chn:", wf.getnchannels())
    print( "width:", wf.getsampwidth())
    print( "sampl rate:", wf.getframerate())
    print( "no. of frame:", wf.getnframes())
    print( "parames:", wf.getparams())
    print( "length(s):", float(wf.getnframes()) / wf.getframerate())

以下は、時間軸で切り出したファイルの生データ、STFTそしてFFTを表示します。
※リアルタイム・スペクトログラムのコードの表示とほぼ同じです。
最後にグラフを保存しています。

def plot_wav(filename,t1,t2):
    ws = wave.open(filename+'.wav', "rb")
    ch = ws.getnchannels()
    width = ws.getsampwidth()
    fr = ws.getframerate()
    fn = ws.getnframes()
    fs = fn / fr
    print(fn,fs)
    origin = ws.readframes(fn)
    data = origin[:]
    sig = np.frombuffer(data, dtype="int16")  /32768.0
    t = np.linspace(0,fs, fn, endpoint=False)
    fig = plt.figure(figsize=(12, 10))
    ax1 = fig.add_subplot(311)
    ax2 = fig.add_subplot(312)
    ax3 = fig.add_subplot(313)

    ax1.set_xlabel('Time [sec]')
    ax1.set_ylabel('Signal')
    ax1.set_ylim(-0.0075,0.0075)
    ax1.set_xlim(t1,t2)
    ax1.plot(t, sig)
    
    nperseg = 1024
    f, t, Zxx = signal.stft(sig, fs=fn, nperseg=nperseg)
    ax2.pcolormesh(fs*t, f/fs/2, np.abs(Zxx), cmap='hsv')
    ax2.set_xlim(t1,t2)
    ax2.set_ylim(20,20000)
    ax2.set_yscale('log')

    freq =fft(sig,int(fn))
    Pyy = np.sqrt(freq*freq.conj())*2/fn
    f = np.arange(20,20000,(20000-20)/int(fn)) #RATE11025,22050;N50,100
    ax3.set_ylim(0,0.000075)
    ax3.set_xlim(20,20000)
    ax3.set_xlabel('Freq[Hz]')
    ax3.set_ylabel('Power')
    ax3.set_xscale('log')
    ax3.plot(f*fr/44100,Pyy)
    
    plt.savefig(filename+'.jpg')
    plt.show()
    
    plt.close()  
    ws.close()

今回は、録音したwavファイルを読み込むこととしました。

if __name__ == '__main__':
    filename=input('input original filename=')
    wf = wave.open(filename+".wav", "r")
    
    printWaveInfo(wf)
    fr = wf.getframerate()
    fn = wf.getnframes()
    fs = float(fn / fr)

関数を呼び出して、画像表示します。
※最初は全体を表示したいので、plot_wav(filename,0,fs)としています。

    plot_wav(filename,0,fs)
    
    # ストリームを開く
    p = pyaudio.PyAudio()
    stream = p.open(format=p.get_format_from_width(wf.getsampwidth()),
                    channels=wf.getnchannels(),
                    rate=wf.getframerate(),
                    output=True)

    # チャンク単位でストリームに出力し音声を再生
    CHANNELS = wf.getnchannels()
    width = wf.getsampwidth()
    RATE = wf.getframerate()
    #fn = wf.readframes(wf.getnframes())
    fn = wf.getnframes()
    fs = float(wf.getnframes()) / wf.getframerate()

ここで、stream.write(data)で読み込んだ録音を出力して確認します。
同時に、framesにデータを格納します。

    frames = []
    for i in range(0, int(RATE / 1024 *fs+0.5)):
        data = wf.readframes(1024)
        frames.append(data)
        stream.write(data)
    print(int(RATE / 1024 * fs+0.5))

表示された音声の切り出す範囲をt1,t2で指定します。
plot_wav(filename,t1,t2)で生データ等を出力します。

    t1=float(input('input t1='))
    t2=float(input('input t2='))
    plot_wav(filename,t1,t2)
        
    loff = wf.getnframes()/1024 #215 #len(frames)
    print(fs,loff,loff*t1/fs,loff*t2/fs)
    #wf.close()

切り出した音声を+out.wavつけて、ファイルに保管します。
※デバッグ用にprint出力しています
wr.writeframes(b''.join(frames[s1:s2]))が肝心なコードでこれでwavファイルに書き出しています。

    wr = wave.open('wav/'+filename+'_out.wav', 'wb')
    wr.setnchannels(CHANNELS)
    wr.setsampwidth(width)  #width=2 ; 16bit
    wr.setframerate(RATE)
    s1=int(loff*(t1)/fs)
    s2=int(loff*(t2)/fs)
    print(fs,loff,s1,s2,t1,t2)
    wr.writeframes(b''.join(frames[s1:s2]))  #int(loff*t2/fs)
    #wr.close()

最後に切り出したwavファイルに基づいて、データ、STFT、そしてFFTを表示して、終了します。

    fn = wr.getnframes()
    fs = float(fn / wr.getframerate())
    print(fn,fs)
    plot_wav('wav/'+filename+'_out',0,fs)
    
    stream.close()
    p.terminate()
    exit()

###結果
「おはよう」は以下の３つのおはようを切り出しました。
※一応、同じような生データ、STFT、FFTになっています



「開けゴマ」も以下のように切り出せました。





###まとめ
・音声の録音データを分離するアプリを作成した
・綺麗に分離することができた。
・時間的には変化していても、時間軸の縮尺を変えると似たような表示を得た

・次回は音声認識に挑戦したいと思う