【Scipy】FFT、STFTとwavelet変換で遊んでみた♬～⑦リアルタイム・スペクトログラム；完成したよ♬

前回の記事は、リアルタイム・スペクトログラムのアプリの高速化について記載したが、今回はさらに高速化して、しかも前回の周波数が間違っていたバグを解消した。。。
###とりあえずの目標
以下のとおり、
※リンクされているものは記事にしたものであり、その記事を前提知識として書いて行くので、参照すると理解し易いと思います。
①Scipy環境を作る
　・環境の確認
②不確定原理について
・FFT変換・逆変換してみる；時間軸が消える
・STFT変換・逆変換してみる；窓関数について
・wavelet変換・逆変換してみる；スペクトログラムの解釈と不確定原理
③音声や地震データや株価や、。。。とにかく一次元の実時系列データに応用する
音声データ入力編
④FFTからwavelet変換まで簡単にたどってみる（上記以外のちょっと理論）
⑤二次元データに応用してみる
⑥天体観測データに応用してみる
⑦リアルタイムにスペクトログラムしてみる
前回の高速化
###コードは以下に置きました
・Scipy-Swan/pyaudio_realtime_last.py
###やったこと
・さらなる高速化とトリガー拾ってサンプリング開始
・パラメータを変更したときの時報のGifアニメーション
###・さらなる高速化とトリガー拾ってサンプリング開始
まず、以下のおまけに掲載のコードでinputとdataの中身が異なるのかというのを検証した。
結論から言うと、inputとdataは全く同じであった。
ということで、わざわざ出力してから再度読み込むなどということはしなくても処理ができることが分かった。
※コード見ればわかるだろうと思うので、説明は省略しています

ということで、サンプリング開始のためのトリガー関数start_mesure()というのを以下のとおり作成しました。
※ちょっと解説付けます

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

サンプリングに必要なものだけを定義します。

def start_measure():
    CHUNK=1024
    RATE=44100 #11025 #22050  #44100
    p=pyaudio.PyAudio()
    input = []
    stream=p.open(format = pyaudio.paInt16,
                  channels = 1,
                  rate = RATE,
                  frames_per_buffer = CHUNK,
                  input = True) 
    input =stream.read(CHUNK)
    sig = np.frombuffer(input, dtype="int16")/32768.0

sigが一定レベル（今は0.001）を超えたら、この関数終了します。

    while True:
        if max(sig) > 0.001:
            break
        input =stream.read(CHUNK)
        sig = np.frombuffer(input, dtype="int16")/32768.0
        #print(max(sig))

サンプリングがすごくシンプルになり、input =stream.read(CHUNK)だけになりました。
一応、関数を抜ける前にこのサンプリングはcloseしておきます。

    stream.stop_stream()
    stream.close()
    p.terminate()
    return input

利用は以下のとおり。

input=start_measure()
print(input)

###サンプリング本体とFFTなどの本体合わせた全体のコード
以下のとおりです。。
※上記と重複している部分は省略します

# -*- coding:utf-8 -*-

import pyaudio
...

def start_measure():
    CHUNK=1024
...
    return

ここから実際のサンプリング開始です。
Nでサンプリングの長さを変更します。
RATEは三種類変更可です。

N=50
CHUNK=1024*N
RATE=44100 #11025 #22050  #44100
p=pyaudio.PyAudio()

stream=p.open(format = pyaudio.paInt16,
              channels = 1,
              rate = RATE,
              frames_per_buffer = CHUNK,
              input = True)

グラフ表示は、ax1;生データ、ax2;STFT、そしてax3;FFTのグラフです。

fig = plt.figure(figsize=(12, 10))
ax1 = fig.add_subplot(311)
ax1.set_xlabel('Time [sec]')
ax1.set_ylabel('Signal')
ax2 = fig.add_subplot(312)
ax2.set_ylabel('Freq[Hz]')
ax2.set_xlabel('Time [sec]')
ax3 = fig.add_subplot(313)
ax3.set_xlabel('Freq[Hz]')
ax3.set_xscale('log')
ax3.set_ylabel('Power')

時間取得して、サンプリング時間、不感時間などを取得しました。
※本番では不要です

start=time.time()
stop_time=time.time()
stp=stop_time

これは重要で、NやRATEを変更したとき自動的に生データ、FFTやSTFTの軸制御するためにfr, fn, fsを定義します。

fr = RATE
fn=51200*N/50  #*RATE/44100
fs=fn/fr
print(fn,fs)

ここから本体。
毎回、start_measure()で閾値超えると、サンプリング開始です。

for s in range(30):
    start_measure()

生データとFFTのグラフは、毎回再描画します。
※軸などは残して測定値だけ描画したかったのですが、そういうのが見つかりませんでした。やろうとすると全部書かないとなので今回はやめときます

    fig.delaxes(ax1)
    fig.delaxes(ax3)
    ax1 = fig.add_subplot(311)
    ax1.set_title('passed time; {:.2f}(sec)'.format(time.time()-start))
    ax1.set_xlabel('Time [sec]')
    ax1.set_ylabel('Signal')
    ax3 = fig.add_subplot(313)

時間計測は本番では不要ですが、あっても割と速いです。

    input = []
    start_time=time.time()
    input = stream.read(CHUNK)
    stop_time=time.time()
    print(stop_time-start_time)

生データのグラフ表示です。

    sig = np.frombuffer(input, dtype="int16")  /32768.0
    t = np.linspace(0,fs, fn, endpoint=False)
    ax1.set_ylim(-0.0075,0.0075)
    ax1.set_xlim(0,fs)
    ax1.plot(t, sig)

STFTの計算及びグラフ表示です。
※$N$や$RATE$を変更したときにうまく表示するようにスケールしてますが、$f_r,fn$に集約したので、ここではな～んだという位シンプルになりました

    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(0,fs)
    ax2.set_ylim(20,20000)
    ax2.set_yscale('log')

FFTの計算及びグラフ表示です。
※$N$や$RATE$を変更したときにうまく表示するようにスケールしてます。ここではax3.plotは若干調整が入ります。

    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*RATE/44100,Pyy)

plt.pause(0.01)で表示していますが、この0.01は最適化していません。
そして、ぱらぱらGifアニメーションのために画像を一枚一枚保存します。

    plt.pause(0.01)
    plt.savefig('out_jihou_test/figure'+str(s)+'.jpg')

for文が完了すると終了します。

stream.stop_stream()
stream.close()
p.terminate()
print( "Stop Streaming")

###パラメータを変更したときの時報のGifアニメーション
時報を徹底的に分析して、周波数のバグを上記のコードで解消しました。
【参考】
・インターネット時報

※基準周波数がオーディオや短波とか。。。なんかすぐ用意できませんでした
ということで、時報をパラメータを変更して取得したので掲載しておきます。
出力は1秒おきにピピピ。。。と鳴って、１０秒ごとにピーンとなり、１分おきに時報でボボボボーンという感じでなっています。
※まあ、具体的には上のサイトをのぞいてみてください









###まとめ
・さらなる高速化を実施した
・周波数のバグを解消してFFTとSTFTのリアルタイム・スペクトログラムが完成した
・トリガーでサンプリング開始するようにした

・音声認識、DL導入を実施したいと思う
###おまけ
コードの下に出力例を置いています。
※iput,data,そしてframesの出力例
inputとdataは全く同じframesは[]がついているだけ異なる。

# -*- coding:utf-8 -*-

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

def start_measure():
    CHUNK=1024
    RATE=44100 #11025 #22050  #44100
    CHANNELS = 1             # 1;monoral 2;ステレオ-
    p=pyaudio.PyAudio()
    WAVE_OUTPUT_FILENAME = "output1.wav"
    FORMAT = pyaudio.paInt16 # int16型
    input = []
    stream=p.open(format = pyaudio.paInt16,
                  channels = 1,
                  rate = RATE,
                  frames_per_buffer = CHUNK,
                  input = True,
                  output = True) # inputとoutputを同時にTrueにする
    input =stream.read(CHUNK)
    frames = []
    frames.append(input)
    wf = wave.open(WAVE_OUTPUT_FILENAME, 'wb')
    wf.setnchannels(CHANNELS)
    wf.setsampwidth(p.get_sample_size(FORMAT))
    wf.setframerate(RATE)
    wf.writeframes(b''.join(frames))
    wf.close()

    wavfile = WAVE_OUTPUT_FILENAME
    wr = wave.open(wavfile, "rb")
    ch = CHANNELS #wr.getnchannels()
    width = p.get_sample_size(FORMAT) #wr.getsampwidth()
    fr = RATE  #wr.getframerate()
    fn = wr.getnframes()
    fs = fn / fr

    origin = wr.readframes(wr.getnframes())
    data = origin[:fn]
    wr.close()
    sig = np.frombuffer(data, dtype="int16")/32768.0
    sig1 = np.frombuffer(input, dtype="int16")/32768.0
    #print(max(sig))
    while True:
        if max(sig1) > 0.001:
            break
        input =stream.read(CHUNK)
        frames = []
        frames.append(input)
        wf = wave.open(WAVE_OUTPUT_FILENAME, 'wb')
        wf.setnchannels(CHANNELS)
        wf.setsampwidth(p.get_sample_size(FORMAT))
        wf.setframerate(RATE)
        wf.writeframes(b''.join(frames))
        wf.close()

        wavfile = WAVE_OUTPUT_FILENAME
        wr = wave.open(wavfile, "rb")
        ch = CHANNELS #wr.getnchannels()
        width = p.get_sample_size(FORMAT) #wr.getsampwidth()
        fr = RATE  #wr.getframerate()
        fn = wr.getnframes()
        fs = fn / fr
        
        origin = wr.readframes(wr.getnframes())
        data = origin[:fn]
        wr.close()
        sig = np.frombuffer(data, dtype="int16")/32768.0
        sig1 = np.frombuffer(input, dtype="int16")/32768.0
        
        print(max(sig),max(sig1))
    stream.stop_stream()
    stream.close()
    p.terminate()
    return data,input,frames

data,input,frames=start_measure()
print("data=",data,"input=",input,"frames=",frames)

以下出力です。

>python pyaudio_test.py
data= b'\x02\x00\x03\x00\x02\x00\x01\x00\x01\x00\x01\x00\x01\x00\x03\x00\x06\x00\t\x00\x0c\x00\x0e\x00\x10\x00\x10\x00\x10\x00\x10\x00\x0f\x00\x0f\x00\x0e\x00\x0e\x00\x0e\x00\x0e\x00\x0c\x00\x0b\x00\t\x00\x05\x00\x01\x00\xfe\xff\xfb\xff\xfb\xff\xfb\xff\xfd\xff\xff\xff\x01\x00\x01\x00\x00\x00\xff\xff\xfb\xff\xf7\xff\xf4\xff\xf1\xff\xf0\xff\xf0\xff\xf2\xff\xf3\xff\xf5\xff\xf7\xff\xf8\xff\xf9\xff\xfa\xff\xfa\xff\xfb\xff\xfc\xff\xfe\xff\xff\xff\．．．\xfd\xff\xf8\xff\xf2\xff\xee\xff\xea\xff\xe9\xff\xe9\xff\xea\xff\xec\xff\xed\xff\xee\xff\xef\xff\xef\xff\xf0\xff\xf1\xff\xf2\xff\xf3\xff\xf5\xff\xf6\xff\xf7\xff\xf9\xff\xfc\xff\xff\xff\x03\x00\x06\x00\x08\x00\x0b\x00\x0c\x00\x0f\x00\x10\x00\x14\x00\x17\x00\x1a\x00\x1d\x00\x1d\x00\x1d\x00\x1c\x00\x19\x00\x16\x00\x15\x00\x14\x00\x14\x00\x16\x00\x16\x00\x16\x00\x16\x00\x13\x00\x0f\x00\x0b\x00\x06\x00\x02\x00\xfd\xff\xfa\xff\xf7\xff\xf4\xff\xf2\xff\xf0\xff\xee\xff\xec\xff\xea\xff\xe8\xff\xe6\xff\xe5\xff\xe5\xff\xe6\xff\xe7\xff\xea\xff\xec\xff\xef\xff\xf3\xff' 
input= b'\x02\x00\x03\x00\x02\x00\x01\x00\x01\x00\x01\x00\x01\x00\x03\x00\x06\x00\t\x00\x0c\x00\x0e\x00\x10\x00\x10\x00\x10\x00\x10\x00\x0f\x00\x0f\x00\x0e\x00\x0e\x00\x0e\x00\x0e\x00\x0c\x00\x0b\x00\t\x00\x05\x00\x01\x00\xfe\xff\xfb\xff\xfb\xff\xfb\xff\xfd\xff\xff\xff\x01\x00\x01\x00\x00\x00\xff\xff\xfb\xff\xf7\xff\xf4\xff\xf1\xff\xf0\xff\xf0\xff\xf2\xff\xf3\xff\xf5\xff\xf7\xff\xf8\xff\xf9\xff\xfa\xff\xfa\xff\xfb\xff\xfc\xff\xfe\xff\xff\xff\．．．\xfd\xff\xf8\xff\xf2\xff\xee\xff\xea\xff\xe9\xff\xe9\xff\xea\xff\xec\xff\xed\xff\xee\xff\xef\xff\xef\xff\xf0\xff\xf1\xff\xf2\xff\xf3\xff\xf5\xff\xf6\xff\xf7\xff\xf9\xff\xfc\xff\xff\xff\x03\x00\x06\x00\x08\x00\x0b\x00\x0c\x00\x0f\x00\x10\x00\x14\x00\x17\x00\x1a\x00\x1d\x00\x1d\x00\x1d\x00\x1c\x00\x19\x00\x16\x00\x15\x00\x14\x00\x14\x00\x16\x00\x16\x00\x16\x00\x16\x00\x13\x00\x0f\x00\x0b\x00\x06\x00\x02\x00\xfd\xff\xfa\xff\xf7\xff\xf4\xff\xf2\xff\xf0\xff\xee\xff\xec\xff\xea\xff\xe8\xff\xe6\xff\xe5\xff\xe5\xff\xe6\xff\xe7\xff\xea\xff\xec\xff\xef\xff\xf3\xff'
frames= [b'\x02\x00\x03\x00\x02\x00\x01\x00\x01\x00\x01\x00\x01\x00\x03\x00\x06\x00\t\x00\x0c\x00\x0e\x00\x10\x00\x10\x00\x10\x00\x10\x00\x0f\x00\x0f\x00\x0e\x00\x0e\x00\x0e\x00\x0e\x00\x0c\x00\x0b\x00\t\x00\x05\x00\x01\x00\xfe\xff\xfb\xff\xfb\xff\xfb\xff\xfd\xff\xff\xff\x01\x00\x01\x00\x00\x00\xff\xff\xfb\xff\xf7\xff\xf4\xff\xf1\xff\xf0\xff\xf0\xff\xf2\xff\xf3\xff\xf5\xff\xf7\xff\xf8\xff\xf9\xff\xfa\xff\xfa\xff\xfb\xff\xfc\xff\xfe\xff\xff\xff\．．．\xfd\xff\xf8\xff\xf2\xff\xee\xff\xea\xff\xe9\xff\xe9\xff\xea\xff\xec\xff\xed\xff\xee\xff\xef\xff\xef\xff\xf0\xff\xf1\xff\xf2\xff\xf3\xff\xf5\xff\xf6\xff\xf7\xff\xf9\xff\xfc\xff\xff\xff\x03\x00\x06\x00\x08\x00\x0b\x00\x0c\x00\x0f\x00\x10\x00\x14\x00\x17\x00\x1a\x00\x1d\x00\x1d\x00\x1d\x00\x1c\x00\x19\x00\x16\x00\x15\x00\x14\x00\x14\x00\x16\x00\x16\x00\x16\x00\x16\x00\x13\x00\x0f\x00\x0b\x00\x06\x00\x02\x00\xfd\xff\xfa\xff\xf7\xff\xf4\xff\xf2\xff\xf0\xff\xee\xff\xec\xff\xea\xff\xe8\xff\xe6\xff\xe5\xff\xe5\xff\xe6\xff\xe7\xff\xea\xff\xec\xff\xef\xff\xf3\xff']