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

[前回の記事](https://qiita.com/MuAuan/items/85b077640901dbb29514)は、とりあえずのリアルタイム・スペクトログラムのアプリについて記載したが、今回はこれを高速化して、実際の様子をGifアニメーションで掲載することとする。
###とりあえずの目標
以下のとおり、
※リンクされているものは記事にしたものであり、その記事を前提知識として書いて行くので、参照すると理解し易いと思います。
①[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)
###やったこと
・高速化；時報をキャッチする
・Gifアニメーションで見る;前回のGifアニメーションの技術を利用して見える化する
###・高速化；時報をキャッチする
この体感を伝えるのに適切なものって。。。ということで、ちょっと音声おかしいけど、参考に掲載の時報をキャッチしてその性能を見せたいと思います
【参考】
・[http://domisan.sakura.ne.jp/timesignal/timesignal_opus.html](インターネット時報(Opus版))

高速化コードの工夫点を以下示します。
※少し長めですが、全部載せます
使っているものは見慣れてきました。

```py
# -*- 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
```
一応、録音する長さCHUNKは1024を単位としてNで可変とし、サンプリングRateも変更してみました。
よくあるアプリだと、このCHUNK＝1024を固定にしてFor文などで回してappendしていますが、一挙に取得する方が高速化には寄与しているようです。
また、参考にパラメータ一覧が掲載されていますが、```pyaudio.paInt16 #16 bit int```を利用します。
【参考】
・[PyAudio Documentation](https://people.csail.mit.edu/hubert/pyaudio/docs/)

```py
N=50
CHUNK=1024*N
RATE=44100 #11025 #22050  #44100
CHANNELS = 1             # 1;monoral 2;ステレオ-
p=pyaudio.PyAudio()
WAVE_OUTPUT_FILENAME = "output.wav"
FORMAT = pyaudio.paInt16 #int16型

stream=p.open(	format = pyaudio.paInt16,
		channels = 1,
		rate = RATE,
		frames_per_buffer = CHUNK,
		input = True,
		output = True) # inputとoutputを同時にTrueにする
```
グラフ描画のためのパラメータです。
ax1(音声生データ描画)、ax2(STFT描画)そしてax3(FFT描画)を描画します。
それぞれ軸のラベルと軸の属性（対数など）そして描画位置を指定しています。

```py
s=1
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()
```
while stream.is_active():で回していますが、for文で回数指定して回す方がいい場合もあると思います。
※今回は終了はctrl-cで終了しています
ax1とax3は上書きすると汚くなるので、一度消去して描画するようにしています。
また、経過時間を表示するようにしました。

```py
while stream.is_active():

    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)
```
別途、標準出力にstream.readにかかる感度時間と不感時間を表示するようにしました。

```py
    start_time=time.time()
    print(start_time-stop_time)
    input = stream.read(CHUNK)
    stop_time=time.time()
```
今回もinputのままで描画することはできなかったので、一度ファイル出力して再度読み込んでいます。
※ここがpython版としては最後の高速化として残っていますが、ご容赦ください
　print文は不感時間が計測できれば高速化のためには不要です（今回残していますが。。）

```py
    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
    print("fn,fs",fn,fs,stop_time-start_time)

    origin = wr.readframes(wr.getnframes())
    data = origin[:fn]
    wr.close()
```
sigは一応規格化しています。
また、グラフの表示範囲が変わると画像が揺れるので、固定としました。
※data読込はint16です。

```py
    sig = np.frombuffer(data, dtype="int16")  /32768.0
    t = np.linspace(0,fs, fn/2, endpoint=False)
    ax1.set_ylim(-0.0075,0.0075)
    ax1.set_xlim(0,fs)
    ax1.plot(t, sig)
```
[STFTは不確定性原理の影響](https://qiita.com/MuAuan/items/504160465e83e556dd3e)も考えて、nperseg=1024としました。
※sampling時間fs*tと周波数範囲f*RATE/N/200に対して、Zxxを表示します。

```py
    nperseg = 1024
    f, t, Zxx = signal.stft(sig, fs=fn/2, nperseg=nperseg)
    ax2.pcolormesh(fs*t, f*RATE/N/200, np.abs(Zxx), cmap='hsv')
    ax2.set_xlim(0,fs)
    ax2.set_ylim(200,20000)
    ax2.set_yscale('log')
```
FFTの周波数範囲が難しいところですが、以下のようにすると自動的に周波数範囲を計算してくれます。

```py
    freq =fft(sig,int(fn/2))
    Pyy = np.sqrt(freq*freq.conj())*2/fn
    f = np.arange(200,RATE*100/50,(RATE*100/50-200)/int(fn/2)) #RATE11025,22050;N50,100
    ax3.set_ylim(0,0.000075)
    ax3.set_xlim(200,20000)
    ax3.set_xlabel('Freq[Hz]')
    ax3.set_ylabel('Power')
    ax3.set_xscale('log')
    ax3.plot(f,Pyy)
```
ある意味、これが高速化で一番肝ですが、表示時間を0.01にしています。
そして、画像を一枚一枚Dirに保存します。
※pngでもOKですが、ファイルが小さいのでjpgにしています

```py
    plt.pause(0.01)
    plt.savefig('out_voice5s/figure'+str(s)+'.jpg')
    s += 1
```
以下のコメントを外すと、スピーカーから音声が出力されてちょっと面白いですが、通常は切っておきます。

```py
    #output = stream.write(input)
```
以下のコードもfor文のときは有効ですが、ctrl-cでは使っていません。

```py
stream.stop_stream()
stream.close()
p.terminate()

print( "Stop Streaming")
```
###・Gifアニメーションで見る
;前回のGifアニメーションの技術を利用して見える化する
以下が時報のGifアニメーションです。
※時報なので一定周波数にピークが出ています⇒４４０Hzだとすると周波数が10倍になっているので修正する必要がある
![out_voice_5s52.gif](https://qiita-image-store.s3.amazonaws.com/0/233744/902563c2-1245-0019-3a83-ba5d105829c1.gif)
さらに時報ではなく、音声「おはようー開けゴマ」を約5s間隔でキャッチしてみる。
※表示は半分位の時間で表示している
![out_voice_5s20.gif](https://qiita-image-store.s3.amazonaws.com/0/233744/e062365c-7f98-d7b0-f680-6349bde13bf8.gif)
もっと高速化してみる、0.5s間隔でキャッチ。
![out_voice_05s32.gif](https://qiita-image-store.s3.amazonaws.com/0/233744/91e0ace4-e948-82b9-2e33-f107a318d2eb.gif)
###まとめ
・リアルタイム・スペクトログラムを高速化した
・Gifアニメーションで見える化した

・【バグ】時報の周波数が１０倍ずれているようなので後日修正する予定である
###おまけ
音声約5s間隔キャッチの時間出力を以下に示す。

```
>python pyaudio_real_time.py
0.013960123062133789
fn,fs 51200 4.6439909297052155 4.517699480056763
C:\Users\user\Anaconda3\lib\site-packages\numpy\core\numeric.py:531: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
0.7087409496307373
fn,fs 51200 4.6439909297052155 3.941642999649048
0.3313875198364258
fn,fs 51200 4.6439909297052155 4.3084776401519775
0.36440467834472656
fn,fs 51200 4.6439909297052155 4.275268316268921
0.36235713958740234
fn,fs 51200 4.6439909297052155 4.287429094314575
0.37495923042297363
fn,fs 51200 4.6439909297052155 4.265578508377075
0.38704490661621094
fn,fs 51200 4.6439909297052155 4.262617826461792
0.4144895076751709
fn,fs 51200 4.6439909297052155 4.226044178009033
0.4065864086151123
fn,fs 51200 4.6439909297052155 4.232300758361816
0.437361478805542
fn,fs 51200 4.6439909297052155 4.213600158691406
0.46187925338745117
fn,fs 51200 4.6439909297052155 4.178019046783447
0.47734522819519043
fn,fs 51200 4.6439909297052155 4.1725780963897705
0.48210930824279785
fn,fs 51200 4.6439909297052155 4.157163858413696
0.5058131217956543
fn,fs 51200 4.6439909297052155 4.134882688522339
0.5725016593933105
fn,fs 51200 4.6439909297052155 4.077564001083374
0.5428729057312012
fn,fs 51200 4.6439909297052155 4.097465753555298
0.5508785247802734
fn,fs 51200 4.6439909297052155 4.098431348800659
0.5634572505950928
fn,fs 51200 4.6439909297052155 4.076149940490723
0.6018803119659424
fn,fs 51200 4.6439909297052155 4.038370370864868
0.5855240821838379
fn,fs 51200 4.6439909297052155 4.064543724060059
Stop Streaming
```