高速フーリエ変換FFTを使って波形データを周波数分析してみた【Numpy】

import pandas as pd
from matplotlib import pyplot as plt
import numpy as np
import japanize_matplotlib

n= 1024 #データ数
Fs = 1024 #サンプリング周波数
dt = 1/Fs #サンプリング周期
Frange = Fs/2.56 #分析周波数
F1, F2 = 50, 250 #正弦波の周波数
gain = 1 #正弦波の振幅
offset = 1 #直流成分
t = 1 #データ上限（秒）

#時間軸データ #０～１秒までを表示
time = np.arange(0, t, dt)

# 単一周波数をもつ正弦波の波形（おさらい）
f1 = gain * np.sin(2*np.pi*F1*time)

# 2種以上の周波数をもつ波形
f2 = gain * np.sin(2*np.pi*F1*time) + gain * 0.8 * np.sin(2*np.pi*F2*time)

# 2種以上の周波数及び直流成分をもつ波形
f3 = gain * np.sin(2*np.pi*F1*time) + gain * 0.8 * np.sin(2*np.pi*F2*time) + offset

# 2種以上の周波数及びランダム要素を持つ波形
f4 = gain * np.sin(2*np.pi*F1*time) + gain * 0.8 * np.sin(2*np.pi*F2*time) + gain * 0.3 * np.random.randn(n)

# 正弦波の列名称をリスト化
dlist=["time", "f1", "f2", "f3", "f4"]
# すべてのデータをDataFrameに格納
df = pd.DataFrame(list(zip(time, f1, f2, f3, f4)), columns = dlist)


#Figureの作成
fig = plt.figure(figsize=(15, 4),tight_layout=True) #figの設定、余分な余白削除

#グラフタイトルの作成
fig.suptitle("正弦波データ", weight=5, fontsize=16) #グラフタイトルの設定、太さ、サイズ

#データ及び軸項目の作成
#波形データのプロット
for i in range(len(dlist)-1): #range()はdfの列数-1を指定
　　maxd = str(round(df[dlist[i + 1]].max(), 2)) #データの最大値抽出
　　mind = str(round(df[dlist[i + 1]].min(), 2)) #データの最小値抽出
　　plt.subplot(1, len(dlist)-1, i+1) #サブプロットの作成、位置　列数は繰り返し回数を指定
　　plt.title(df.columns[i + 1] + " 最大値：" + maxd + " 最小値：" + mind, loc='left') #グラフタイトルの設定、位置
　　plt.plot(time, df.iloc[:,i + 1]) #グラフデータの指定（x, y）
　　plt.xlim(0,0.1) #横軸の最小値、最大値
　　plt.ylim(-3, 3) #横軸の最小値、最大値
　　plt.yticks(np.arange(-3, 3.1, step=0.5)) #軸目盛りを指定値から等間隔に分割
　　plt.grid() #目盛線の表示
　　plt.xlabel("Time[sec]") #横軸ラベル
　　plt.ylabel("振幅[－]") #縦軸ラベル

#周波数軸
Fdf = pd.DataFrame(np.fft.fftfreq(n, d = dt), columns = ["frequency"]) #データ数、サンプリング周期を指定する。

#高速フーリエ変換
for f in range(len(df.columns)-1): #range()はdfの列数-1を指定
　　temp = np.fft.fft(df.iloc[:, f+1]) #tempをFFT
　　temp = abs(temp * 2 / n) #交流成分の正規化（2倍してサンプル数で割り算して絶対値を取る）
　　temp[0] = abs(temp[0] / 2) #直流成分の正規化（サンプル数（2）で割り算して絶対値を取る）
　　Fdf[df.columns[f + 1]] = temp #FdfのDataFrameにFFTを追加

#print(Fdf)

#Figureの作成
fig = plt.figure(figsize=(15, 4),tight_layout=True) #figの設定、余分な余白削除

#グラフタイトルの作成
fig.suptitle("FFT パワースペクトルデータ", weight=5, fontsize=16) #グラフタイトルの設定、太さ、サイズ

#波形データのプロット
for i in range(len(Fdf.columns)-1): #range()はdfの列数-1を指定
　　plt.subplot(1, len(Fdf.columns)-1, i + 1) #サブプロットの作成、位置　列数は繰り返し回数を指定
　　plt.title(Fdf.columns[i + 1], loc='left') #グラフタイトルの設定、位置
　　plt.plot(Fdf.iloc[:, 0][0:int(n/2.56)], Fdf.iloc[:, i + 1][0:int(n/2.56)]) #グラフデータの指定（x, y）
　　plt.xlim(0, Frange) #横軸の最小値、最大値
　　plt.grid() #目盛線の表示
　　plt.xlabel("Frequency[Hz]") #横軸ラベル
　　plt.ylabel("振幅[－]") #縦軸ラベル

[0:int(n/2.56)]

plt.xlim(-1, 1) #横軸の最小値、最大値

plt.xlim(49, 50) #横軸の最小値、最大値

#FFT最大値と周波数の取得
for m in range(len(Fdf.columns)-1):
　　Fmax = max(Fdf.iloc[:, m+1]) #FFTの最大値取得
　　maxHz = np.argmax(Fdf.iloc[:, m+1]) #FFTの最大値周波数位置の取得
　　fftmax = Fdf.columns[m+1] + " peak max :" + str(round(Fmax, 2)) + " " + "@" + str(maxHz) + "Hz"
　　print(fftmax)