【Python実践】波形の「正体」を特定しよう！ FFT（フーリエ変換）を使った4つの応用例～ノイズ除去・異常検知・時間周波数解析、AI前処理～

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# 日本語フォント設定
plt.rcParams['font.family'] = 'Meiryo'

def generate_test_data(duration=1.0, fs=1000, freq=50, noise_std=0.5):
    """
    テストデータ生成関数
    FFT版では「冒頭の無音区間」は必須ではありませんが、比較のため同じデータ構造を作ります。
    """
    t = np.linspace(0, duration, int(fs * duration))
    
    # 正解データ（信号）
    clean = np.sin(2 * np.pi * freq * t) 
    
    # ノイズ
    np.random.seed(42)
    noise = np.random.normal(0, noise_std, len(t))
    
    # 入力データ
    data = clean + noise
    
    return t, clean, data

def fft_noise_filter(data, fs, threshold_ratio=0.3):
    """
    FFTによる周波数フィルタリング（閾値処理）
    
    Parameters
    ----------
    data : array_like
        入力信号データ。
    fs : int
        サンプリング周波数。
    threshold_ratio : float
        閾値の係数（0.0〜1.0）。
        「最大振幅の何割以下をノイズとみなすか」を決めます。
        0.3 なら、最大パワーの30%以下の成分をすべてカットします。

    Returns
    -------
    cleaned_signal : ndarray
        ノイズ除去後の信号データ。
    """
    n = len(data)
    
    # 1. FFTで周波数領域へ
    # 信号全体を一発で変換します
    F = np.fft.fft(data)
    freq = np.fft.fftfreq(n, d=1/fs)
    
    # 2. 振幅スペクトルを計算
    amp = np.abs(F)
    
    # 3. ノイズの閾値を決定（ここがポイント）
    # ここでは「最大ピークの〇〇%」を閾値としますが、
    # 「平均値 * 3」など統計的に決める方法もあります。
    max_amp = np.max(amp)
    threshold = max_amp * threshold_ratio
    
    # 4. フィルタリング（閾値以下の成分をゼロにする）
    # Boolean maskを作成: 閾値より大きい場所だけ True
    filter_mask = amp > threshold
    
    # マスクを適用（小さい成分はゼロになる）
    F_filtered = F * filter_mask
    
    # 5. IFFTで波形に戻す
    cleaned_signal = np.fft.ifft(F_filtered)
    
    # 実部だけを取り出して返す
    return cleaned_signal.real

def plot_results(t, clean, noisy, denoised):
    """
    結果比較プロット関数
    """
    plt.figure(figsize=(10, 8))

    # --- 1段目：正解データ ---
    plt.subplot(3, 1, 1)
    plt.plot(t, clean, color='red', linestyle='-', label='正解データ')
    plt.title("正解データ（本来の信号）")
    plt.ylabel("振幅")
    plt.grid(True, alpha=0.3)
    plt.legend(loc='upper right')

    # --- 2段目：正解データにノイズを加算 ---
    plt.subplot(3, 1, 2)
    plt.plot(t, noisy, color='gray', alpha=0.8, label='入力データ')
    plt.title("正解データにノイズを加算")
    plt.ylabel("振幅")
    plt.grid(True, alpha=0.3)
    plt.legend(loc='upper right')

    # --- 3段目：ノイズ除去結果 ---
    plt.subplot(3, 1, 3)
    plt.plot(t, denoised, color='blue', label='除去後データ')
    plt.title("ノイズ除去結果（FFTフィルタリング後）")
    plt.xlabel("時間 [s]")
    plt.ylabel("振幅")
    plt.grid(True, alpha=0.3)
    plt.legend(loc='upper right')

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    # パラメータ設定
    FS = 1000           # サンプリング周波数

    # 1. データの準備
    t, clean_sig, noisy_sig = generate_test_data(fs=FS, freq=50, noise_std=0.5)
    
    # 2. ノイズ除去の実行
    # 最大振幅の30%以下の「弱い成分」をすべてカットします
    denoised_sig = fft_noise_filter(noisy_sig, FS, threshold_ratio=0.3)
    
    # 3. 結果の描画
    plot_results(t, clean_sig, noisy_sig, denoised_sig)

import numpy as np
import matplotlib.pyplot as plt

# 日本語フォント設定
plt.rcParams['font.family'] = 'Meiryo'

def generate_data(duration=2.0, fs=1000, has_anomaly=False):
    """
    データ生成関数
    has_anomaly=True の場合、50Hzの正常音に加えて、
    300Hzの異常音が「一瞬だけ」混入したデータを生成します。
    """
    t = np.linspace(0, duration, int(fs * duration))
    
    # 1. 正常信号（50Hzの低い揺れ）
    sig = np.sin(2 * np.pi * 50 * t) * 0.5
    
    # 2. 異常信号（300Hzの高い揺れ）
    # FFTは「全体」を見るため、一瞬の混入でも「成分」として検出できます
    if has_anomaly:
        anomaly_period = (t >= 1.0) & (t <= 1.2) # 1.0~1.2秒の間だけ発生
        sig[anomaly_period] += np.sin(2 * np.pi * 300 * t[anomaly_period]) * 1.0
    
    # 3. ノイズ（全体にうっすら）
    np.random.seed(42)
    noise = np.random.normal(0, 0.2, len(t))
    
    return t, sig + noise

def fft_anomaly_score(data, fs, target_band=(200, 400)):
    """
    FFTを使って、特定周波数帯のエネルギー量を計算する関数
    """
    n = len(data)
    
    # 1. 窓関数の適用（スペクトル漏れ防止）
    window = np.hanning(n)
    data_windowed = data * window

    # 2. FFT実行
    F = np.fft.fft(data)
    freq = np.fft.fftfreq(n, d=1/fs)
    
    # 振幅スペクトル（データ点数で割って正規化）
    amp = np.abs(F) / (n / 2)
    
    # 3. 正の周波数部分だけ取り出す
    pos_mask = freq > 0
    freq = freq[pos_mask]
    amp = amp[pos_mask]
    
    # 4. 監視バンド内のインデックスを抽出
    min_f, max_f = target_band
    target_idx = np.where((freq >= min_f) & (freq <= max_f))
    
    # 5. バンド内の最大エネルギー（または合計）をスコアとする
    # ここでは「その帯域で最も強く出ている成分の高さ」をスコアにします
    if len(target_idx[0]) == 0:
        score = 0
    else:
        score = np.max(amp[target_idx])
        
    return freq, amp, score

def plot_comparison(t, normal_data, anomaly_data, normal_res, anomaly_res, threshold):
    """
    正常 vs 異常 の比較結果を描画
    """
    n_freq, n_amp, n_score = normal_res
    a_freq, a_amp, a_score = anomaly_res
    
    plt.figure(figsize=(10, 8))
    
    # --- 1段目：波形比較 ---
    plt.subplot(2, 1, 1)
    plt.plot(t, normal_data, color='green', alpha=0.6, label='正常データ')
    plt.plot(t, anomaly_data, color='red', alpha=0.6, label='異常データ(一瞬だけ混入)')
    plt.title("1. 元の波形（重なっていて、目視では違いが分かりにくい）")
    plt.ylabel("振幅")
    plt.legend(loc='upper right')
    plt.grid(True, alpha=0.3)
    
    # --- 2段目：FFTスペクトル比較 ---
    plt.subplot(2, 1, 2)
    
    # 正常データのスペクトル
    plt.plot(n_freq, n_amp, color='green', alpha=0.5, label=f'正常スペクトル (Score: {n_score:.3f})')
    
    # 異常データのスペクトル
    plt.plot(a_freq, a_amp, color='red', alpha=0.8, label=f'異常スペクトル (Score: {a_score:.3f})')
    
    # 監視バンドのエリア表示
    plt.axvspan(200, 400, color='yellow', alpha=0.2, label='監視バンド (200-400Hz)')
    
    # 閾値ライン
    plt.axhline(threshold, color='blue', linestyle='--', label=f'判定閾値 ({threshold:.3f})')
    
    plt.title("2. FFTスペクトル比較（300Hzに明確な『異常ピーク』が出現）")
    plt.xlabel("周波数 [Hz]")
    plt.ylabel("振幅")
    plt.xlim(0, 500) # 見やすいように500Hzまで拡大
    plt.legend(loc='upper right')
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    FS = 1000
    
    # 1. データの生成
    # 正常データ（学習用：基準を作るためのデータ）
    t, normal_data = generate_data(duration=2.0, fs=FS, has_anomaly=False)
    # 異常データ（テスト用：300Hzの異常を混ぜたもの）
    _, anomaly_data = generate_data(duration=2.0, fs=FS, has_anomaly=True)
    
    # 2. 基準値（閾値）の決定
    # まず、正常データだけでスコアを出し「普段の状態」を把握する
    normal_res = fft_anomaly_score(normal_data, FS, target_band=(200, 400))
    
    # 「正常データのスコア」に少し余裕（マージン）を持たせた値を閾値にする
    # （本来は多数の正常データから統計的に決めますが、今回は簡易的に設定）
    normal_score = normal_res[2]
    threshold = normal_score + 0.05
    
    # 3. 異常検知の実行
    # 決定した閾値を使って、異常データ（未知のデータ）の判定を行う
    anomaly_res = fft_anomaly_score(anomaly_data, FS, target_band=(200, 400))   
    
    print(f"正常データのスコア: {normal_res[2]:.4f}")
    print(f"異常データのスコア: {anomaly_res[2]:.4f}")
    print(f"判定閾値: {threshold:.4f}")
    
    if anomaly_res[2] > threshold:
        print(">> 判定結果: 【異常検知】（スペクトルに未知のピークがあります）")
    else:
        print(">> 判定結果: 正常")

    # 4. 描画
    plot_comparison(t, normal_data, anomaly_data, normal_res, anomaly_res, threshold)

import numpy as np
import matplotlib.pyplot as plt

# 日本語フォント設定
plt.rcParams['font.family'] = 'Meiryo'

def plot_spectrum(data, fs, title="FFTスペクトル解析"):
    """
    FFTを行い、周波数スペクトル（成分表）を描画する関数
    
    Parameters
    ----------
    data : array_like
        解析したい1次元の波形データ。
    fs : int
        サンプリング周波数 (Hz)。
    """
    n = len(data)
    
    # 1. FFT実行（高速フーリエ変換）
    # 時間軸のデータを周波数軸のデータに変換します
    F = np.fft.fft(data)
    
    # 2. 横軸（周波数）を作成
    freq = np.fft.fftfreq(n, d=1/fs)
    
    # 3. 縦軸（振幅）を計算
    # np.absで絶対値（大きさ）にし、データ点数 N/2 で割って正規化します
    # (正規化しないと、データ数が増えるほど値が巨大になってしまうため)
    amplitude = np.abs(F) / (n / 2)
    
    # 4. 表示用に「右半分（正の周波数）」だけを取り出す
    # FFTの結果は左右対称になるため、0Hz〜fs/2Hz の範囲だけ見ればOKです
    mask = freq > 0
    freq_pos = freq[mask]
    amp_pos = amplitude[mask]
    
    # 5. 描画
    plt.figure(figsize=(10, 6))
    
    plt.plot(freq_pos, amp_pos, color='blue')
    
    plt.title(title)
    plt.xlabel('周波数 [Hz]')
    plt.ylabel('振幅 (信号の強さ)')
    plt.grid(True, alpha=0.3)
    
    # 見やすいように範囲を制限（今回は500Hzまで）
    plt.xlim(0, 500)
    
    plt.tight_layout()
    plt.show()

# ==========================================
# テストデータの生成と実行
# ==========================================

def generate_composite_data(duration=1.0, fs=1000):
    """
    テストデータ生成：
    耳で聞くと「ポー」という和音に聞こえるが、波形は複雑なデータ。
    1. 50Hz の低い音 (振幅 1.0)
    2. 120Hz の高い音 (振幅 0.5)
    3. ランダムノイズ
    """
    t = np.linspace(0, duration, int(fs * duration))
    
    # 信号1: 50Hz (強い)
    sig1 = 1.0 * np.sin(2 * np.pi * 50 * t)
    
    # 信号2: 120Hz (中くらい)
    sig2 = 0.5 * np.sin(2 * np.pi * 120 * t)
    
    # ノイズ (弱い)
    np.random.seed(42)
    noise = np.random.normal(0, 0.2, len(t))
    
    return t, sig1 + sig2 + noise

if __name__ == "__main__":
    # パラメータ設定
    FS = 1000
    
    # 1. データ生成
    t, data = generate_composite_data(duration=1.0, fs=FS)
    
    # 2. 元波形の確認
    plt.figure(figsize=(10, 3))
    plt.plot(t[:500], data[:500], color='gray', linewidth=0.8) # 見やすくするため最初の0.5秒だけ表示
    plt.title("元の波形（波が重なり合って、何Hzか分からない）")
    plt.xlabel("時間 [s]")
    plt.tight_layout()
    plt.show()
    
    # 3. FFTスペクトル解析の実行
    plot_spectrum(data, FS, title="【解析結果】FFTスペクトル")

import numpy as np
import matplotlib.pyplot as plt

# 日本語フォント設定
plt.rcParams['font.family'] = 'Meiryo'

def generate_composite_data(duration=1.0, fs=1000):
    """
    テストデータ生成関数
    AIに学習させたい「特徴」を含んだ波形を作ります。
    ここでは「正常なモーター音（50Hz）」に「異常な金属音（300Hz）」が混ざった状態を模倣します。
    """
    t = np.linspace(0, duration, int(fs * duration))
    
    # 信号1: ベースの回転音 (50Hz)
    sig1 = 1.0 * np.sin(2 * np.pi * 50 * t)
    
    # 信号2: 異常音 (300Hz)
    # これの有無や強さが、AIにとっての判断材料になります
    sig2 = 0.5 * np.sin(2 * np.pi * 300 * t)
    
    # ノイズ
    np.random.seed(42)
    noise = np.random.normal(0, 0.2, len(t))
    
    return t, sig1 + sig2 + noise

def extract_fft_features(data, fs):
    """
    音響データをAI（LightGBMなど）学習用に前処理する関数
    FFTを行い、スペクトルの形状から「数値特徴量」を算出します。
    
    Parameters
    ----------
    data : ndarray
        入力波形データ（1次元）。
    fs : int
        サンプリング周波数。
        
    Returns
    -------
    features : dict
        抽出された特徴量を格納した辞書。
        これをpd.DataFrameなどに変換してAIに入力します。
    """
    n = len(data)
    
    # 1. FFT実行
    F = np.fft.fft(data)
    freq = np.fft.fftfreq(n, d=1/fs)
    amp = np.abs(F) / (n / 2) # 振幅スペクトル
    
    # 正の領域だけ抽出
    mask = freq > 0
    freq = freq[mask]
    amp = amp[mask]
    
    # 2. 特徴量抽出 (Feature Engineering)
    # ここで「波形の性格」を表す数値を計算します
    
    # (A) ピーク周波数（一番強い成分は何Hzか？）
    peak_idx = np.argmax(amp)
    peak_freq = freq[peak_idx]
    peak_amp = amp[peak_idx]
    
    # (B) スペクトル重心（音の高さの平均的な重心）
    # 重心 = Σ(周波数 * 振幅) / Σ(振幅)
    spectral_centroid = np.sum(freq * amp) / np.sum(amp)
    
    # (C) 全エネルギー（音の総量）
    total_energy = np.sum(amp)
    
    # (D) 特定帯域の強度（例: 200Hz以上の高周波成分の割合）
    high_freq_mask = freq > 200
    high_freq_energy = np.sum(amp[high_freq_mask])
    
    # 辞書形式でまとめる
    features = {
        "Peak_Freq_Hz": peak_freq,       # ピーク周波数
        "Peak_Amp": peak_amp,            # ピーク強度
        "Centroid_Hz": spectral_centroid,# 重心周波数
        "Total_Energy": total_energy,    # 全エネルギー
        "High_Freq_Energy": high_freq_energy # 高周波成分量
    }
    
    return features

# ==========================================
# 実行サンプル
# ==========================================
if __name__ == "__main__":
    # パラメータ設定
    FS = 1000
    DURATION = 1.0
    
    # 1. テストデータ生成
    t, wave_data = generate_composite_data(duration=DURATION, fs=FS)
    
    # 2. AI用特徴量抽出の実行
    # 波形データ(数千個の数字)が、わずか5個の「意味のある数字」に圧縮されます
    features = extract_fft_features(wave_data, FS)
    
    # 3. 結果確認（これがAIへの入力データになります）
    print("-" * 40)
    print("AI入力用データ (特徴量ベクトル)")
    print("-" * 40)
    for key, val in features.items():
        print(f"{key}: {val:.4f}")
    print("-" * 40)
    print("※ これを CSVの1行として保存し、LightGBMなどで学習させます。")