Python
scipy
numpy
matplotlib
フーリエ変換

【Python】ただひたすらにフーリエ変換【備忘録】

ただひたすらに(1次元)フーリエ変換.

意外と何やってるかわからなくなるので備忘録も兼ねて.

フーリエ変換(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波

01.png



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

もちろん周波数も振幅も変わらない

02.png



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波

周波数も振幅も変わらない

03.png



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

04.png

周波数については,ほぼ正確にわかる

一方で振幅は誤差がありそう.



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

波の形状は複雑だが周波数は求められている.

05.png



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

06.png



Code-07

# parameters

y = np.zeros(N)*dt
y[50000] =5
y[N-50000]=5

yf = ifft(y)*(N/2) # 逆フーリエ変換&規格化


逆フーリエ変換

周波数・振幅はあってる.cos波になった.

07.png

周波数軸から時間軸に直すときは,(対称性を考慮して)2点打つ.



Code-08

# parameters

y = np.zeros(N)*dt
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点打つ.

08.png



Code-09

from scipy.special import jv

y = jv(0, t)


第一種ベッセル関数を試してみる.

完全に周期的な関数ではないので,ピークは確認されたが,すそを引いた周波数の結果が得られた.

09.png



Code-10

y = t - np.floor(t)


$y = x - |x|$

これも一種の周期関数

「整数値」ごとに周期的なピーク

10.png



Code-11

y = np.ceil(t) - np.round(t)


$y = |x+1| - round(x)$

奇数ごとに周期的なピーク

11.png



Code-12

y = np.zeros(N)*dt

y[::10000] =1

周期1ごとに値1

ちょっと変わった感じの関数が得られた.

12.png


周波数解析に使われるフーリエ変換をひたすらにやりました。

いつか誰かの役に立てばと思います。

以上です。


参考文献

Discrete Fourier transforms (scipy.fftpack)

https://docs.scipy.org/doc/scipy/reference/fftpack.html