はじめに

フーリエ変換は、信号処理、画像解析、物理、生命科学、音声解析、センサデータ解析など、さまざまな場面で出てきます。

ただ、初めて学ぶと、

フーリエ変換は、時間変化するデータを周波数成分に分けるものです

と言われても、少し分かったような、分からないような気分になります。

この記事では、難しい証明ではなく、Pythonでダミーデータを作りながら、

フーリエ変換は、混ざった波の中から「どの周波数がどれくらい含まれているか」を見ている

という感覚をつかむことを目指します。

Google Colabでそのまま動くコードを使います。

この記事でやること

この記事では、次の流れでフーリエ変換を見ていきます。

2 Hz、8 Hz、20 Hz の3つの波を作る
それらを足し合わせて「混ざった波」を作る
時間波形だけを見ると、何が混ざっているか分かりにくいことを確認する
FFTで周波数ごとの成分に分ける
ノイズが入っても、主要な周波数成分が見えることを確認する
おまけとして、FFTを使った簡単な低周波成分の取り出しを見る

まず用語を整理する

波とは何か

ここでは、時間とともに周期的に上下するデータを「波」と呼ぶことにします。

たとえば、単純な正弦波は次のように書けます。

$$
x(t) = A \sin(2\pi f t)
$$

ここで、

$x(t)$ は時刻 $t$ における信号の値
$A$ は振幅
$f$ は周波数
$t$ は時間

です。

振幅とは何か

振幅は、波の大きさです。

振幅が大きいほど、波の上下の変化が大きくなります。

周波数とは何か

周波数は、1秒間に何回繰り返すかを表します。

単位は Hz です。

たとえば、

2 Hz は、1秒間に2回繰り返す
8 Hz は、1秒間に8回繰り返す
20 Hz は、1秒間に20回繰り返す

という意味です。

時間領域とは何か

横軸を時間にして、信号の変化を見る見方を「時間領域」と呼びます。

普段よく見る波形グラフは、だいたい時間領域のグラフです。

周波数領域とは何か

横軸を周波数にして、どの周波数がどれくらい含まれているかを見る見方を「周波数領域」と呼びます。

フーリエ変換は、時間領域のデータを周波数領域へ変換する道具です。

スペクトルとは何か

周波数ごとの成分の大きさを並べたものを、ここでは「スペクトル」と呼びます。

ざっくり言えば、

スペクトルを見ると、信号の中にどんな周期の波が入っているかが分かる

ということです。

フーリエ変換の直感

今回作る信号は、次のように3つの波を足し合わせたものです。

$$
x(t)=
1.0 \sin(2\pi \cdot 2t)
+
0.6 \sin(2\pi \cdot 8t)
+
0.3 \sin(2\pi \cdot 20t)
$$

つまり、この信号には、

2 Hz の波が振幅 1.0
8 Hz の波が振幅 0.6
20 Hz の波が振幅 0.3

だけ入っています。

しかし、これらを足し合わせると、時間波形だけを見ても「2 Hz、8 Hz、20 Hzが入っている」とは一目では分かりにくくなります。

フーリエ変換を使うと、この混ざった波から、

2 Hz が強い
8 Hz も入っている
20 Hz も少し入っている

という情報を取り出せます。

離散フーリエ変換の式

コンピュータで扱うデータは、連続的な関数ではなく、サンプリングされた数列です。

そのため、実際には離散フーリエ変換、Discrete Fourier Transform、略してDFTを使います。

DFTは次のように書けます。

$$
X_k =
\sum_{n=0}^{N-1}
x_n
\exp\left(
-2\pi i \frac{kn}{N}
\right)
$$

ここで、

$x_n$ は時系列データ
$N$ はデータ点数
$X_k$ は周波数成分
$i$ は虚数単位

です。

この式だけ見ると少し怖いですが、直感的には、

いろいろな周波数の波をものさしとして用意し、元の信号とどれくらい合うかを調べている

と考えると分かりやすいです。

この記事では、DFTを高速に計算するFFT、Fast Fourier Transformを使います。

Pythonでは np.fft.rfft() を使うと、実数の時系列データに対するFFTを簡単に計算できます。

Google Colabで実行するコード

# ============================================
# Fourier transform minimal demo
# mixed waves -> frequency components
# ============================================

import os
import numpy as np
import matplotlib.pyplot as plt

# 図を保存するフォルダ
os.makedirs("fig", exist_ok=True)

# 再現性のために乱数を固定
rng = np.random.default_rng(42)

# ============================================
# 1. ダミーデータを作る
# ============================================

fs = 200          # サンプリング周波数 [Hz]
duration = 2.0    # 測定時間 [s]

t = np.arange(0, duration, 1/fs)
N = len(t)

print(f"Number of data points: {N}")
print(f"Sampling frequency: {fs} Hz")
print(f"Duration: {duration} s")
print(f"Frequency resolution: {fs/N:.2f} Hz")

# 3つの波を作る
f1, A1 = 2, 1.0    # 2 Hz, 振幅 1.0
f2, A2 = 8, 0.6    # 8 Hz, 振幅 0.6
f3, A3 = 20, 0.3   # 20 Hz, 振幅 0.3

wave1 = A1 * np.sin(2*np.pi*f1*t)
wave2 = A2 * np.sin(2*np.pi*f2*t)
wave3 = A3 * np.sin(2*np.pi*f3*t)

# 混ざった波
x = wave1 + wave2 + wave3

# ============================================
# 2. FFT用の関数を定義する
# ============================================

def one_sided_amplitude_spectrum(signal, fs):
    """
    実数の時系列信号から、片側振幅スペクトルを計算する関数。

    signal : 1次元の時系列データ
    fs     : サンプリング周波数 [Hz]
    """
    N = len(signal)

    # 実数信号用のFFT
    fft_values = np.fft.rfft(signal)

    # 各FFT成分に対応する周波数
    freq = np.fft.rfftfreq(N, d=1/fs)

    # 振幅スペクトルに変換
    # 片側スペクトルなので、DC成分とナイキスト周波数以外を2倍する
    amp = np.abs(fft_values) / N * 2
    amp[0] /= 2
    if N % 2 == 0:
        amp[-1] /= 2

    return freq, amp


def print_top_frequencies(freq, amp, title, n=6):
    """
    振幅が大きい周波数成分を表示する関数。
    0 Hz成分は除いて表示する。
    """
    top_idx = np.argsort(amp[1:])[-n:][::-1] + 1

    print()
    print(title)
    print("-" * len(title))

    for i in top_idx:
        print(f"{freq[i]:6.1f} Hz : amplitude = {amp[i]:.3f}")


def save_and_show(filename):
    """
    図を保存して表示する関数。
    """
    plt.tight_layout()
    plt.savefig(f"fig/{filename}", dpi=160, bbox_inches="tight")
    plt.show()


# ============================================
# 3. 成分の波を描く
# ============================================

plt.figure(figsize=(10, 4))
plt.plot(t, wave1, label="2 Hz")
plt.plot(t, wave2, label="8 Hz")
plt.plot(t, wave3, label="20 Hz")

plt.title("Three sine waves used as components")
plt.xlabel("time [s]")
plt.ylabel("amplitude")
plt.legend()
plt.grid(True, alpha=0.3)

save_and_show("fig1_components.png")


# ============================================
# 4. 混ざった波を描く
# ============================================

plt.figure(figsize=(10, 4))
plt.plot(t, x)

plt.title("Mixed signal in the time domain")
plt.xlabel("time [s]")
plt.ylabel("amplitude")
plt.grid(True, alpha=0.3)

save_and_show("fig2_mixed_signal.png")


# ============================================
# 5. 混ざった波をFFTする
# ============================================

freq, amp = one_sided_amplitude_spectrum(x, fs)

print_top_frequencies(
    freq,
    amp,
    title="Clean signal: large frequency components",
    n=3
)

plt.figure(figsize=(10, 4))
plt.plot(freq, amp)
plt.scatter([f1, f2, f3], [A1, A2, A3], zorder=3, label="true components")

plt.xlim(0, 50)
plt.title("Amplitude spectrum of the mixed signal")
plt.xlabel("frequency [Hz]")
plt.ylabel("amplitude")
plt.legend()
plt.grid(True, alpha=0.3)

save_and_show("fig3_spectrum_clean.png")


# ============================================
# 6. ノイズを加えてみる
# ============================================

noise = rng.normal(loc=0, scale=0.4, size=N)
x_noisy = x + noise

plt.figure(figsize=(10, 4))
plt.plot(t, x_noisy, label="noisy signal")
plt.plot(t, x, linewidth=2, alpha=0.8, label="original mixed signal")

plt.title("Mixed signal with noise")
plt.xlabel("time [s]")
plt.ylabel("amplitude")
plt.legend()
plt.grid(True, alpha=0.3)

save_and_show("fig4_noisy_signal.png")


# ============================================
# 7. ノイズ入り信号をFFTする
# ============================================

freq_noisy, amp_noisy = one_sided_amplitude_spectrum(x_noisy, fs)

print_top_frequencies(
    freq_noisy,
    amp_noisy,
    title="Noisy signal: large frequency components",
    n=6
)

plt.figure(figsize=(10, 4))
plt.plot(freq_noisy, amp_noisy)

plt.xlim(0, 50)
plt.title("Amplitude spectrum of the noisy signal")
plt.xlabel("frequency [Hz]")
plt.ylabel("amplitude")
plt.grid(True, alpha=0.3)

save_and_show("fig5_spectrum_noisy.png")


# ============================================
# 8. おまけ：FFTで低周波成分だけ取り出す
# ============================================

# FFTする
fft_noisy = np.fft.rfft(x_noisy)
freq_for_filter = np.fft.rfftfreq(N, d=1/fs)

# 10 Hzより高い成分を0にする
fft_low = fft_noisy.copy()
fft_low[freq_for_filter > 10] = 0

# 逆FFTで時間波形に戻す
x_low = np.fft.irfft(fft_low, n=N)

plt.figure(figsize=(10, 4))
plt.plot(t, x_noisy, alpha=0.5, label="noisy signal")
plt.plot(t, x_low, linewidth=2, label="FFT low-pass reconstruction <= 10 Hz")

plt.title("Simple FFT-based low-pass reconstruction")
plt.xlabel("time [s]")
plt.ylabel("amplitude")
plt.legend()
plt.grid(True, alpha=0.3)

save_and_show("fig6_lowpass_reconstruction.png")

実行結果の例

上のコードを実行すると、次のような出力が得られます。

Number of data points: 400
Sampling frequency: 200 Hz
Duration: 2.0 s
Frequency resolution: 0.50 Hz

Clean signal: large frequency components
----------------------------------------
   2.0 Hz : amplitude = 1.000
   8.0 Hz : amplitude = 0.600
  20.0 Hz : amplitude = 0.300

Noisy signal: large frequency components
----------------------------------------
   2.0 Hz : amplitude = 1.000
   8.0 Hz : amplitude = 0.579
  20.0 Hz : amplitude = 0.287
  88.0 Hz : amplitude = 0.082
  58.0 Hz : amplitude = 0.074
  13.0 Hz : amplitude = 0.073

きれいな信号では、作った通りに 2 Hz、8 Hz、20 Hz の成分が検出されています。

ノイズを足した信号でも、2 Hz、8 Hz、20 Hz の成分が大きく出ています。

図1：3つの成分波

図1では、3つの正弦波を表示しています。

2 Hz の波はゆっくり上下します。
8 Hz の波はそれより速く上下します。
20 Hz の波はさらに細かく上下します。

周波数が高いほど、同じ時間の中で何度も振動していることが分かります。

図2：混ざった波

図2は、3つの波を足し合わせたものです。

$$
x(t)=
1.0 \sin(2\pi \cdot 2t)
+
0.6 \sin(2\pi \cdot 8t)
+
0.3 \sin(2\pi \cdot 20t)
$$

時間波形を見ると、何か周期的な構造があることは分かります。

しかし、これを見ただけで、

2 Hz
8 Hz
20 Hz

が入っていると読み取るのは、少し大変です。

ここでフーリエ変換を使います。

図3：混ざった波の振幅スペクトル

図3では、横軸が周波数、縦軸が振幅です。

この図を見ると、2 Hz、8 Hz、20 Hz に鋭いピークがあります。

これは、

この信号には、2 Hz、8 Hz、20 Hz の成分が含まれている

ということを意味します。

さらに、ピークの高さは振幅に対応しています。

今回のデータは、

2 Hz の振幅が 1.0
8 Hz の振幅が 0.6
20 Hz の振幅が 0.3

になるように作りました。

スペクトルでも、おおむねその通りの高さになっています。

つまり、時間波形では分かりにくかった混ざり方が、周波数領域ではかなり分かりやすく見えます。

フーリエ変換は何を見ているのか

ここまでを見ると、フーリエ変換は次のように考えられます。

信号の中に、どの周波数の波がどれくらい含まれているかを調べる方法

時間領域では、信号は「時間に沿った値の変化」として見えます。

一方、周波数領域では、信号は「周波数ごとの成分の集まり」として見えます。

同じデータでも、見る方向を変えると、見え方が大きく変わります。

図4：ノイズを加えた波

現実のデータには、たいていノイズが含まれます。

図4では、混ざった波にランダムノイズを加えています。

時間波形を見ると、元のなめらかな波よりもかなりギザギザしています。

この状態で、もとの周波数成分を見つけられるでしょうか。

図5：ノイズ入り信号の振幅スペクトル

図5を見ると、ノイズによって小さなピークがたくさん出ています。

しかし、それでも 2 Hz、8 Hz、20 Hz の成分は大きく残っています。

このように、時間波形ではノイズで見えにくくなっていても、周波数領域では主要な周期成分が見えることがあります。

もちろん、ノイズが非常に大きい場合や、信号が時間とともに変化する場合には、FFTだけで簡単に解釈できないこともあります。

それでも、

時間波形では見えにくい構造を、周波数の視点から見直す

という意味で、フーリエ変換はとても強力です。

図6：おまけ：低周波成分だけ取り出す

最後に、おまけとして、FFTを使って 10 Hz 以下の成分だけを残してみました。

今回の信号には、

2 Hz
8 Hz
20 Hz

の成分が含まれています。

10 Hz以下だけを残すと、2 Hzと8 Hzは残りますが、20 Hzの成分は除かれます。

そのため、図6の再構成波形は、ノイズ入り信号よりなめらかになります。

ただし、ここで大事なのは、

なめらかになったからといって、必ず真の信号に近づいたとは限らない

という点です。

今回の場合、20 Hz成分も人工的に入れた本物の信号です。
10 Hzより高い成分を消すという操作は、ノイズだけでなく、20 Hzの信号成分も消しています。

フィルタ処理では、

何を信号とみなし、何をノイズとみなすか

を考えることが重要です。

サンプリング周波数と周波数分解能

FFTを使うときには、サンプリング周波数と測定時間も重要です。

今回のコードでは、

fs = 200
duration = 2.0

としました。

つまり、1秒あたり200点で、2秒間のデータを作っています。

データ点数は、

$$
N = 200 \times 2 = 400
$$

です。

周波数分解能は、おおまかに次のように表せます。

$$
\Delta f = \frac{f_s}{N}
$$

今回の場合は、

$$
\Delta f = \frac{200}{400} = 0.5 \ \mathrm{Hz}
$$

です。

これは、FFTの周波数軸が 0.5 Hz 刻みになる、という意味です。

また、サンプリング周波数が $f_s$ のとき、原理的に表現できる最大周波数はおおむね、

$$
\frac{f_s}{2}
$$

です。

これをナイキスト周波数と呼びます。

今回のサンプリング周波数は 200 Hz なので、ナイキスト周波数は 100 Hz です。

FFTの結果を読むときの注意

FFTは便利ですが、結果を読むときにはいくつか注意点があります。

1. 振幅スペクトルだけでは位相情報を見ていない

FFTの結果は、本当は複素数です。

この記事では分かりやすさのために、np.abs() を使って大きさだけを見ています。

つまり、振幅スペクトルでは「どの周波数がどれくらい含まれているか」は見ていますが、波のずれ、つまり位相の情報は直接見ていません。

2. 測定時間が短いと周波数分解能が粗くなる

測定時間が短いと、近い周波数を分けにくくなります。

今回の測定時間は2秒なので、周波数分解能は0.5 Hzです。

より細かく周波数を見たい場合は、より長い時間のデータが必要になります。

3. 現実のデータではピークがきれいに出るとは限らない

今回の例では、2 Hz、8 Hz、20 Hzの波がきれいに入るように作っています。

現実のデータでは、

ノイズ
ドリフト
測定時間の短さ
周波数が時間とともに変わる
窓関数の影響
サンプリング周波数の不足

などにより、スペクトルがもっと複雑になります。

4. 周波数成分があることと、原因が分かることは別

FFTで特定の周波数ピークが見つかったとしても、それだけで原因が分かるわけではありません。

FFTはあくまで、

その周波数成分がデータに含まれている

ことを示します。

その周波数が何に由来するのかは、実験条件や測定対象の知識と合わせて考える必要があります。

まとめ

この記事では、ダミーデータを使ってフーリエ変換の基本的な考え方を見ました。

重要なポイントは次の通りです。

時間領域では、信号を時間変化として見る
周波数領域では、信号を周波数成分の集まりとして見る
フーリエ変換は、混ざった波を周波数ごとの成分に分ける道具である
FFTを使うと、コンピュータで高速にフーリエ変換を計算できる
ノイズがあっても、主要な周期成分がスペクトル上で見えることがある
ただし、FFTの結果は測定条件や前処理に強く影響される

フーリエ変換は、数式だけを見ると少し難しく見えます。

しかし、実際に波を作り、混ぜて、FFTして、スペクトルを見ると、

混ざった波を、周波数ごとに分けて見る方法

という感覚がつかみやすくなります。

おまけ：試してみると理解が深まる変更

余裕があれば、コードの値を変えて遊んでみると理解が深まります。

たとえば、次のような変更がおすすめです。

1. 周波数を変えてみる

f1, A1 = 3, 1.0
f2, A2 = 12, 0.6
f3, A3 = 30, 0.3

ピークの位置が変わることを確認できます。

2. 振幅を変えてみる

f1, A1 = 2, 0.2
f2, A2 = 8, 1.0
f3, A3 = 20, 0.5

ピークの高さが変わることを確認できます。

3. ノイズを大きくしてみる

noise = rng.normal(loc=0, scale=1.0, size=N)

ノイズが大きくなると、スペクトルの小さなピークが増え、主要な成分が見えにくくなることがあります。

4. 測定時間を長くしてみる

duration = 5.0

測定時間を長くすると、周波数分解能が細かくなります。

5. 低周波成分だけを残すしきい値を変えてみる

fft_low[freq_for_filter > 5] = 0

または、

fft_low[freq_for_filter > 30] = 0

のように変えて、再構成される波形がどう変わるかを見てみてください。