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


ただひたすらに（１次元）フーリエ変換．

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

フーリエ変換（Fourier transform; FT）は、時間(or空間)軸と周波数軸の変換．

基本のコードは以下．


```python: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](https://qiita-image-store.s3.amazonaws.com/0/192073/220990d7-438a-c6a7-e45a-5aa5f5cc5d37.png)

----


```python: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](https://qiita-image-store.s3.amazonaws.com/0/192073/22ac8388-7864-862f-9524-eca46a5045c8.png)

----

```ruby: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](https://qiita-image-store.s3.amazonaws.com/0/192073/0147a72b-8a19-5a16-4a79-767cbad14540.png)

----

```ruby: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](https://qiita-image-store.s3.amazonaws.com/0/192073/658d8738-367e-f345-63f9-b14ac4e2b6f6.png)

周波数については，ほぼ正確にわかる
一方で振幅は誤差がありそう．

----

```ruby: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](https://qiita-image-store.s3.amazonaws.com/0/192073/e9fb9431-39fa-1c2b-87f4-819a9f447884.png)

----

```ruby: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](https://qiita-image-store.s3.amazonaws.com/0/192073/30cd29a5-d52b-419e-a1bb-6c6d56159783.png)

----

```ruby:Code-07
# parameters
y = np.zeros(N)*dt
y[50000] =5
y[N-50000]=5

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

**逆フーリエ変換**
周波数・振幅はあってる．cos波になった．
![07.png](https://qiita-image-store.s3.amazonaws.com/0/192073/8fafd68c-fd19-582b-32f5-af4750084bc8.png)

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

----

```ruby: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](https://qiita-image-store.s3.amazonaws.com/0/192073/622edfb7-0635-3316-a325-48250b8112f9.png)

----

```ruby:Code-09
from scipy.special import jv

y = jv(0, t)
```

第一種ベッセル関数を試してみる．
完全に周期的な関数ではないので，ピークは確認されたが，すそを引いた周波数の結果が得られた．
![09.png](https://qiita-image-store.s3.amazonaws.com/0/192073/73a7d90b-7299-3e82-a6b6-d5c124982d5c.png)


----


```ruby:Code-10
y = t - np.floor(t)
```

$y = x - |x|$
これも一種の周期関数
「整数値」ごとに周期的なピーク
![10.png](https://qiita-image-store.s3.amazonaws.com/0/192073/ced28a62-8f70-d39a-a14a-4d69cb9d2ec8.png)

----


```ruby:Code-11
y = np.ceil(t) - np.round(t)
```

$y = |x+1| - round(x)$
奇数ごとに周期的なピーク
![11.png](https://qiita-image-store.s3.amazonaws.com/0/192073/aa9b3864-c2c1-b7df-d6ea-f6566d6e9599.png)


----


```ruby:Code-12
y = np.zeros(N)*dt
y[::10000] =1
```

周期1ごとに値1
ちょっと変わった感じの関数が得られた．
![12.png](https://qiita-image-store.s3.amazonaws.com/0/192073/385ba840-4046-814b-15a7-53dd4a9f6fb2.png)


----

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



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

以上です。


# 参考文献
Discrete Fourier transforms (scipy.fftpack)
https://docs.scipy.org/doc/scipy/reference/fftpack.html