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

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

ということで、録音された音声は以下のとおりです。
「おはよう」
![ohayo.jpg](https://qiita-image-store.s3.amazonaws.com/0/233744/b072b866-2659-3729-cfcb-f9203c0191b9.jpeg)
「開けゴマ」
![hirakegoma.jpg](https://qiita-image-store.s3.amazonaws.com/0/233744/96155208-2959-23be-cb86-fa256c911b38.jpeg)
###・音声の分離
上記の音声は何度か「おはよう」と「開けゴマ」を繰り返しているので、これを分離します。
コードは以下のとおりです。
流れとしては、
①録音の全体を読み
②生データ、STFT、そしてFFTを表示する
③切り出したい時間軸を選ぶ
を繰り返す。
使うのは以下のとおり、見慣れたものです。

```py
#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ファイルを読み込み、諸量を出力します。

```py
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を表示します。
※リアルタイム・スペクトログラムのコードの表示とほぼ同じです。
最後にグラフを保存しています。

```py
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ファイルを読み込むこととしました。

```py
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)としています。

```py
    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にデータを格納します。

```py
    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)で生データ等を出力します。

```py
    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ファイルに書き出しています。

```py
    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を表示して、終了します。

```py
    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になっています
![ohayo_out.jpg](https://qiita-image-store.s3.amazonaws.com/0/233744/afa69763-0856-b15d-77b2-4fdc634ce86b.jpeg)
![ohayo_out.jpg](https://qiita-image-store.s3.amazonaws.com/0/233744/19fc2f5d-ad22-5763-9ebc-671cecd6b257.jpeg)
![ohayo_out.jpg](https://qiita-image-store.s3.amazonaws.com/0/233744/f4bc0cb4-73cd-8bc3-746d-62709f8e2276.jpeg)
「開けゴマ」も以下のように切り出せました。
![hirakegoma_out.jpg](https://qiita-image-store.s3.amazonaws.com/0/233744/7ab29272-e081-d778-d899-6dd528f89e58.jpeg)
![hirakegoma_out.jpg](https://qiita-image-store.s3.amazonaws.com/0/233744/c2f96e1d-099b-5361-1b3e-85b4266e5217.jpeg)
![hirakegoma_out.jpg](https://qiita-image-store.s3.amazonaws.com/0/233744/24100ba9-2d8d-d584-a2b3-0de7a7d57610.jpeg)
![hirakegoma_out.jpg](https://qiita-image-store.s3.amazonaws.com/0/233744/710b483b-a698-95c7-58b8-5571c034f851.jpeg)
![hirakegoma_out.jpg](https://qiita-image-store.s3.amazonaws.com/0/233744/b31560b8-e4a6-ae18-2755-fba8fbf176e4.jpeg)
###まとめ
・音声の録音データを分離するアプリを作成した
・綺麗に分離することができた。
・時間的には変化していても、時間軸の縮尺を変えると似たような表示を得た

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