はじめに
フーリエ変換は、信号処理、画像解析、物理、生命科学、音声解析、センサデータ解析など、さまざまな場面で出てきます。
ただ、初めて学ぶと、
フーリエ変換は、時間変化するデータを周波数成分に分けるものです
と言われても、少し分かったような、分からないような気分になります。
この記事では、難しい証明ではなく、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
のように変えて、再構成される波形がどう変わるかを見てみてください。