FFT、STFTそしてwavelet変換はいろいろデータ処理する時はお世話になりそうな技術です。
しかも、それなりに慣れてこないと使いこなすことはできない技術の一つだと思います。
そして、今回はScipy使ってやろうと考えていたが、環境もそれなりに難易度高いと思います。
###とりあえずの目標
以下のとおりの予定で数回に分けて書こうと思う。
※理解が違っていたら、都度ここも直します。。。
①Scipy環境を作る
・環境の確認
②不確定原理について
・FFT変換・逆変換してみる；時間軸が消える
・STFT変換・逆変換してみる；窓関数について
・wavelet変換・逆変換してみる；スペクトログラムの解釈と不確定原理
③音声や地震データや株価や、。。。とにかく一次元の実時系列データに応用する
④FFTからwavelet変換まで簡単にたどってみる（上記以外のちょっと理論）
⑤二次元データに応用してみる
⑥天体観測データに応用してみる
⑦リアルタイムにスペクトログラムしてみる
###今回やったこと
①Scipy環境を作る ・環境の確認
・FFT変換の例題
・STFT変換の例題
・wavelet変換の例題
###①Scipy環境を作る
まず、Scipyの関数が良く出来ている。
だから、浮気せずに順に掲載されているものを動かしていこうと思う。
ところが、STFT変換が通らない。
ということで、例によって環境構築から。
まず、最初にやるべきはScipyのバージョンチェック
STFTをやりたいのであれば、'0.19.0'以上であることが必要です。
>python
>>> import scipy
>>> scipy.__version__
'0.18.0'
【参考】
・How to check the version of scipy
・module 'scipy.signal' has no attribute 'stft'
「stft func is added in the SciPy 0.19.0, check this post to upgrade your module to the newsest version.」
ということで、ウワンは以下のようにconda環境でScipyのインストールをやりました。
※実はWindows環境でも'pip install scipy==0.19.1'としましたが、途中でエラーが出てインストールできませんでした
>conda install scipy==0.19.1
Solving environment: done
## Package Plan ##
environment location: C:\Users\tosio\Anaconda3
added / updated specs:
- scipy==0.19.1
The following packages will be downloaded:
package | build
---------------------------|-----------------
openssl-1.0.2p | hfa6e2cd_0 5.4 MB
scikit-learn-0.20.0 | py35heebcf9a_1 5.1 MB
numpy-1.11.3 | py35h4a99626_4 3.4 MB
intel-openmp-2019.1 | 144 1.7 MB
icc_rt-2019.0.0 | h0cc432a_1 9.4 MB
scipy-0.19.1 | py35h1e2829e_3 12.6 MB
mkl-2018.0.3 | 1 178.1 MB
ca-certificates-2018.03.07 | 0 155 KB
numexpr-2.6.8 | py35h9ef55f4_0 129 KB
------------------------------------------------------------
Total: 215.9 MB
The following NEW packages will be INSTALLED:
ca-certificates: 2018.03.07-0
icc_rt: 2019.0.0-h0cc432a_1
intel-openmp: 2019.1-144
The following packages will be UPDATED:
mkl: 11.3.3-1 --> 2018.0.3-1
numexpr: 2.6.1-np111py35_0 --> 2.6.8-py35h9ef55f4_0
numpy: 1.11.1-py35_1 --> 1.11.3-py35h4a99626_4
openssl: 1.0.2j-vc14_0 --> 1.0.2p-hfa6e2cd_0
scikit-learn: 0.17.1-np111py35_1 --> 0.20.0-py35heebcf9a_1
scipy: 0.18.1-np111py35_0 --> 0.19.1-py35h1e2829e_3
Proceed ([y]/n)? y
Downloading and Extracting Packages
openssl-1.0.2p | 5.4 MB | ############################################################################ | 100%
scikit-learn-0.20.0 | 5.1 MB | ############################################################################ | 100%
numpy-1.11.3 | 3.4 MB | ############################################################################ | 100%
intel-openmp-2019.1 | 1.7 MB | ############################################################################ | 100%
icc_rt-2019.0.0 | 9.4 MB | ############################################################################ | 100%
scipy-0.19.1 | 12.6 MB | ############################################################################ | 100%
mkl-2018.0.3 | 178.1 MB | ############################################################################ | 100%
ca-certificates-2018 | 155 KB | ############################################################################ | 100%
numexpr-2.6.8 | 129 KB | ############################################################################ | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done
####・環境の確認
環境確認のために、関数の説明ページにある例題を動かしてみます。
####・FFT変換の例題
from scipy.fftpack import fft, ifft
import numpy as np
import matplotlib.pyplot as plt
x = np.arange(5)
y =fft(x)
xy = ifft(y)
s=np.allclose(ifft(fft(x)), x, atol=1e-15) # within numerical accuracy.
print(x,xy,s)
結果は以下のとおり戻っていて、np.allclose関数がTrueを返しています。
※例題がつまらないのでおまけにそれらしいのを追記しました
>python fft_ex.py
[0 1 2 3 4] [ 0.+0.j 1.+0.j 2.+0.j 3.+0.j 4.+0.j] True
####・STFT変換の例題
FFT変換の例題と比較すると、格段に面白い例題になっています。
※訳文付きで載せます
【参考】
・scipy.signal.istft
#scipy.signal.istft example
#https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.istft.html
#
import numpy as np #added by author
from scipy import signal
import matplotlib.pyplot as plt
#Generate a test signal, a 2 Vrms sine wave at 50Hz corrupted by 0.001 V**2/Hz of white noise sampled at 1024 Hz.
#テスト信号、1024 Hzでサンプリングされた0.001 V ** 2 / Hzのホワイトノイズで破損した50 Hzの2 Vrmsの正弦波を生成します
fs = 1024
N = 10*fs
nperseg = 512
amp = 2 * np.sqrt(2)
noise_power = 0.001 * fs / 2
time = np.arange(N) / float(fs)
carrier = amp * np.sin(2*np.pi*50*time)
noise = np.random.normal(scale=np.sqrt(noise_power),
size=time.shape)
x = carrier + noise
#Compute and plot the STFT’s magnitude.
#STFTの振幅を計算してプロットします
f, t, Zxx = signal.stft(x, fs=fs, nperseg=nperseg)
plt.figure()
plt.pcolormesh(t, f, np.abs(Zxx), vmin=0, vmax=amp)
plt.ylim([f[1], f[-1]])
plt.title('STFT Magnitude')
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.yscale('log')
plt.show()
#Zero the components that are 10% or less of the carrier magnitude, then convert back to a time series via inverse STFT
#キャリア振幅の10％以下の成分をゼロにしてから、逆STFTを介して時系列に変換し直す
Zxx = np.where(np.abs(Zxx) >= amp/10, Zxx, 0)
_, xrec = signal.istft(Zxx, fs)
#Compare the cleaned signal with the original and true carrier signals.
#きれいにされた信号を元のそして本当の搬送波信号と比較
plt.figure()
plt.plot(time, x, time, xrec, time, carrier)
plt.xlim([2, 2.1])
plt.xlabel('Time [sec]')
plt.ylabel('Signal')
plt.legend(['Carrier + Noise', 'Filtered via STFT', 'True Carrier'])
plt.show()
#Note that the cleaned signal does not start as abruptly as the original, since some of the coefficients of the transient were also removed:
#トランジェントの係数の一部も削除されているため、クリーン化された信号は元の信号ほど急激には開始されません。
plt.figure()
plt.plot(time, x, time, xrec, time, carrier)
plt.xlim([0, 0.1])
plt.xlabel('Time [sec]')
plt.ylabel('Signal')
plt.legend(['Carrier + Noise', 'Filtered via STFT', 'True Carrier'])
plt.show()
【補足】このグラフの見どころは、緑の線が最初徐々に振幅が大きくなっていて実際のキャリアとずれているところです。
###・wavelet変換の例題
分かりやすさのために二行追記してみた。
つまり周波数によって支配的な振動項が変わっているのが一目で分かる。
※このグラフは縦：横＝周波数：時間のグラフで、時間経過とともにスペクトラムが変化する様子を表しています。
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from scipy import signal
import matplotlib.pyplot as plt
t = np.linspace(-1, 1, 200, endpoint=False)
sig = np.cos(2 * np.pi * 7 * t) + signal.gausspulse(t - 0.4, fc=2)
#**以下二行追加**
plt.plot(t,sig)
plt.pause(3)
#******
widths = np.arange(1, 31)
cwtmatr = signal.cwt(sig, signal.ricker, widths)
plt.imshow(cwtmatr, extent=[-1, 1, 1, 31], cmap='PRGn', aspect='auto',
vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())
plt.show()
###まとめ
・Scipyの環境を構築して、FFT、STFT、そしてCWTのそれぞれのサンプルを動かしてみた
・さて、いろいろできるので次回が楽しみ。
・とりあえずの目標を書いていこうと思う
###おまけ
・FFT変換の例
スケールはともかく、時間領域（vs振幅）の図から、周波数領域（vsパワー）の図に変換される。どちらも、周波数や時間（それぞれ共役量）の記憶がなくなる（積分される）。
from scipy.fftpack import fft, ifft
import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(0, 100, 10000, endpoint=False)
sig = np.cos(2 * np.pi * 4 * t)+np.cos(2 * np.pi * 8 * t) + np.cos(2 * np.pi * 15 * t) #*np.exp(-0.1*t) *5
plt.plot(t, sig)
plt.axis([0, 5, -5, 5])
plt.show()
freq =fft(sig,1024)
Pyy = np.abs(freq)/1025 #freq*freq.conj(freq)/1025
f = np.arange(1024)
plt.plot(f,Pyy)
plt.axis([0, 200, 0, 1])
plt.show()
以下のとおり、逆変換で戻りました。
t1=np.linspace(0, 10, 1024, endpoint=False)
sig1=ifft(freq)
plt.plot(t1, sig1)
plt.axis([0, 5, -5, 5])
plt.show()
