Help us understand the problem. What is going on with this article?

【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[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

Kuma_T
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away