FFT(Fast Fourier Transformation)
file:/Users/bob/Github/TeamNishitani/jupyter_num_calc/fft
https://github.com/daddygongon/jupyter_num_calc/tree/master/notebooks_python
cc by Shigeto R. Nishitani 2017-19
https://github.com/daddygongon/jupyter_num_calc/tree/master/notebooks_python
cc by Shigeto R. Nishitani 2017-19
FFTの応用¶
Fast Fourier Transformation(FFT)高速フーリエ変換(あるいはデジタル(離散)フーリエ変換(DFT))は，周波数分解やフィルターを初め，画像処理などの多くの分野で使われている．基本となる考え方は，直交基底による関数の内挿法である．最初にその応用例を見た後，どのような理屈でFFTが動いているかを解説する．
いくつかFFTを解説する面白い動画があります．ただ，少し高度かも．．．
- 初音ミクを三角関数で描いてみた
- [A visual introduction.
周波数分解¶
はじめの例は，周波数分解．先ずは，非整合な波を二つ用意しておく．
これを重ねあわせた波を作る．
ゆっくり変化する波に，激しく変化する波が重なっていることが読み取れる．これにFFTを掛ける その強さを求めて，周波数で表示すると，
もとの2つの周波数に対応するところにピークができているのが確認できる．広がりは，誤差のせい． logplotでも良い．
python code¶
scipyにあるfft, ifftを使う．
In [1]:
%matplotlib inline
from scipy.fftpack import fft # 古いライブラリではscipyのみで, fftpack -> fftかも
import matplotlib.pyplot as plt
import numpy as np
def func(x):
return np.sin(x/(13))+np.sin(x/(2/2))
x = np.linspace(0, 256, 256) #0から2πまでの範囲を100分割したnumpy配列
plt.plot(x, func(x), color = 'b')
plt.plot(x, np.sin(x/(13)), color = 'r', linewidth=0.8)
plt.plot(x, np.sin(x/(2/2)), color = 'r', linewidth=0.8)
plt.grid()
plt.show()
yy = func(x)
out = fft(yy)
def spectrum_power(x):
re, im = x.real, x.imag
return np.sqrt(re**2+im**2)
plt.plot(x,spectrum_power(out))
plt.xlim(0,128)
plt.show()
plt.plot(x,spectrum_power(out))
plt.xlim(0,128)
plt.yscale('log')
plt.show()
/Users/bob/opt/anaconda3/lib/python3.8/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.23.4
warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
高周波フィルター¶
次の例は，高周波フィルター．たとえば次のようなローレンツ関数を考える． $$ {\it f1}\, := \,t\mapsto 10+ \frac{40000}{380+ \left( t-128 \right) ^{2} } $$
これに(測定機器などからの)白色ノイズ(y_noise)が載ると，計測結果は次のようになる．
In [7]:
%matplotlib inline
from scipy.fftpack import fft, ifft
import matplotlib.pyplot as plt
import numpy as np
def func(x):
a,b,c,d=10, 40000, 380, 128
return a+b/(c+(x-d)**2)
xdata = np.linspace(0, 256, 256)
y = func(xdata)
y_noise = 10 * np.random.normal(size=xdata.size)
ydata = y + y_noise
plt.plot(xdata, ydata, 'b-', label='data')
plt.grid()
plt.show()
これに高周波フィルターを掛けるとノイズが消えるが，その様子を示そう．
先ずは，FFTを掛ける． これは次のようなスペクトル強度分布(spectral power distribution)(log scaleで)をもっている．
In [37]:
out = fft(ydata)
def spectral_power(x):
re, im = x.real, x.imag
return np.sqrt(re**2+im**2)
plt.plot(x,spectral_power(out))
plt.xlim(0,128)
plt.yscale('log')
plt.show()
低周波の部分に，ゆっくりとした変化を表す成分が固まっている．次のような三角フィルターを用意する．これは，低周波ほど影響を大きく