FFTの周波数分解能の違いをグラフ描画してみた【Python・numpy】

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 = 50.25 #正弦波の周波数
gain = 1 #正弦波の振幅
offset = 0 #直流成分
t = 4 #データ上限（秒）

#時間軸データ #t秒まで、dt秒周期で作成
time = np.arange(0, t, dt)
# 信号を生成（周波数10の正弦波+周波数20の正弦波+ランダムノイズ）
f = gain * np.sin(2*np.pi*F1*time)

#周波数軸
Ffreq = np.fft.fftfreq(n, d = dt) #データ数、サンプリング周期を指定する。
#高速フーリエ変換
FFT_f = np.fft.fft(f[:n])
FFT_f = abs(FFT_f * 2 / n) #交流成分の正規化（2倍してサンプル数で割り算して絶対値を取る）
FFT_f[0] = abs(FFT_f[0] / 2) #直流成分の正規化（サンプル数（2）で割り算して絶対値を取る）

#FFT最大値と周波数の取得
Fmax = max(FFT_f[:int(n/2.56)]) #FFTの最大値取得
maxHz = Ffreq[np.argmax(FFT_f[:int(n/2.56)])] #FFTの最大値周波数位置の取得
fftmax = "peak:" + str(round(Fmax, 2)) + " " + "@" + str(maxHz) + "Hz"

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

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

#データ及び軸項目の作成
#正弦波のプロット
plt.subplot(3,1,1) #サブプロットの作成、位置
plt.title("時刻歴グラフ　" + "FFT分析データ数：最初の" + str(n) + "データ " + str(t) + "sec波形", loc='left') #グラフタイトルの設定、位置
plt.plot(time,f) #グラフデータの指定（x, y）
plt.xlim(0,0.1) #横軸の最小値、最大値
plt.grid() #目盛線の追加
plt.xlabel("Time[sec]") #横軸ラベル
plt.ylabel("振幅[-]") #縦軸ラベル

#FFTのプロット
plt.subplot(3,1,2) #サブプロットの作成、位置
plt.title("パワースペクトルグラフ　" + fftmax, loc='left') #グラフタイトルの設定、位置
plt.plot(Ffreq[:int(n/2.56)], FFT_f[:int(n/2.56)]) #グラフデータの指定（x, y）
plt.xlim(0, Frange) #横軸の最小値、最大値
plt.grid() #目盛線の追加
plt.xlabel("Frequency[Hz]") #横軸ラベル
plt.ylabel("振幅[-]") #縦軸ラベル

#FFTのプロット
plt.subplot(3,1,3) #サブプロットの作成、位置
plt.title("パワースペクトルグラフ　" + fftmax, loc='left') #グラフタイトルの設定、位置
plt.plot(Ffreq[:int(n/2.56)], FFT_f[:int(n/2.56)], marker="o") #グラフデータの指定（x, y）
plt.xticks(np.arange(0, Frange + 1, step = 1)) #軸目盛りを指定値から等間隔に分割
plt.xlim(45, 55) #横軸の最小値、最大値
plt.grid() #目盛線の追加
plt.xlabel("Frequency[Hz]") #横軸ラベル
plt.ylabel("振幅[-]") #縦軸ラベル