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

```python
import numpy as np
import FFT_lib as fft

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

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

# サンプルデータの周波数と振幅を設定
freq1, freq2 = 6.3, 25        # Hz
amp1, amp2 = 3.0, 1.5         # 振幅

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

# 窓関数を生成
window = fft.generate_hanning_window(N)

# 生成した信号に窓関数を適用
mixed_signal_window = fft.apply_window(mixed_signal, window)

# 生成した信号と窓関数をプロット
t_short = 2.0  # 最初の2.0秒だけ表示
fft.plot_testData(t, 2.0, signal1, signal2, mixed_signal, window, mixed_signal_window)
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, mixed_signal, frequency, amplitude)

# 窓関数をかけた場合のFFT解析を実行
frequency_window, amplitude_window, fft_result_window = fft.execute_fft(
    mixed_signal_window, sampling_rate
)

# FFT解析の結果をプロット（窓あり・なしを比較）
fft.plot_fftResult_window(
    t,
    mixed_signal,
    mixed_signal_window,
    frequency,
    amplitude,
    frequency_window,
    amplitude_window
)
print("窓適用後の FFT をプロットしました。")

# === 畳み込み定理の検証 ===
freqs, Y_windowed, conv_scaled, X, H = fft.verify_convolution_theorem(mixed_signal, window, sampling_rate)

# 誤差の簡易評価
abs_diff = np.max(np.abs(Y_windowed - conv_scaled))
rel_diff = abs_diff / (np.max(np.abs(Y_windowed)) + 1e-12)
print(f"最大絶対差 (FFT(windowed) vs (1/N)*circular_conv): {abs_diff:.3e}")
print(f"相対差（最大振幅比）: {rel_diff:.3e}")

# 比較プロット
fft.plot_convolution_verification(freqs, Y_windowed, conv_scaled, X, H, sampling_rate, show_phase=False)

print("畳み込み定理の検証プロットを表示しました。")
```

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

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

# 3つの周波数が混合した信号を生成
def generate_testData(t, freq1, freq2, amp1, amp2):
    """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)

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

    return signal1, signal2, mixed_signal


def plot_testData(t, t_short, signal1, signal2, mixed_signal, window, mixed_signal_window):
    """時間域での信号をプロット"""
    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(),
                mixed_signal[mask].min(), mixed_signal_window[mask].min())
    y_max = max(signal1[mask].max(), signal2[mask].max(),
                mixed_signal[mask].max(), mixed_signal_window[mask].max())

    axes[0, 0].plot(t_short_list, signal1[mask], 'b-', linewidth=2)
    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_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)
    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, window, 'r-', linewidth=2)
    axes[1, 0].set_title(r'$window: h(t)$', fontweight='bold')
    axes[1, 0].grid(True, alpha=0.3)
    axes[1, 0].set_xlim(0, t_short)

    axes[1, 1].plot(t_short_list, mixed_signal[mask], 'purple', linewidth=2, label='y(t)')
    axes[1, 1].plot(t_short_list, mixed_signal_window[mask],
                    color='orange', linestyle='--', linewidth=2, label='y(t)h(t)')
    axes[1, 1].set_title(r'$y = a_1 \sin(\omega_1 t) + a_2 \sin(\omega_2 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)
    axes[1, 1].legend()

    plt.tight_layout()
    plt.show()


def execute_fft(mixed_signal, sampling_rate):
    """信号にFFTを適用して周波数解析"""
    N = len(mixed_signal)

    fft_result = fft(mixed_signal)
    frequencies = fftfreq(N, 1 / sampling_rate)
    amplitude_spectrum = np.abs(fft_result)

    pos_mask = frequencies >= 0
    frequencies_positive = frequencies[pos_mask]
    amplitude_positive = amplitude_spectrum[pos_mask]

    amplitude_single_sided = amplitude_positive * 2 / N
    amplitude_single_sided[0] = amplitude_positive[0] / N

    return frequencies_positive, amplitude_single_sided, fft_result


def generate_hanning_window(N):
    """ハニング窓を生成"""
    return np.hanning(N)


def apply_window(signal, window):
    """窓を信号に適用"""
    if len(signal) != len(window):
        raise ValueError("signal and window must have the same length")
    return signal * window


def plot_fftResult(t, mixed_signal, frequency, amplitude):
    """FFT結果を可視化"""
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 10))

    ax1.plot(t, mixed_signal, 'navy', linewidth=2)
    ax1.set_title('Time Domain Signal', fontweight='bold')
    ax1.grid(True, alpha=0.3)

    ax2.plot(frequency, amplitude, 'navy', linewidth=2, marker='o', markersize=3)
    ax2.set_title('Frequency Domain Signal', fontweight='bold')
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()


def plot_fftResult_window(t, mixed_signal, mixed_signal_window,
                          frequency, amplitude,
                          frequency_window, amplitude_window):
    """窓あり・なしFFT比較"""
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 10))

    ax1.plot(t, mixed_signal, '--', label='y(t)')
    ax1.plot(t, mixed_signal_window, label='y(t)h(t)')
    ax1.grid(True, alpha=0.3)
    ax1.legend()

    ax2.