昨夜は③音声データ入力編を記事にしたが、そこでリアルタイムのスペクトログラムのイメージが出来たので、早速予定を変更してこちらを先に実施することとした。
###とりあえずの目標
以下のとおり、
※リンクされているものは記事にしたものであり、その記事を前提知識として書いて行くので、参照すると理解し易いと思います。
①Scipy環境を作る
・環境の確認
②不確定原理について
・FFT変換・逆変換してみる；時間軸が消える
・STFT変換・逆変換してみる；窓関数について
・wavelet変換・逆変換してみる；スペクトログラムの解釈と不確定原理
③音声や地震データや株価や、。。。とにかく一次元の実時系列データに応用する
音声データ入力編
④FFTからwavelet変換まで簡単にたどってみる（上記以外のちょっと理論）
⑤二次元データに応用してみる
⑥天体観測データに応用してみる
⑦リアルタイムにスペクトログラムしてみる
###今回やったこと
「⑦リアルタイムにスペクトログラムしてみる」をやってみた
・コード解説
・得られたスペクトログラム
###・コード解説
import pyaudio
import wave
from scipy.fftpack import fft, ifft
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from swan import pycwt
中間の出力用にWAVE_OUTPUT_FILENAMEを使う。
測定Rateは44.1kHzではなく、半分の22.1kHzとしている。
録音時間は少し長めの5sとしている。
また、録音はモノラルとした。
CHUNK = 1024
FORMAT = pyaudio.paInt16 # int16型
CHANNELS = 1 # 1;monoral 2;ステレオ
RATE = 22100 # 22.1kHz 44.1kHz
RECORD_SECONDS = 5 # 5秒録音
WAVE_OUTPUT_FILENAME = "output2.wav"
リアルタイム感を出すために、Whileで継続するようにした。
測定開始の合図として、print("* recording")を入れている。
s=1
while True:
p = pyaudio.PyAudio()
stream = p.open(format=FORMAT,
channels=CHANNELS,
rate=RATE,
input=True,
frames_per_buffer=CHUNK)
print("* recording")
5s測定を継続して、dataにため込む。
frames = []
for i in range(0, int(RATE / CHUNK * RECORD_SECONDS)):
data = stream.read(CHUNK)
frames.append(data)
print("* done recording")
stream.stop_stream()
stream.close()
p.terminate()
一度、ファイルに出力する。
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 = wr.getnchannels()
width = wr.getsampwidth()
fr = wr.getframerate()
fn = wr.getnframes()
fs = fn / fr
origin = wr.readframes(wr.getnframes())
data = origin[:fn]
wr.close()
元の音声画像を出力する。
amp = max(data)
plt.figure(figsize=(12, 10))
# ステレオ前提 > monoral
sig = np.frombuffer(data, dtype="int16") /32768.0
t = np.linspace(0,fs, fn/2, endpoint=False)
plt.subplot(311)
plt.xlim([0, 5])
plt.plot(t, sig)
plt.pause(1)
#plt.close()
縦軸横軸が一致するように、ファクターで合わせている。
plt.pcolormeshはcmap='hsv'とすることにより、俄然見やすくなった。
※plt.pcolormesh(X, Y, Z, cmap='hsv')という書式で、(X,Y)面にZ(X,Y)の大きさに応じてcmapのルールで色付けする。
nperseg = 256
sig = np.frombuffer(data, dtype="int16")/32768.0
print('fs',fs, fn)
f, t, Zxx = signal.stft(sig, fs=fs*fn/50, nperseg=nperseg)
plt.subplot(312)
plt.pcolormesh(t, 5*f, np.abs(Zxx), cmap='hsv')
plt.ylim([10*f[1], 10*f[-1]])
plt.xlim([0, 5])
plt.title('STFT Magnitude')
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.yscale('log')
plt.pause(1)
#plt.close()
元画像、STFTのスペクトログラム、そしてFFT画像を出力する。
freq =fft(sig,int(fn/2))
Pyy = np.sqrt(freq*freq.conj())*2/fn
f = np.arange(int(fn/2))
plt.subplot(313)
plt.plot(f,Pyy)
plt.axis([100, max(f)/2, 0,0.00005]) #max(Pyy)])
plt.xscale('log')
plt.pause(1)
#plt.close()
plt.savefig('figure'+str(s)+'.png')
wavelet変換のスペクトログラムも同時に出力しようとしたが、出来ないので今回は別途出力することとした。
※やはりパラメータは出力画像が一致するようにファクターで合わせている。
Fs = 1/0.0002
omega0 = 5 #0.2 #1 #2 #8
# (1) Freqを指定してcwt
x = np.linspace(0,fs, fn/2, endpoint=False)
freqs=np.arange(10,2000,2.5)
r=pycwt.cwt_f(sig,freqs,Fs,pycwt.Morlet(omega0))
rr = np.abs(r)
plt.rcParams['figure.figsize'] = (10, 6)
fig = plt.figure()
ax1 = fig.add_axes([0.1, 0.75, 0.7, 0.2])
ax2 = fig.add_axes([0.1, 0.1, 0.7, 0.60], sharex=ax1)
ax3 = fig.add_axes([0.83, 0.1, 0.03, 0.6])
ax1.plot(x, sig, 'k')
img = ax2.imshow(np.flipud(rr), extent=[0, 5,100, 20000], aspect='auto', cmap='hsv')
twin_ax = ax2
twin_ax.set_yscale('log')
twin_ax.set_xlim(0, 5)
twin_ax.set_ylim(100, 20000)
ax2.tick_params(which='both', labelleft=False, left=False)
twin_ax.tick_params(which='both', labelleft=True, left=True, labelright=False)
fig.colorbar(img, cax=ax3)
plt.pause(1)
plt.savefig('figure_'+str(s)+'.png')
s += 1
###・得られたスペクトログラム
おはよう
開けゴマ
###まとめ
・リアルタイムなFFT,STFT,そしてwavelet変換ができるようになった
・STFTまではかなりリアルタイム的に動くが、wavelet変換は待ち時間が発生している。
・全体に高速化を実施して文字通りなリアルタイム測定ができるようにしたい
・スペクトログラム的に、STFT変換とwavelet変換の違いが見えにくい状況である