Pythonで学習するFFT解析（高速フーリエ変換解析）➀

import numpy as np
import FFT_lib as fft

# サンプリング設定
sampling_rate = 1000  # 1000Hz（1秒に1000回測定）
duration = 2.0        # 2秒間
N = int(sampling_rate * duration)  # 総サンプル数

# 時間軸を作成
t = np.linspace(0, duration, N, endpoint=False)

# サンプルデータの周波数と振幅を設定
freq1, freq2, freq3 = 10, 30, 60         # Hz
amp1, amp2, amp3 = 3.0, 1.5, 0.8         # 振幅

# 3つの周波数が混合した信号を生成
signal1, signal2, signal3, mixed_signal = fft.generate_testData(t, freq1, freq2, freq3, amp1, amp2, amp3)
print(f"テストデータが作成されました：y(t) = {amp1}*sin(2π{freq1}t) + {amp2}*sin(2π{freq2}t) + {amp3}*sin(2π{freq3}t)")

# 生成した信号をプロット
t_short = 0.5  # 最初の0.5秒だけ表示
fft.plot_testData(t, 0.5, signal1, signal2, signal3, mixed_signal)
print("時間領域での信号をプロットしました。")

# FFT解析を実行
frequency, amplitude, fft_result = fft.execute_fft(mixed_signal, sampling_rate)
print("FFT解析を実行しました。")
print(f"周波数分解能: {frequency[1] - frequency[0]:.2f} Hz")
print(f"最大周波数: {frequency[-1]:.0f} Hz （ナイキスト周波数）")

# FFT解析の結果をプロット
fft.plot_fftResult(t, t_short, mixed_signal, frequency, amplitude)

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq

# グラフを美しくする設定
plt.rcParams['figure.figsize'] = (12, 8)
plt.rcParams['font.size'] = 11

# 3つの周波数が混合した信号を生成
def generate_testData(t,freq1,freq2,freq3,amp1,amp2,amp3):
    """3つの周波数が混合した信号を生成
    Takes  : 
        t: 時間軸 [sec]
        freq1, freq2, freq3: 周波数1, 2, 3 [Hz]
        amp1, amp2, amp3: 振幅1, 2, 3 [ ]
    Returns:
        signal1, signal2, signal3: 個別の正弦波信号 [ ]
        mixed_signal: 3つの信号を混合したもの [ ]
    """
    signal1 = amp1 * np.sin(2 * np.pi * freq1 * t) 
    signal2 = amp2 * np.sin(2 * np.pi * freq2 * t) 
    signal3 = amp3 * np.sin(2 * np.pi * freq3 * t) 

    # 3つの信号を混合
    mixed_signal = signal1 + signal2 + signal3

    return signal1, signal2, signal3, mixed_signal

def plot_testData(t, t_short, signal1, signal2, signal3, mixed_signal):
    """時間域での信号をプロット
    Takes  : t: 時間軸
             t_short: 最初の何秒間を表示するか
             signal1, signal2, signal3: 個別の正弦波信号
             mixed: 3つの信号を混合したもの
    Returns: None
    """
    fig, axes = plt.subplots(2, 2, figsize=(15, 10))

    mask = t <= t_short
    t_short_list = t[mask]

    # 縦軸範囲を全プロットで統一
    y_min = min(signal1[mask].min(), signal2[mask].min(), signal3[mask].min(), mixed_signal[mask].min())
    y_max = max(signal1[mask].max(), signal2[mask].max(), signal3[mask].max(), mixed_signal[mask].max())

    # 個別信号
    axes[0, 0].plot(t_short_list, signal1[mask], 'b-', linewidth=2, label='10Hz')
    axes[0, 0].set_title(r'$y_1 = a_1 \sin(\omega_1 t)$', fontweight='bold')
    axes[0, 0].grid(True, alpha=0.3)
    axes[0, 0].set_ylabel('amplitude')
    axes[0, 0].set_xlim(0, t_short)
    axes[0, 0].set_ylim(y_min, y_max)
    
    axes[0, 1].plot(t_short_list, signal2[mask], 'g-', linewidth=2, label='30Hz')
    axes[0, 1].set_title(r'$y_2 = a_2 \sin(\omega_2 t)$', fontweight='bold')
    axes[0, 1].grid(True, alpha=0.3)
    axes[0, 1].set_xlim(0, t_short)
    axes[0, 1].set_ylim(y_min, y_max)

    axes[1, 0].plot(t_short_list, signal3[mask], 'r-', linewidth=2, label='60Hz')
    axes[1, 0].set_title(r'$y_3 = a_3 \sin(\omega_3 t)$', fontweight='bold')
    axes[1, 0].grid(True, alpha=0.3)
    axes[1, 0].set_ylabel('amplitude')
    axes[1, 0].set_xlim(0, t_short)
    axes[1, 0].set_ylim(y_min, y_max)

    # 混合信号
    axes[1, 1].plot(t_short_list, mixed_signal[mask], 'purple', linewidth=2)
    axes[1, 1].set_title(r'$y = a_1 \sin(\omega_1 t) + a_2 \sin(\omega_2 t) + a_3 \sin(\omega_3 t)$', fontweight='bold')
    axes[1, 1].grid(True, alpha=0.3)
    axes[1, 1].set_xlim(0, t_short)
    axes[1, 1].set_ylim(y_min, y_max)

    plt.tight_layout()
    plt.show()

def execute_fft(mixed_signal, sampling_rate):
    """信号にFFTを適用して周波数解析
    Takes  : signal: 入力信号 [ ]
                sampling_rate: サンプリング周波数 [Hz]
    Returns: frequencies_positive: 正の周波数成分 [Hz]
                amplitude_single_sided: 片側振幅スペクトラム [ ]
                fft_result: FFTの複素数結果 [ ]
    """
    N = len(mixed_signal)

    # FFT実行
    fft_result = fft(mixed_signal)

    # 周波数軸を作成
    frequencies = fftfreq(N, 1 / sampling_rate)

    # 振幅スペクトラムを計算（複素数の絶対値）
    amplitude_spectrum = np.abs(fft_result)

    # 正の周波数のみ取得（実信号の場合、負の周波数は冗長）
    positive_freq_mask = frequencies >= 0  # DC成分(0Hz)も含む
    frequencies_positive = frequencies[positive_freq_mask]
    amplitude_positive = amplitude_spectrum[positive_freq_mask]

    # 片側振幅スペクトラムに変換（DC成分以外は2倍）
    amplitude_single_sided = amplitude_positive * 2 / N
    amplitude_single_sided[0] = amplitude_positive[0] / N  # DC成分のみ2倍しない

    return frequencies_positive, amplitude_single_sided, fft_result

def plot_fftResult(t, t_short, mixed_signal, frequency, amplitude):
    """FFT結果を美しく可視化（片側振幅スペクトラム）"""
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 10))

    # 上段：時間領域
    time_mask = t <= t_short
    ax1.plot(t[time_mask], mixed_signal[time_mask], 'navy', linewidth=2)
    ax1.set_title('Time Domain Signal (Test Data)', fontsize=14,fontweight='bold')
    ax1.set_xlabel('Time [sec]')
    ax1.set_ylabel('Amplitude')
    ax1.set_xlim(0, t_short)
    ax1.grid(True, alpha=0.3)

    # 下段：周波数領域
    ax2.plot(frequency, amplitude, 'crimson', linewidth=2)
    ax2.set_title('Frequency Domain Signal (FFT result)', fontsize=14, fontweight='bold')
    ax2.set_xlabel('Frequency [Hz]')
    ax2.set_ylabel('Amplitude')
    ax2.set_xlim(0, frequency[-1])
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()