0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

はじめに

時系列データを扱うとき、よく出てくる言葉にサンプリング周波数があります。

例えば、センサー、音声、心電図、顕微鏡画像の時系列、物理計測データなどでは、連続的に変化している現象を、一定時間ごとに測定してデータ化します。

このとき、

どれくらい細かい時間間隔で測るか

がとても重要です。

サンプリングが粗すぎると、本当は速く振動している信号が、まるでゆっくり振動しているように見えてしまうことがあります。

この現象をエイリアシング、または折り返しと呼びます。

この記事では、Pythonでダミーデータを作りながら、

  • サンプリングとは何か
  • サンプリング周波数とは何か
  • ナイキスト周波数とは何か
  • エイリアシングとは何か
  • FFTで見ても、サンプリング後には元の高周波を区別できないこと

を、図で確認します。


この記事のゴール

この記事の結論を先に書くと、次のようになります。

サンプリング周波数が足りないと、高い周波数の信号が低い周波数の信号に見えてしまう。

より具体的には、この記事では

本当は7 Hzの正弦波なのに、10 Hzでサンプリングすると3 Hzの波に見えてしまう

という例を扱います。


用語の整理

まず、この記事で使う用語を整理します。

信号

時間とともに変化する量を、ここでは信号と呼びます。

今回は、次のような単純な正弦波を考えます。

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

ここで、$t$ は時間、$f$ は信号の周波数です。

例えば $f=7$ Hz なら、1秒間に7回振動する波です。


サンプリング

連続的な信号を、一定の時間間隔で測ることをサンプリングと呼びます。

例えば、0.1秒ごとに値を測るなら、

$$
t = 0, 0.1, 0.2, 0.3, \ldots
$$

のような時刻で信号を記録することになります。


サンプリング周波数

1秒間に何回サンプリングするかを、サンプリング周波数と呼びます。

記号としては $f_s$ をよく使います。

サンプリング間隔を $\Delta t$ とすると、

$$
f_s = \frac{1}{\Delta t}
$$

です。

例えば、0.1秒ごとに測るなら、

$$
f_s = \frac{1}{0.1} = 10
$$

なので、サンプリング周波数は10 Hzです。


ナイキスト周波数

サンプリング周波数 $f_s$ の半分を、ナイキスト周波数と呼びます。

$$
f_N = \frac{f_s}{2}
$$

例えば、サンプリング周波数が10 Hzなら、

$$
f_N = \frac{10}{2} = 5
$$

なので、ナイキスト周波数は5 Hzです。

重要なのは、ざっくり言うと、

サンプリング周波数 $f_s$ で測ったデータから、$f_s/2$ より高い周波数を正しく区別することはできない

ということです。

この $f_s/2$ が、サンプリング後に扱える周波数の上限になります。


Google Colabで実行する準備

以下のコードはGoogle Colabでそのまま実行できます。

実行すると、figures フォルダにPNG画像が保存されます。

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

# 図の保存先
os.makedirs("figures", exist_ok=True)

1. まずは7 Hzの信号を作る

今回は、1秒間に7回振動する信号を作ります。

$$
x(t) = \sin(2\pi \cdot 7 t)
$$

です。

ここでは、かなり細かい時間間隔で作った波形を「本当の連続波形に近いもの」として扱います。

厳密にはコンピュータ上では完全な連続時間信号は作れませんが、十分細かく作れば、可視化用の基準として使えます。


2. 十分細かくサンプリングした場合と、粗くサンプリングした場合を比べる

7 Hzの信号を、

  • 50 Hzでサンプリングする場合
  • 10 Hzでサンプリングする場合

で比べます。

50 Hzなら、1秒間に50点を取ります。
10 Hzなら、1秒間に10点しか取りません。

# 真の信号として扱う高解像度の波形
f_signal = 7.0     # 信号の周波数 [Hz]
duration = 1.0     # 表示する時間 [s]
fs_ref = 2000      # 可視化用の細かいサンプリング周波数 [Hz]

t_ref = np.arange(0, duration, 1/fs_ref)
x_ref = np.sin(2*np.pi*f_signal*t_ref)

# 2種類のサンプリング周波数
fs_good = 50
fs_bad = 10

t_good = np.arange(0, duration, 1/fs_good)
x_good = np.sin(2*np.pi*f_signal*t_good)

t_bad = np.arange(0, duration, 1/fs_bad)
x_bad = np.sin(2*np.pi*f_signal*t_bad)

# 図を作成
fig, axes = plt.subplots(2, 1, figsize=(10, 7), sharex=True)

axes[0].plot(t_ref, x_ref, label="true signal: 7 Hz")
axes[0].scatter(t_good, x_good, s=35, label="sampled at 50 Hz")
axes[0].set_title("Enough sampling: 7 Hz signal sampled at 50 Hz")
axes[0].set_ylabel("Amplitude")
axes[0].legend()
axes[0].grid(alpha=0.3)

axes[1].plot(t_ref, x_ref, label="true signal: 7 Hz")
axes[1].scatter(t_bad, x_bad, s=50, label="sampled at 10 Hz")
axes[1].plot(t_bad, x_bad, marker="o", linestyle="--", label="line connecting samples")
axes[1].set_title("Sparse sampling: 7 Hz signal sampled at 10 Hz")
axes[1].set_xlabel("Time [s]")
axes[1].set_ylabel("Amplitude")
axes[1].legend()
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.savefig("figures/fig1_sampling_comparison.png", dpi=180)
plt.show()

図1:十分なサンプリングと粗いサンプリングの比較

fig1_sampling_comparison.png

上の図では、50 Hzでサンプリングした場合、点が十分に細かく並んでおり、7 Hzの波形をそれなりに追えています。

一方、10 Hzでサンプリングした場合は、点の数が少なすぎます。

点だけを見たり、点を線でつないだりすると、本当の7 Hzの波ではなく、もっとゆっくりした波に見えてしまいます。

これがエイリアシングの入口です。


3. 10 Hzサンプリングでは7 Hzが3 Hzに見える

サンプリング周波数が10 Hzのとき、ナイキスト周波数は

$$
f_N = \frac{10}{2} = 5 \ \mathrm{Hz}
$$

です。

つまり、10 Hzでサンプリングしたデータから正しく扱えるのは、基本的には5 Hz以下の成分です。

今回の信号は7 Hzです。
7 Hzはナイキスト周波数5 Hzを超えています。

そのため、サンプリング後のデータでは、7 Hzの信号が3 Hzの信号のように見えてしまいます。

なぜ3 Hzかというと、10 Hzを基準にして折り返されるからです。

$$
7 \ \mathrm{Hz}
\rightarrow
|10 - 7| = 3 \ \mathrm{Hz}
$$

次の図では、7 Hzの真の波形と、3 Hzに見える波形を重ねて表示します。

fs_bad = 10
f_alias = abs(fs_bad - f_signal)  # 10 - 7 = 3 Hz

t_ref = np.arange(0, duration, 1/fs_ref)
x_true = np.sin(2*np.pi*f_signal*t_ref)

# 7 Hzを10 Hzでサンプリングした場合、3 Hzの波に位相を反転させたものと同じ点を通る
x_alias = -np.sin(2*np.pi*f_alias*t_ref)

t_sample = np.arange(0, duration, 1/fs_bad)
x_sample = np.sin(2*np.pi*f_signal*t_sample)

plt.figure(figsize=(10, 5))
plt.plot(t_ref, x_true, label="true signal: 7 Hz")
plt.plot(t_ref, x_alias, linestyle="--", label="alias candidate: 3 Hz (phase-flipped)")
plt.scatter(t_sample, x_sample, s=70, zorder=3, label="samples at 10 Hz")

plt.title("The same samples can look like a lower-frequency wave")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("figures/fig2_alias_illusion.png", dpi=180)
plt.show()

図2:7 Hzの信号が3 Hzの信号のように見える

fig2_alias_illusion.png

青い線が本当の7 Hzの信号です。
点は、10 Hzでサンプリングした値です。

破線は、3 Hzの信号です。正確には位相が反転した3 Hzの波ですが、サンプリングされた点だけを見ると、7 Hzの信号と区別できません。

つまり、サンプリング後のデータだけを見ると、

これは7 Hzだったのか?
それとも3 Hzだったのか?

を判断できなくなります。

この「高い周波数が低い周波数に見えてしまう」現象が、エイリアシングです。


4. エイリアシングを式で見る

サンプリング周波数を $f_s$、本当の信号周波数を $f$ とします。

サンプリング後に見える周波数は、次のように書けます。

$$
f_{\mathrm{alias}} =\left|\left( f + \frac{f_s}{2} \right) \bmod f_s- \frac{f_s}{2}
\right|
$$

少し難しく見えますが、意味はシンプルです。

サンプリング周波数を超えた成分は、ナイキスト周波数を境に折り返して見える

ということです。

10 Hzでサンプリングする場合、観測できる周波数の上限は5 Hzです。

そのため、

  • 1 Hzは1 Hzに見える
  • 3 Hzは3 Hzに見える
  • 5 Hzは5 Hz付近に見える
  • 7 Hzは3 Hzに見える
  • 9 Hzは1 Hzに見える
  • 11 Hzは1 Hzに見える

というように、高い周波数が低い周波数へ折り返されます。

これを図にしてみます。

fs = 10.0

def alias_frequency(f, fs):
    """
    周波数 f がサンプリング周波数 fs で観測されたとき、
    どの周波数として見えるかを返す。
    """
    return np.abs((f + fs/2) % fs - fs/2)

freqs = np.linspace(0, 30, 1201)
aliases = alias_frequency(freqs, fs)

plt.figure(figsize=(10, 5))
plt.plot(freqs, aliases, label="observed alias frequency")
plt.axvline(fs/2, linestyle="--", label="Nyquist frequency = 5 Hz")
plt.axhline(fs/2, linestyle=":", label="maximum observable frequency")
plt.scatter([7], [alias_frequency(7, fs)], s=80, zorder=3, label="7 Hz -> 3 Hz")

plt.title("Aliasing map when sampling frequency is 10 Hz")
plt.xlabel("True signal frequency [Hz]")
plt.ylabel("Observed frequency after sampling [Hz]")
plt.ylim(-0.1, 5.5)
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("figures/fig3_aliasing_map.png", dpi=180)
plt.show()

図3:10 Hzサンプリングでのエイリアシングマップ

fig3_aliasing_map.png

この図では、横軸が本当の信号周波数、縦軸がサンプリング後に見える周波数です。

サンプリング周波数が10 Hzなので、ナイキスト周波数は5 Hzです。

5 Hzを超えた信号は、そのまま5 Hz以上として見えるのではなく、低い周波数へ折り返されます。

今回の例では、7 Hzが3 Hzに見えています。


5. FFTを使っても、サンプリング後には区別できない

ここで重要なのは、

FFTを使えば元の7 Hzを取り戻せるわけではない

という点です。

サンプリング後のデータだけを見ると、7 Hzの信号と3 Hzの信号は同じように見えます。

そのため、FFTをしても、7 Hzではなく3 Hzのピークとして現れます。

次のコードでは、

  • 本当に3 Hzの信号を10 Hzでサンプリングした場合
  • 本当は7 Hzの信号を10 Hzでサンプリングした場合

を比較します。

fs_fft = 10
duration_fft = 5.0

t_fft = np.arange(0, duration_fft, 1/fs_fft)

def compute_amplitude_spectrum(f_true, fs, t):
    """
    周波数 f_true の正弦波を作り、FFTの振幅スペクトルを返す。
    """
    x = np.sin(2*np.pi*f_true*t)
    x = x - np.mean(x)

    spectrum = np.fft.rfft(x)
    freq_axis = np.fft.rfftfreq(len(x), d=1/fs)

    amp = np.abs(spectrum) / len(x) * 2
    return freq_axis, amp

freq_axis, amp_3hz = compute_amplitude_spectrum(3.0, fs_fft, t_fft)
_, amp_7hz = compute_amplitude_spectrum(7.0, fs_fft, t_fft)

peak_3hz = freq_axis[np.argmax(amp_3hz[1:]) + 1]
peak_7hz = freq_axis[np.argmax(amp_7hz[1:]) + 1]

print(f"True 3 Hz signal -> FFT peak: {peak_3hz:.1f} Hz")
print(f"True 7 Hz signal -> FFT peak: {peak_7hz:.1f} Hz")

plt.figure(figsize=(10, 5))
plt.plot(freq_axis, amp_3hz, marker="o", label="true 3 Hz sampled at 10 Hz")
plt.plot(freq_axis, amp_7hz, marker="s", linestyle="--", label="true 7 Hz sampled at 10 Hz")
plt.axvline(fs_fft/2, linestyle=":", label="Nyquist frequency = 5 Hz")

plt.xlim(0, fs_fft/2)
plt.title("FFT cannot distinguish 3 Hz and aliased 7 Hz after sampling")
plt.xlabel("Frequency [Hz]")
plt.ylabel("Amplitude spectrum")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("figures/fig4_fft_3hz_vs_7hz.png", dpi=180)
plt.show()

図4:3 Hzの信号と7 Hzの信号がFFTでは同じ3 Hzに見える

fig4_fft_3hz_vs_7hz.png

このコードを実行すると、次のような出力になります。

True 3 Hz signal -> FFT peak: 3.0 Hz
True 7 Hz signal -> FFT peak: 3.0 Hz

本当に3 Hzの信号も、7 Hzがエイリアシングした信号も、FFTでは同じ3 Hzのピークになります。

つまり、サンプリング後のデータに対してどれだけ解析をしても、サンプリング時点で失われた情報は基本的には復元できません。

これはとても重要です。


6. サンプリング周波数を変えるとどうなるか

最後に、本当は7 Hzの信号を、いろいろなサンプリング周波数で測ってみます。

そして、FFTで見えるピーク周波数を調べます。

7 Hzの信号を正しく扱うには、少なくとも2倍より大きいサンプリング周波数が必要です。

$$
f_s > 2f
$$

今回なら、

$$
f_s > 14 \ \mathrm{Hz}
$$

が目安です。

ただし、実際の測定では「ぎりぎり2倍」ではなく、もっと余裕を持たせることが多いです。

f_signal = 7.0
duration_est = 5.0

sampling_rates = np.array([8, 10, 12, 13, 15, 20, 50])

rows = []

for fs_i in sampling_rates:
    t_i = np.arange(0, duration_est, 1/fs_i)
    x_i = np.sin(2*np.pi*f_signal*t_i)
    x_i = x_i - np.mean(x_i)

    spectrum = np.fft.rfft(x_i)
    freq_i = np.fft.rfftfreq(len(x_i), d=1/fs_i)
    amp_i = np.abs(spectrum) / len(x_i) * 2

    # DC成分を除いて最大ピークを探す
    idx = np.argmax(amp_i[1:]) + 1
    estimated_peak = freq_i[idx]

    rows.append({
        "sampling frequency [Hz]": fs_i,
        "Nyquist frequency [Hz]": fs_i / 2,
        "FFT peak [Hz]": estimated_peak
    })

df = pd.DataFrame(rows)
display(df)

plt.figure(figsize=(9, 5))
plt.plot(df["sampling frequency [Hz]"], df["FFT peak [Hz]"], marker="o", label="FFT peak after sampling")
plt.axhline(f_signal, linestyle="--", label="true frequency = 7 Hz")
plt.axvline(2*f_signal, linestyle=":", label="2 x true frequency = 14 Hz")

plt.title("Estimated frequency depends on sampling frequency")
plt.xlabel("Sampling frequency [Hz]")
plt.ylabel("Estimated peak frequency [Hz]")
plt.ylim(0, 8)
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("figures/fig5_sampling_rate_vs_estimated_frequency.png", dpi=180)
plt.show()

図5:サンプリング周波数と推定されるピーク周波数

fig5_sampling_rate_vs_estimated_frequency.png

Colabでは、例えば次のような表が表示されます。

sampling frequency [Hz] Nyquist frequency [Hz] FFT peak [Hz]
8 4.0 1.0
10 5.0 3.0
12 6.0 5.0
13 6.5 6.0
15 7.5 7.0
20 10.0 7.0
50 25.0 7.0

サンプリング周波数が14 Hzより低いと、7 Hzの信号を正しく見られません。

15 Hz以上になると、FFTのピークが7 Hzに戻ります。

ただし、14 Hzちょうどのような「ぎりぎり」の条件では、位相やノイズの影響を強く受けることがあります。
実際の測定では、理論上の下限ぴったりを狙うのではなく、十分な余裕を持たせるのが安全です。


7. 実データ解析で気をつけたいこと

今回の例は、ノイズのない単純な正弦波でした。

実際のデータでは、もっと面倒です。

例えば、

  • 複数の周波数成分が混ざっている
  • ノイズが入っている
  • 信号の周波数が時間とともに変化する
  • サンプリング間隔が完全には一定でない
  • 測定前のフィルタ処理が入っている
  • センサーや装置の応答特性がある

といったことがあります。

そのため、実データでは、

FFTをしたらピークが見えたので、それが本当の周波数である

とは、すぐには言えません。

特に、サンプリング周波数が十分でない場合、高周波成分が低周波成分として混ざっている可能性があります。


8. エイリアシングを防ぐには

エイリアシングを防ぐための基本は、次の2つです。

1. 十分に高いサンプリング周波数で測る

観測したい最大周波数を $f_{\max}$ とすると、少なくとも

$$
f_s > 2 f_{\max}
$$

が必要です。

ただし実用上は、ノイズやフィルタ設計、ピーク検出、波形の形状なども考える必要があるため、2倍ぎりぎりではなく、もっと高いサンプリング周波数を使うことが多いです。


2. サンプリング前に高周波成分を取り除く

サンプリング前に、ナイキスト周波数を超える高周波成分を取り除くフィルタを使うことがあります。

これをアンチエイリアシングフィルタと呼びます。

ポイントは、

サンプリング後に頑張るのではなく、サンプリング前に折り返しを防ぐ

ということです。

サンプリングされてしまった後では、低周波に見えている成分が、本当に低周波だったのか、高周波が折り返したものだったのかを区別できない場合があります。


まとめ

この記事では、サンプリング周波数とエイリアシングについて、Pythonで作ったダミーデータを使って確認しました。

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

  • サンプリングとは、連続的な信号を一定時間ごとに測ること
  • サンプリング周波数 $f_s$ は、1秒間に何回測るかを表す
  • ナイキスト周波数は $f_s/2$
  • ナイキスト周波数を超える成分は、低い周波数に折り返して見える
  • 10 Hzでサンプリングすると、7 Hzの信号は3 Hzの信号のように見える
  • サンプリング後にFFTをしても、エイリアシングした成分は元に戻せない
  • 実データでは、十分なサンプリング周波数とアンチエイリアシングが重要

一言でまとめると、

データ解析の前に、そもそも正しく測れているかを確認することが大事

です。

信号処理や時系列解析では、解析手法そのものよりも、サンプリング条件が結果を大きく左右することがあります。


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

余裕があれば、以下を変更して遊んでみると理解が深まります。

  • f_signal = 7.08.09.0 に変えてみる
  • fs_bad = 101215 に変えてみる
  • 信号にノイズを加えてみる
  • 複数の正弦波を足してみる
  • 表示時間 duration を長くしてみる

例えば、信号にノイズを加えるなら、次のようにできます。

noise_level = 0.2
x_noisy = np.sin(2*np.pi*f_signal*t_ref) + noise_level * np.random.randn(len(t_ref))

エイリアシングは、ノイズがなくても起きます。
ノイズがあると、さらに見分けにくくなります。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?