Pythonで信号処理(1)：フーリエ変換

#はじめに
Pythonで行う信号処理の勉強資料。

#使用するライブラリー
numpy, matplotlib
scipyのfftpack

#実行内容

1. 特定周波数の正弦波を生成する。(numpy)
2. FFTを実行する。(scipy)
3. 実行結果をプロットする。(matplotlib)

#プログラムコード

##1.正弦波の生成部分

N = 2**20 # data number
dt = 0.0001 # data step [s]
f1 = 5 # frequency[Hz]
A1 = 1 # Amplitude
p1 = 90*pi/180  # phase [rad]

#波形形成
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) 

##2.フーリエ変換部分

# 離散フーリエ変換&規格化
yf = fft(y)/(N/2) 

# y : numpy 配列
# N : サンプリング数

全体コード

fft.py

import numpy as np
from scipy.fftpack import fft
import matplotlib.pyplot as plt

#円周率π
pi = np.pi

# parameters

N = 2**20 # data number
dt = 0.0001 # data step [s]
f1 = 5 # frequency[Hz]
A1 = 1 # Amplitude
p1 = 90*pi/180 # 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)

#フーリエ変換
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.xlabel("frequency")
plt.ylabel("amplitude")

plt.tight_layout()
plt.savefig("01")

plt.show()

#3.実行結果
上：時間領域のグラフ：位相が90度ずれていることを確認。
下：周波数領域のグラフ：5Hzの周波数が確認。

#まとめ
1.PythonでFFTを実行した。

##感想
これまで信号処理には、matlabを使ってきたが、pythonのscipyのほうも簡単。

#参考資料