ただひたすらに（１次元）フーリエ変換．
意外と何やってるかわからなくなるので備忘録も兼ねて．
フーリエ変換（Fourier transform; FT）は、時間(or空間)軸と周波数軸の変換．
基本のコードは以下．
Code-01
import numpy as np
from scipy.fftpack import fft
import matplotlib.pyplot as plt
# parameters
N = 2**20 # data number
dt = 0.0001 # data step [s]
f1, f2 = 5, 8 # frequency[Hz]
A1, A2 = 5, 0 # Amplitude
p1, p2 = 0, 0 # phase
t = np.arange(0, N*dt, dt) # time
freq = np.linspace(0, 1.0/dt, N) # frequency step
y = A1*np.sin(2*np.pi*f1*t + p1) + A2*np.sin(2*np.pi*f2*t + p2)
yf = fft(y)/(N/2) # 離散フーリエ変換&規格化
plt.figure(2)
plt.subplot(211)
plt.plot(t, y)
plt.xlim(0, 1)
plt.xlabel("time")
plt.ylabel("amplitude")
plt.subplot(212)
plt.plot(freq, np.abs(yf))
plt.xlim(0, 10)
#plt.ylim(0, 5)
plt.xlabel("frequency")
plt.ylabel("amplitude")
plt.tight_layout()
plt.savefig("01")
plt.show()
周期5Hz sin波
Code-02
# parameters
N = 2**20 # data number
dt = 0.0001 # data step [s]
f1, f2 = 5, 8 # frequency[Hz]
A1, A2 = 5, 0 # Amplitude
p1, p2 = 2, 0 # phase
y = A1*np.sin(2*np.pi*f1*t + p1) + A2*np.sin(2*np.pi*f2*t + p2)
周期5 位相シフト2
もちろん周波数も振幅も変わらない
Code-03
# parameters
f1, f2 = 5, 8 # frequency[Hz]
A1, A2 = 5, 0 # Amplitude
p1, p2 = 0, 0 # phase
y = A1*np.cos(2*np.pi*f1*t + p1) + A2*np.sin(2*np.pi*f2*t + p2)
周期5 cos波
周波数も振幅も変わらない
Code-04
# parameters
f1, f2 = 2, 8 # frequency[Hz]
A1, A2 = 10, 0 # Amplitude
p1, p2 = 0, 0 # phase
y = A1*np.sin(2*np.pi*f1*t + p1) + A2*np.sin(2*np.pi*f2*t + p2)
周波数2 振幅10
周波数については，ほぼ正確にわかる
一方で振幅は誤差がありそう．
Code-05
# parameters
f1, f2 = 5, 8 # frequency[Hz]
A1, A2 = 5, 3 # Amplitude
p1, p2 = 0, 0 # phase
y = A1*np.sin(2*np.pi*f1*t + p1) + A2*np.sin(2*np.pi*f2*t + p2)
周波数5 振幅5 & 周波数8 振幅3
波の形状は複雑だが周波数は求められている．
Code-06
# parameters
f1, f2 = 5, 0 # frequency[Hz]
A1, A2 = 5, 0 # Amplitude
y = A1*np.sin(2*np.pi*f1*t)**2
周波数5 $\sin^2$波
周波数5でも，2乗したら実質の周波数は2倍の10
Code-07
# parameters
y = np.zeros(N)*dt
y[50000] =5
y[N-50000]=5
yf = ifft(y)*(N/2) # 逆フーリエ変換&規格化
逆フーリエ変換
周波数・振幅はあってる．cos波になった．
周波数軸から時間軸に直すときは，(対称性を考慮して)2点打つ．
Code-08
# parameters
y = np.zeros(N)*dt
y[int(f1/dt)] =5
y[N-int(f1/dt)] =5
y[int(f2/dt)] =3
y[N-int(f2/dt)] =3
yf = ifft(y)*(N/2) # 逆フーリエ変換&規格化
2点打つ．
Code-09
from scipy.special import jv
y = jv(0, t)
第一種ベッセル関数を試してみる．
完全に周期的な関数ではないので，ピークは確認されたが，すそを引いた周波数の結果が得られた．
Code-10
y = t - np.floor(t)
$y = x - |x|$
これも一種の周期関数
「整数値」ごとに周期的なピーク
Code-11
y = np.ceil(t) - np.round(t)
$y = |x+1| - round(x)$
奇数ごとに周期的なピーク
Code-12
y = np.zeros(N)*dt
y[::10000] =1
周期1ごとに値1
ちょっと変わった感じの関数が得られた．
周波数解析に使われるフーリエ変換をひたすらにやりました。
いつか誰かの役に立てばと思います。
以上です。
参考文献
Discrete Fourier transforms (scipy.fftpack)
https://docs.scipy.org/doc/scipy/reference/fftpack.html