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

import numpy as np
import FFT_lib as fft

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

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

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

#ホワイトノイズの標準偏差を定義（ホワイトノイズ、ピンクノイズに使用）
noise_std = 5.0

# 2つの周波数が混合した信号を生成
signal1, signal2, mixed_signal = fft.generate_testData(t, freq1, freq2, amp1, amp2)

"""------以下、ノイズは3つのうち1つを選択すること------"""
#混合信号にホワイトノイズを追加（必要に応じてコメントアウト）
noise, mixed_signal_withNoise = fft.add_whiteNoise(mixed_signal, noise_std)

#混合信号にピンクノイズを追加（必要に応じてコメントアウト）
#noise, mixed_signal_withNoise = fft.add_pinkNoise(sampling_rate, mixed_signal, noise_std)

#混合信号にインパルスノイズを追加（必要に応じてコメントアウト）
#noise, mixed_signal_withNoise = fft.add_impulseNoise(sampling_rate, mixed_signal)
"""-----------------------------------------------"""

print(f"テストデータが作成されました：y(t) = {amp1}*sin(2π{freq1}t) + {amp2}*sin(2π{freq2}t) + w(t)")

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

#ノイズ成分のみのfFT解析を実行
frequency_noise, amplitude_noise, fft_result_noise = fft.execute_fft(noise, sampling_rate)

#ノイズ成分のみのFFT解析の結果をプロット
fft.plot_fftResult(t, t_short, noise, frequency_noise, amplitude_noise)
print("ノイズ成分のみのFFT解析結果をプロットしました。")

# テストデータ全体のFFT解析を実行
frequency, amplitude, fft_result = fft.execute_fft(mixed_signal_withNoise, 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_withNoise, frequency, amplitude)
print("FFT解析結果をプロットしました。")

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq, rfft, rfftfreq, irfft
# グラフを美しくする設定
plt.rcParams['figure.figsize'] = (12, 8)
plt.rcParams['font.size'] = 11

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

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

    return signal1, signal2, mixed_signal

def add_whiteNoise(original_singal, noise_std):
    """元の信号にホワイトノイズを追加
    Takes  : original_singal: 元の信号 [ ]
             noise_std: ホワイトノイズの標準偏差 [ ]
    Returns: white_noise: 生成したホワイトノイズ [ ]
             noisy_signal: ホワイトノイズを追加した信号 [ ]
    """
    #ホワイトノイズを生成
    white_noise = np.random.normal(0, noise_std, size=original_singal.shape)
    #ホワイトノイズを追加
    noisy_signal = original_singal + white_noise

    return white_noise, noisy_signal

def add_pinkNoise(sampling_rate,original_singal, noise_std):
    """元の信号にピンクノイズを追加
    Takes  : sampling_rate: サンプリング周波数 [Hz]
             original_singal: 元の信号 [ ]
             noise_std: ピンクノイズの標準偏差 [ ]
    Returns: pink_noise: 生成したピンクノイズ [ ]
             noisy_signal: ピンクノイズを追加した信号 [ ]
    """
    N = len(original_singal)
    #ホワイトノイズを生成
    white_noise = np.random.normal(0, noise_std,N)
    # FFT
    W = rfft(white_noise)
    # 周波数軸を作成
    freqs = rfftfreq(N, 1/sampling_rate)
    # 1/f のスペクトルにスケーリング（0Hzを除外）
    W[1:] /= np.sqrt(freqs[1:])
    
    # 逆FFTで時系列に戻す
    pink_noise = irfft(W)
    
    # ピンクノイズを追加
    noisy_signal = original_singal + pink_noise

    return pink_noise, noisy_signal

def add_impulseNoise(sampling_rate, original_singal):
    """元の信号にインパルスノイズを追加
    Takes  : sampling_rate: サンプリング周波数 [Hz]
             original_singal: 元の信号 [ ]
    Returns: impulse_noise: 生成したインパルスノイズ [ ]
             noisy_signal: インパルスノイズを追加した信号 [ ]
    """
    N = len(original_singal)
    # インパルスノイズを生成
    impulse_noise = np.zeros(N)
    impulse_noise[0] = sampling_rate

    # インパルスノイズを追加
    noisy_signal = original_singal + impulse_noise
    
    return impulse_noise, noisy_signal
    

def plot_testData_withNoise(t, t_short, signal1, signal2, white_noise, 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_max_abs = max(np.abs(signal1[mask].max()), np.abs(signal2[mask].max()), np.abs(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_max_abs, y_max_abs)

    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_max_abs, y_max_abs)

    #ノイズ成分
    axes[1, 0].plot(t_short_list, white_noise[mask], 'r-', linewidth=2, label='60Hz')
    axes[1, 0].set_title(r'$w(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_max_abs, y_max_abs)

    # 混合信号
    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) + w(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_max_abs, y_max_abs)

    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
    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

    return frequencies_positive, amplitude_single_sided, fft_result


def plot_fftResult(t, t_short, mixed_signal, frequency, amplitude):
    """FFT結果を美しく可視化（片側振幅スペクトラム）
    Takes  : t: 時間軸
             t_short: 最初の何秒間を表示するか
             mixed_signal: 3つの信号を混合したもの
             frequency: 周波数軸 [Hz]
             amplitude: 片側振幅スペクトラム [ ]
    Returns: None
    """
    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()