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

【第1回】
概要と簡単な実装編

本記事では、FFTの概要の説明と、pythonによってテストデータに対して実際にFFT解析を実行します。FFTの基礎を具体例を交えて視覚的に理解できます。

目次

FFT解析の概要

FFT解析とは?

FFT解析という言葉は、エンジニアならだれもが一度は聞いたことはあると思います。そもそもFFTとはFast Fourier Transformの略であり、日本語で言うと、高速フーリエ変換となります。日本語からもわかるように、フーリエ変換を高速で行うものです。では、フーリエ変換とは何か?というところですが、これに関しては基礎の基礎なので大学の教科書に任せるとして、本記事では触れません。「フーリエ変換は時間領域のデータを周波数領域に変換するもの」ぐらいのざっくりした説明でとどめておきたいと思います。

FFTの基礎となる理論

時間領域のデータを周波数領域に変換するための最も基礎となる理論が、(連続時間で定義された)フーリエ変換です。しかし、現実世界での実際のデジタル信号処理では、連続的な信号ではなく、離散的にサンプリングされたデータを扱います。ここで登場するのが離散フーリエ変換、もとい、DFT(Discrete Fourier Transform)**です。
FFTとは、このDFTの計算をアルゴリズムを工夫して計算量を減らして速くしたものであり、数学的にDFTと同じ結果が高速に得られます。

  • 連続フーリエ変換
    • 連続時間で定義された最も基礎となる理論
  • 離散フーリエ変換(DFT):
    • 離散時間で定義されたフーリエ変換
  • 高速フーリエ変換(FFT):
    • DFTと数学的には同じだが、計算量を減らして高速にしたもの

 

もう少し詳細を詳しく説明します。
まず、連続時間におけるフーリエ変換は以下の式で定義されます。

\[X(\omega) = \int_{-\infty}^{\infty} x(t) \cdot e^{-j\omega t} dt\]

ここで:

\[
\begin{array}{ll}
x(t) & :時間領域の連続信号\\
X(\omega) & :周波数領域の連続信号 \\
j & :虚数単位
\end{array}
\]

この変換により、時間軸上の連続信号を周波数軸上の連続スペクトルに変換できます。
このフーリエ変換を離散化した、DFTは以下の数式で表されます:

\[X[k] = \sum_{n=0}^{N-1} x[n] \cdot \exp\left(-j\frac{2\pi kn}{N}\right)\]

ここで:

\[
\begin{array}{ll}
x[n] & :時間領域の離散信号(n = 0, 1, \ldots, N-1) \\
X[k] & :周波数領域の離散信号(k = 0, 1, \ldots, N-1) \\
N & :データ点数 \\
j & :虚数単位
\end{array}
\]

この変換により、時間軸上のN個のデータ点を、周波数軸上のN個の成分に変換できます。

実例を使ってFFT解析を理解する

内容

見出し2-2

見出し2-3

内容

Python実装を確認する

メインプログラム(main.py)

aaa

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)

関数ファイル(FFT_lib.py)

ああああ

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  : 
        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()

次回予告

次回は、Pythonコードの内容を細かく説明します。

よかったらシェアしてね!
  • URLをコピーしました!
  • URLをコピーしました!

コメント

コメント一覧 (1件)

A WordPress Commenter へ返信する コメントをキャンセル

CAPTCHA


目次