サンプリング定理とエイリアシング:ナイキスト周波数とアンチエイリアシングフィルタ

サンプリング定理とナイキスト周波数の理論・Python実装を scipy.signal.decimate / resample / butter で解説。エイリアシング(折返し雑音)が起きる仕組み、アンチエイリアシングフィルタ設計、帯域通過標本化、sinc 再構成、オーバーサンプリングと整数倍/有理倍リサンプリングまで体系化。

はじめに

マイクで録音された音声、心電図センサーが計測した生体信号、レーダーが受信したRF波――これらはすべてアナログ信号です。現代のデジタル信号処理システムはこれらを扱うために、まずアナログ信号を離散的な数値列に変換する必要があります。この変換プロセスが標本化(サンプリング)であり、その根幹をなす理論がナイキスト-シャノンのサンプリング定理です。

サンプリング定理は「どれだけの速さでサンプルを取得すれば、元のアナログ信号を完全に復元できるか」という問いに対する厳密な答えを与えます。この定理を正しく理解しないと、エイリアシングと呼ばれる深刻な信号歪みが生じ、50Hzの電源ノイズが10Hzに見えたり、高速回転するプロペラがゆっくり逆回転しているように見えたりする現象が起こります。

本記事では、サンプリングの数学的モデルからナイキスト-シャノン定理の証明、エイリアシングの仕組み、アンチエイリアシングフィルタの設計、そしてPythonによる完全な実装まで体系的に解説します。関連する周波数解析の基礎については FFTの仕組みとPython実装 も合わせてご参照ください。

アナログ信号のデジタル化(A/D変換)

サンプリングの必要性

アナログ-デジタル変換器(ADC)は、連続時間信号 \(x(t)\) を等間隔の時刻 \(t = nT_s\) (\(n \in \mathbb{Z}\) )でサンプリングし、離散値列 \(x[n] = x(nT_s)\) を生成します。ここで \(T_s\) はサンプリング周期、その逆数 \(f_s = 1/T_s\) がサンプリング周波数(単位: Hz)です。

実際のADCは量子化も行いますが、本記事では理想的なサンプリング(量子化誤差なし)を仮定し、純粋にサンプリング周波数の選択が信号品質に与える影響に焦点を当てます。

インパルス列による数学的モデル

理想サンプリングは、アナログ信号 \(x(t)\) にディラックのインパルス列(シャー関数)\(s(t)\) を乗算する操作として数学的にモデル化されます。

\[s(t) = \sum_{n=-\infty}^{\infty} \delta(t - nT_s) \tag{1}\]

サンプリングされた信号 \(x_s(t)\) は次のようになります。

\[x_s(t) = x(t) \cdot s(t) = \sum_{n=-\infty}^{\infty} x(nT_s)\, \delta(t - nT_s) \tag{2}\]

この表現の重要な点は、\(x_s(t)\) が実数値列 \(x[n] = x(nT_s)\) を完全に保持していることです。連続時間信号から離散時間信号への情報の「圧縮」です。

サンプリングのスペクトルへの影響

式 \((2)\) をフーリエ変換すると、乗算は畳み込みになります。時間領域での積はスペクトル領域での畳み込みに対応するため:

\[X_s(f) = X(f) * S(f) \tag{3}\]

インパルス列のフーリエ変換は再びインパルス列です。

\[S(f) = f_s \sum_{k=-\infty}^{\infty} \delta(f - kf_s) \tag{4}\]

これを式 \((3)\) に代入すると:

\[X_s(f) = f_s \sum_{k=-\infty}^{\infty} X(f - kf_s) \tag{5}\]

式 \((5)\) が示す本質:サンプリング後のスペクトルは、元のスペクトル \(X(f)\) をサンプリング周波数 \(f_s\) の整数倍だけシフトしたものをすべて重ね合わせた(スペクトルの複製)ものです。

ナイキスト-シャノンのサンプリング定理

定理の陳述

定理(Nyquist-Shannon):帯域制限信号 \(x(t)\) (\(|f| > f_{\max}\) で \(X(f) = 0\) )は、サンプリング周波数が次の条件を満たす場合に、サンプル列 \(\{x[n]\}\) から完全に復元できる。

\[f_s > 2 f_{\max} \tag{6}\]

この条件を満たす最小サンプリング周波数 \(2f_{\max}\) をナイキストレートと呼びます。

証明のスケッチ

式 \((5)\) から、スペクトル複製が互いに重ならない条件を考えます。\(k=0\) の複製(元のスペクトル)の帯域幅は \([-f_{\max}, f_{\max}]\) です。隣の複製(\(k=\pm 1\) )は \([f_s - f_{\max}, f_s + f_{\max}]\) に位置します。

これらが重なり合わないための条件は:

\[f_{\max} < f_s - f_{\max}\] \[\Rightarrow f_s > 2f_{\max} \tag{7}\]

この条件が満たされれば、\(X_s(f)\) から理想ローパスフィルタ(カットオフ周波数 \(f_s/2\) 、ゲイン \(1/f_s\) )を適用することで \(X(f)\) を完全に取り出せます。

sinc補間による完全復元

時間領域での完全復元はsinc補間で表されます。

\[x(t) = \sum_{n=-\infty}^{\infty} x[n] \cdot \text{sinc}(f_s t - n) \tag{8}\]

ここで \(\text{sinc}(u) = \sin(\pi u) / (\pi u)\) はsinc関数です。各サンプル \(x[n]\) は、時刻 \(n/f_s\) を中心とするsinc関数で重み付けされ、これらの無限和が元の連続信号を再構成します。

式 \((8)\) の意味:デジタル信号からアナログ信号を復元する理想D/A変換は、sinc関数を基底とする補間です。これが周波数領域での理想LPFと等価であることはフーリエ変換から明らかです。

エッジケース:なぜ不等号は「以上」ではなく「超える」なのか

式 \((6)\) の条件 \(f_s > 2f_{\max}\) が狭義不等号(\(\geq\) ではなく \(>\) )である理由は、しばしば見過ごされる重要な注意点です。\(f_s = 2f_{\max}\) ちょうど(臨界サンプリング)の場合に何が起こるかを、最大周波数成分の位相に注目して確認します。

\(f_{\max} = 100\) Hz の余弦波と正弦波を、ちょうど \(f_s = 2f_{\max} = 200\) Hz でサンプリングすると:

import numpy as np

fmax = 100
fs = 2 * fmax  # 臨界サンプリング(等号ちょうど)

n = np.arange(8)
t = n / fs

x_cos = np.cos(2 * np.pi * fmax * t)
x_sin = np.sin(2 * np.pi * fmax * t)
print("cos samples:", np.round(x_cos, 6))
print("sin samples:", np.round(x_sin, 10))

実行結果:

cos samples: [ 1. -1.  1. -1.  1. -1.  1. -1.]
sin samples: [ 0.  0. -0.  0. -0. -0. -0. -0.]

余弦波は \(\{+1, -1, +1, -1, \dots\}\) という有効な(情報を保持した)サンプル列になりますが、正弦波は全サンプルが厳密に \(0\) になります。これは \(\sin(2\pi f_{\max} \cdot n/(2f_{\max})) = \sin(\pi n) = 0\) (すべての整数 \(n\) で)という恒等式そのものです。位相が90度ずれているだけの、同じ周波数・同じ振幅の信号が、片方は完全に復元可能でもう片方は跡形もなく消失します。これは深刻な問題です。信号の位相を知らない限り、\(f_s = 2f_{\max}\) での標本化は元信号の振幅を正しく保証しません。この位相依存の縮退があるからこそ、サンプリング定理は狭義不等号 \(f_s > 2f_{\max}\) を要求し、実務上は安全マージンを取ってナイキストレートより十分高い周波数でサンプリングします。

エイリアシング

エイリアシングの発生メカニズム

\(f_s < 2f_{\max}\) の場合、スペクトルの複製が重なり合います(エイリアシング)。式 \((5)\) で \(k=0\) と \(k=-1\) の複製を考えると:

\[X_{\text{alias}}(f) \ni X(f) + X(f + f_s) \tag{9}\]

この重なりは情報の不可逆な混合であり、もはや元の信号を分離して復元することはできません。

エイリアス周波数の計算

周波数 \(f\) の成分がサンプリング周波数 \(f_s\) でサンプリングされた場合、観測されるエイリアス周波数は次の式で計算されます。

\[f_\text{alias} = \left| f - \text{round}\left(\frac{f}{f_s}\right) \cdot f_s \right| \tag{10}\]

例として、\(f_s = 1000\) Hz で \(f = 1300\) Hz の信号をサンプリングすると:

\[f_\text{alias} = |1300 - \text{round}(1300/1000) \times 1000| = |1300 - 1000| = 300 \text{ Hz}\]

1300Hzの信号が300Hzとして観測されます。

ウェゴンホイール効果

映画やビデオで高速回転する車輪が逆回転しているように見える「ウェゴンホイール効果(Wagon Wheel Effect)」は、エイリアシングの直感的な例です。カメラのフレームレート(サンプリング周波数)より速く回転すると、ナイキスト条件が破られ、折り返しが起きて逆方向に見えます。

具体的に、フレームレート \(f_s = 24\) fps のカメラで 1 回転/フレーム より速い 1.1 回転/フレーム(\(f = 26.4\) Hz)の車輪を撮影した場合:

\[f_\text{alias} = |26.4 - \text{round}(26.4/24) \times 24| = |26.4 - 24| = 2.4 \text{ Hz}\]

車輪は実際より100倍遅い 2.4 Hz(1/2.4秒で1回転)でゆっくり回転しているように見えます。

ナイキスト周波数とナイキストレートの違い

この2つの用語は混同されやすいため、明確に区別します。

用語定義意味
ナイキスト周波数 \(f_N\)\(f_N = f_s / 2\)与えられたサンプリング周波数で表現できる最大周波数
ナイキストレート \(f_{NR}\)\(f_{NR} = 2 f_{\max}\)信号を完全に復元するために必要な最小サンプリング周波数
\[f_N = \frac{f_s}{2} \tag{11}\]

ナイキスト周波数はシステムの「能力」を表します。サンプリング周波数44.1 kHz(CDの規格)のシステムは最大22.05 kHzまでの信号を表現できます。人間の可聴域(約20 kHz)をカバーするために選ばれた値です。

ナイキストレートは信号の「要求」を表します。最高周波数20 kHzの音声信号を完全に復元したければ、少なくとも40 kHz以上でサンプリングする必要があります。

エイリアシングを防ぐ条件は「ナイキスト周波数 ≥ 信号の最高周波数」、すなわち「サンプリング周波数 ≥ ナイキストレート」と等価です。

アンチエイリアシングフィルタの設計

なぜADC前にLPFが必要か

実際のアナログ信号は有限帯域ではありません。センサーノイズ、EMI、電源ハムなど、意図しない高周波成分が常に存在します。これらをサンプリング前に除去しないと、ナイキスト周波数以上の全成分がエイリアシングを引き起こします。

**アンチエイリアシングフィルタ(AAF)**は、ADCの直前に配置するアナログローパスフィルタです。カットオフ周波数 \(f_c\) を \(f_s/2\) 以下に設定し、エイリアシングの原因となる高周波成分を事前に除去します。

カットオフ周波数の選択

理想的なAAFは遷移帯域が無限に急峻なブリック・ウォールフィルタですが、現実には実現不可能です。実用的な設計では次の考慮が必要です。

  1. 通過域上限 \(f_p\) :保持したい信号の最高周波数(例: 音声なら20 kHz)
  2. 阻止域下限 \(f_{stop}\) :エイリアシングを確実に防ぐ周波数(\(= f_s - f_p\) )
  3. 遷移帯域 \([f_p, f_{stop}]\) :この帯域でどれだけ急峻に減衰できるかがフィルタの次数を決定する

オーバーサンプリングによって \(f_s/2\) と \(f_{\max}\) の間のマージンを大きくすれば、より緩やかな(低次数の)アナログフィルタで十分になります。これはCDオーディオで \(f_{\max} = 20\) kHz に対して \(f_s = 44.1\) kHz(マージン: 4.1 kHz)が選ばれた設計理由の一つです。

バターワースフィルタやチェビシェフフィルタの詳細な設計については バターワースフィルタの設計 をご参照ください。

実務上の落とし穴:クロックジッタと非帯域制限信号

理論と実装の間には次のようなギャップが常に存在します。

  • サンプリングクロックジッタ:実際のADCのサンプリング間隔は \(T_s\) ちょうどではなく、微小な揺らぎ(ジッタ)\(\Delta t\) を含みます。振幅 \(A\) ・周波数 \(f\) の信号に対するジッタ起因の実効ノイズは概ね \(A \cdot 2\pi f \cdot \Delta t_{\text{rms}}\) で見積もられ、高周波・高振幅であるほどジッタの影響が拡大します。高精度ADC(16bit以上)では、ジッタが数十ps程度でも高域SNRを支配する要因になり得ます。
  • 有限時間信号は厳密には帯域制限されない:不確定性原理により、時間的に有限な長さを持つ信号は周波数領域で無限に広がるスペクトルを持ちます(両者を同時に有限にはできません)。したがって「帯域制限信号」はあくまで実用上の近似(無視できるレベルまでエネルギーが減衰した帯域を \(f_{\max}\) とみなす)であり、AAFは理論的な意味でも常に必要とされます。
  • 窓関数によるスペクトル漏れ:有限長で観測した信号のDFTは、真のスペクトルに窓関数の周波数応答が畳み込まれるため、帯域外に見かけ上のエネルギーが漏れ出します( 窓関数とPSDの理論 参照)。これもAAFで抑えきれない高周波成分が観測されうる一因です。

帯域通過標本化(サブナイキストサンプリング)

これまでの議論は信号が \([0, f_{\max}]\) に存在する低域信号(baseband signal)を前提としていました。しかし無線通信やレーダーでは、信号のエネルギーが \([f_L, f_H]\) (\(f_L > 0\) )という狭い帯域に集中している帯域通過信号(bandpass signal)を扱うことが多くあります。この場合、ナイキストレート \(2f_H\) より大幅に低い周波数で標本化しても、情報を失わずに復元できることがあります。これを帯域通過標本化(bandpass sampling、undersampling とも呼ばれる)といいます。

一般化されたサンプリング定理の導出

帯域幅を \(B = f_H - f_L\) とします。式 \((5)\) のスペクトル複製の考え方を帯域通過信号に適用すると、\(X(f)\) の複製(間隔 \(f_s\) )どうしが重ならない条件は、単純な \(f_s > 2f_H\) だけでなく、もっと緩い条件でも成立し得ます。

複製後の帯域 \([f_L - kf_s,\ f_H - kf_s]\) (\(k\) は整数)のいずれかが元の帯域 \([f_L, f_H]\) と重ならず、かつベースバンド(\(0\) 〜\(f_s/2\) )の同じ位置に一意に写像されればよいので、次の条件を満たす整数 \(n\) (\(1 \leq n \leq \lfloor f_H / B \rfloor\) )が存在すれば復元可能です。

\[ \frac{2f_H}{n} \leq f_s \leq \frac{2f_L}{n - 1} \quad (n \geq 2), \qquad f_s \geq 2f_H \ \ (n = 1) \tag{14} \]

\(n\) は「元の帯域の何番目のスペクトル複製がベースバンドに位置するか」を表す整数です。\(n=1\) の場合が通常のナイキスト条件 \(f_s \geq 2f_H\) に一致し、\(n \geq 2\) では、通常よりずっと低い \(f_s\) で標本化しても、\([f_L, f_H]\) の情報がベースバンド上の別の位置にエイリアシングして「移動」するだけで、失われないことを意味します。

数値例:RF帯域 [20, 25] MHzのサブナイキストサンプリング

\(f_L = 20\) MHz、\(f_H = 25\) MHz(帯域幅 \(B = 5\) MHz)のRF信号を考えます。\(n_{\max} = \lfloor f_H / B \rfloor = \lfloor 25/5 \rfloor = 5\) なので、\(n = 1, \dots, 5\) について式 \((14)\) を計算すると:

import numpy as np

fL, fH = 20e6, 25e6
B = fH - fL
n_max = int(np.floor(fH / B))

for n in range(1, n_max + 1):
    lo = 2 * fH / n
    if n == 1:
        print(f"n={n}: fs >= {lo/1e6:.3f} MHz(通常のナイキスト条件)")
    else:
        hi = 2 * fL / (n - 1)
        print(f"n={n}: fs in [{lo/1e6:.3f}, {hi/1e6:.3f}] MHz")

実行結果:

n=1: fs >= 50.000 MHz(通常のナイキスト条件)
n=2: fs in [25.000, 40.000] MHz
n=3: fs in [16.667, 20.000] MHz
n=4: fs in [12.500, 13.333] MHz
n=5: fs in [10.000, 10.000] MHz

通常のナイキスト条件(\(n=1\) )では \(f_s \geq 50\) MHz が必要ですが、\(n=4\) の帯域 \([12.5, 13.333]\) MHz を選べば、理論下限に近い \(f_s \approx 13\) MHz(\(2B = 10\) MHz にかなり近い)で標本化しても情報が失われません。ADCのサンプリングレートを\(1/4\) 近くまで下げられることになり、消費電力・データレートの観点で大きな利点になります。これがソフトウェア無線(SDR)のIF(中間周波数)サンプリングで広く使われる技術です。

複素(I/Q)サンプリングとの関係

上記の実数信号の場合、スペクトルは正負の周波数に対称(エルミート対称)なので、帯域幅 \(B\) の情報を表現するには実質的に \(2B\) の帯域を消費してしまいます。信号を複素ベースバンド信号(I/Q信号)としてダウンコンバートしてから標本化すれば、スペクトルは片側(\([0, B]\) または \([-B/2, B/2]\) )のみとなり、必要なサンプリングレートは \(f_s > B\) まで下げられます。これは実数サンプリングの半分で済み、無線通信の受信機(SDR、5G基地局など)で複素サンプリングが標準的に使われる理由です。

Pythonによる実装

実装1:エイリアシングの可視化

まず、エイリアシングが発生する状況を時間領域と周波数領域の両方で可視化します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, sosfilt, resample, freqz

# --- パラメータ設定 ---
f_signal = 130       # 信号周波数 [Hz]
fs_high   = 1000     # 十分高いサンプリング周波数(エイリアシングなし)
fs_low    = 200      # 不十分なサンプリング周波数(エイリアシング発生)
duration  = 0.1      # 表示時間 [s]

# --- 元の信号(高サンプリングレートで生成)---
t_high = np.arange(0, duration, 1 / fs_high)
x_orig = np.sin(2 * np.pi * f_signal * t_high)

# --- 低サンプリングレートでサンプリング ---
t_low = np.arange(0, duration, 1 / fs_low)
x_alias = np.sin(2 * np.pi * f_signal * t_low)

# エイリアス周波数を計算式で確認
f_alias = abs(f_signal - round(f_signal / fs_low) * fs_low)
print(f"信号周波数: {f_signal} Hz")
print(f"サンプリング周波数: {fs_low} Hz (ナイキスト: {fs_low//2} Hz)")
print(f"エイリアス周波数: {f_alias} Hz")
# => エイリアス周波数: 70 Hz

# --- プロット(時間領域) ---
fig, axes = plt.subplots(2, 1, figsize=(10, 6))

axes[0].plot(t_high * 1000, x_orig, 'b-', label=f'元信号 ({f_signal} Hz)', lw=1.5)
axes[0].plot(t_low * 1000, x_alias, 'ro', ms=8, label=f'低fs={fs_low}Hzでのサンプル')
# エイリアス波形を重ねて表示
t_fine = np.linspace(0, duration, 2000)
axes[0].plot(t_fine * 1000,
             np.sin(2 * np.pi * f_alias * t_fine),
             'r--', lw=1.5, label=f'エイリアス波形 ({f_alias} Hz)')
axes[0].set_xlabel('Time [ms]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('エイリアシング:時間領域')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# --- FFTによる周波数領域の比較 ---
N_high = len(x_orig)
N_low  = len(x_alias)

X_high = np.fft.rfft(x_orig)
X_low  = np.fft.rfft(x_alias)

freqs_high = np.fft.rfftfreq(N_high, 1 / fs_high)
freqs_low  = np.fft.rfftfreq(N_low,  1 / fs_low)

axes[1].plot(freqs_high, 2 / N_high * np.abs(X_high),
             'b-', label=f'fs={fs_high} Hz(エイリアシングなし)')
axes[1].plot(freqs_low, 2 / N_low * np.abs(X_low),
             'r-', label=f'fs={fs_low} Hz(エイリアシングあり)')
axes[1].axvline(fs_low / 2, color='k', ls='--', label=f'ナイキスト周波数 ({fs_low//2} Hz)')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('エイリアシング:周波数領域')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

実行結果:

信号周波数: 130 Hz
サンプリング周波数: 200 Hz (ナイキスト: 100 Hz)
エイリアス周波数: 70 Hz

130Hzの信号が200Hzでサンプリングされることで70Hzとして現れることが実測値でも確認できます。時間領域(上)では低サンプリングレートの点列(赤丸)が70Hzの波形(赤破線)にきれいに乗っており、周波数領域(下)でもfs=1000Hzでは130Hzにピークがあるのに対し、fs=200Hzでは70Hzにピークが出現しています。

エイリアシングの時間領域・周波数領域比較(130Hzの信号がfs=200Hzで70Hzに折り返る)

実装2:スペクトル折り返しの全貌

サンプリング周波数に対して様々な周波数の信号がどの周波数に折り返されるかを系統的に示します。

import numpy as np
import matplotlib.pyplot as plt

fs = 1000   # サンプリング周波数 [Hz]
fn = fs / 2  # ナイキスト周波数

# 様々な入力周波数に対するエイリアス周波数を計算
f_input = np.linspace(0, 3 * fs, 3000)
f_alias = np.abs(f_input - np.round(f_input / fs) * fs)

plt.figure(figsize=(10, 5))
plt.plot(f_input, f_alias, 'b-', lw=2)
plt.axhline(fn, color='r', ls='--', label=f'ナイキスト周波数 = {fn:.0f} Hz')
plt.fill_between([0, fn], [0, 0], [fn, fn],
                 alpha=0.1, color='green', label='エイリアスなし領域')
plt.xlabel('入力周波数 [Hz]')
plt.ylabel('観測周波数(エイリアス周波数)[Hz]')
plt.title(f'周波数折り返し特性 (fs = {fs} Hz)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(0, fn * 1.1)
plt.tight_layout()
plt.show()

周波数折り返し特性:三角波状のパターンがfs=1000Hzを周期に繰り返す

このグラフの三角波状パターンがナイキスト周波数での「折り返し」を視覚的に示します。\(0\) から \(f_s/2 = 500\) Hz の間の信号(緑の領域)はそのまま観測されますが、\(f_s/2\) を超えると折り返しが始まり、\(f_s = 1000\) Hz でエイリアス周波数はいったん \(0\) Hz に戻り、\(3f_s/2 = 1500\) Hz で再び \(f_s/2 = 500\) Hz から折り返すという、周期 \(f_s\) の周期的なパターンが繰り返されます。これは式 \((5)\) のスペクトル複製が周波数軸上で \(f_s\) 間隔に並ぶことの直接的な帰結です。

実装3:アンチエイリアシングフィルタの効果

ADC前にアンチエイリアシングフィルタを適用した場合と適用しない場合を比較します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, sosfilt, sosfreqz

np.random.seed(42)

# --- 帯域外成分を含む信号の生成 ---
fs_orig = 10000   # 元の高サンプリングレート [Hz]
fs_target = 1000  # ダウンサンプリング後のレート [Hz]
M = fs_orig // fs_target  # ダウンサンプリング比

duration = 0.5
t = np.arange(0, duration, 1 / fs_orig)

# 帯域内信号(100Hz)+ 帯域外信号(600Hz: ダウン後にエイリアシング)
x = np.sin(2 * np.pi * 100 * t) + 0.5 * np.sin(2 * np.pi * 600 * t)

# --- アンチエイリアシングフィルタの設計(バターワース LPF)---
# カットオフ: ナイキスト周波数の0.9倍
fc = (fs_target / 2) * 0.9  # = 450 Hz
sos = butter(8, fc, btype='low', fs=fs_orig, output='sos')
x_filtered = sosfilt(sos, x)

# --- ダウンサンプリング(AAFなし vs あり)---
x_down_noaaf = x[::M]          # AAFなし
x_down_aaf   = x_filtered[::M] # AAFあり

# --- 周波数応答の比較 ---
N_down = len(x_down_noaaf)
freqs_down = np.fft.rfftfreq(N_down, 1 / fs_target)

X_noaaf = np.fft.rfft(x_down_noaaf)
X_aaf   = np.fft.rfft(x_down_aaf)

# --- 400Hz(エイリアス周波数)での振幅を数値で確認 ---
idx400 = np.argmin(np.abs(freqs_down - 400))
amp_noaaf_400 = 2 / N_down * np.abs(X_noaaf[idx400])
amp_aaf_400   = 2 / N_down * np.abs(X_aaf[idx400])
print(f"400Hzでの振幅(AAFなし): {amp_noaaf_400:.4f}")
print(f"400Hzでの振幅(AAFあり): {amp_aaf_400:.6f}")
print(f"エイリアス周波数(400Hz)での抑制量: {20*np.log10(amp_noaaf_400/amp_aaf_400):.2f} dB")

fig, axes = plt.subplots(3, 1, figsize=(10, 10))

# フィルタの周波数応答(sos形式には sosfreqz を使う。freqz(sos, ...) はb,a係数専用でsos非対応)
w, h = sosfreqz(sos, worN=2048, fs=fs_orig, whole=False)
idx600 = np.argmin(np.abs(w - 600))
print(f"フィルタの600Hzでのゲイン: {20*np.log10(np.abs(h[idx600])+1e-12):.2f} dB")
axes[0].plot(w, 20 * np.log10(np.abs(h) + 1e-12), 'b-', lw=2)
axes[0].axvline(fc, color='r', ls='--', label=f'カットオフ fc = {fc:.0f} Hz')
axes[0].axvline(fs_target / 2, color='g', ls='--',
                label=f'ナイキスト周波数 = {fs_target//2} Hz')
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_title('アンチエイリアシングフィルタの周波数応答(8次バターワース)')
axes[0].set_xlim(0, 2000)
axes[0].set_ylim(-80, 5)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# AAFなしのスペクトル
axes[1].plot(freqs_down, 20 * np.log10(2 / N_down * np.abs(X_noaaf) + 1e-12),
             'r-', lw=1.5)
axes[1].axvline(fs_target / 2, color='k', ls='--', alpha=0.5)
axes[1].set_title('ダウンサンプリング後のスペクトル(AAFなし)— 600Hzが400Hzにエイリアス')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Magnitude [dB]')
axes[1].set_xlim(0, fs_target / 2)
axes[1].grid(True, alpha=0.3)

# AAFありのスペクトル
axes[2].plot(freqs_down, 20 * np.log10(2 / N_down * np.abs(X_aaf) + 1e-12),
             'b-', lw=1.5)
axes[2].axvline(fs_target / 2, color='k', ls='--', alpha=0.5)
axes[2].set_title('ダウンサンプリング後のスペクトル(AAFあり)— エイリアスを約21dB抑制')
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_ylabel('Magnitude [dB]')
axes[2].set_xlim(0, fs_target / 2)
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

実行結果:

400Hzでの振幅(AAFなし): 0.5000
400Hzでの振幅(AAFあり): 0.045634
エイリアス周波数(400Hz)での抑制量: 20.79 dB
フィルタの600Hzでのゲイン: -20.46 dB

600Hzの成分はダウンサンプリング後のナイキスト周波数 \(f_s/2 = 500\) Hz を超えているため、AAFなしでは \(|600 - 1000| = 400\) Hzにエイリアシングします。実測では、AAFなしの場合の400Hz成分の振幅は理論値どおり \(0.5\) (元の600Hz成分の振幅と一致)ですが、8次バターワースフィルタ(カットオフ450Hz)を適用すると600Hzでのフィルタゲインが \(-20.46\) dB まで減衰しており、その結果400Hzでの残留振幅は \(0.045634\) まで下がります。これはおよそ \(20.79\) dB の抑制に相当し、完全にゼロにはならないものの、実用上十分なレベルまでエイリアス成分が抑えられていることが数値で確認できます。フィルタの次数を上げる、あるいはオーバーサンプリング比を増やしてカットオフとナイキスト周波数のマージンを広げることで、さらなる抑制が可能です。

アンチエイリアシングフィルタの効果:フィルタ周波数応答(上)、AAFなしのスペクトル(中)、AAFありのスペクトル(下)

実装4:scipy.signal.resampleによるリサンプリング

scipy.signal.resample はFFTベースのリサンプリングを実装しており、アップサンプリングとダウンサンプリングの両方に利用できます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import resample, butter, sosfilt

# --- 元信号の生成 ---
fs_orig = 1000   # 元のサンプリングレート
duration = 1.0
t_orig = np.arange(0, duration, 1 / fs_orig)
x_orig = (np.sin(2 * np.pi * 50 * t_orig)
          + 0.3 * np.sin(2 * np.pi * 120 * t_orig))

# --- ダウンサンプリング(1000 Hz → 400 Hz、デシメーション比 M=2.5)---
fs_down = 400
n_down = int(len(x_orig) * fs_down / fs_orig)
x_down = resample(x_orig, n_down)  # FFTベース:内部でAAFを適用
t_down = np.linspace(0, duration, n_down)

# --- アップサンプリング(元の1000 Hz → 4000 Hz)---
fs_up = 4000
n_up = int(len(x_orig) * fs_up / fs_orig)
x_up = resample(x_orig, n_up)      # 補間によるアップサンプリング
t_up = np.linspace(0, duration, n_up)

# --- スペクトル比較 ---
def spectrum(x, fs):
    N = len(x)
    X = np.fft.rfft(x)
    f = np.fft.rfftfreq(N, 1 / fs)
    return f, 2 / N * np.abs(X)

fig, axes = plt.subplots(3, 1, figsize=(10, 9))

f_orig, A_orig = spectrum(x_orig, fs_orig)
f_down, A_down = spectrum(x_down, fs_down)
f_up,   A_up   = spectrum(x_up,   fs_up)

axes[0].plot(f_orig, A_orig, 'b-')
axes[0].set_title(f'元信号のスペクトル (fs = {fs_orig} Hz)')
axes[0].set_ylabel('Amplitude')
axes[0].set_xlim(0, fs_orig / 2)
axes[0].grid(True, alpha=0.3)

axes[1].plot(f_down, A_down, 'r-')
axes[1].set_title(f'ダウンサンプリング後 (fs = {fs_down} Hz)')
axes[1].set_ylabel('Amplitude')
axes[1].set_xlim(0, fs_down / 2)
axes[1].grid(True, alpha=0.3)

axes[2].plot(f_up, A_up, 'g-')
axes[2].set_title(f'アップサンプリング後 (fs = {fs_up} Hz)')
axes[2].set_ylabel('Amplitude')
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_xlim(0, fs_up / 2)
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"元信号:             {len(x_orig)} サンプル ({fs_orig} Hz)")
print(f"ダウンサンプリング: {len(x_down)} サンプル ({fs_down} Hz)")
print(f"アップサンプリング: {len(x_up)} サンプル ({fs_up} Hz)")

実行結果:

元信号:             1000 サンプル (1000 Hz)
ダウンサンプリング: 400 サンプル (400 Hz)
アップサンプリング: 4000 サンプル (4000 Hz)

scipy.signal.resample はFFTによってスペクトルを計算し、必要なサンプル数にトリミング(ダウンサンプリング)またはゼロ埋め(アップサンプリング)してIDFFTを適用します。この手法は内部的にアンチエイリアシングを処理しますが、任意比のリサンプリングが可能という利点があります。

リサンプリングの比較:元信号(上)、400Hzへのダウンサンプリング後(中)、4000Hzへのアップサンプリング後(下)のスペクトル

図から、ダウンサンプリング後も50Hzと120Hzの2つの成分が正しく保持されている(ナイキスト周波数200Hzを下回っているため)ことが視覚的に確認できます。アップサンプリング後も同じ2つの成分がそのままの周波数位置に保たれており、スペクトル内容を変えずにサンプル数だけが変化していることがわかります。

実装5:帯域通過標本化(サブナイキストサンプリング)の数値検証

前節で導出した式 \((14)\) が実際に成り立つことを、\([20, 25]\) MHz 帯にある2つのトーン(22 MHz, 24 MHz)を \(f_s = 13\) MHz(\(n=4\) の帯域内)で標本化して確認します。

import numpy as np

# --- 帯域通過信号の生成:[20, 25] MHz帯に2トーン ---
f1, f2 = 22e6, 24e6       # 帯域内の2つのトーン周波数
fs = 13e6                 # サブナイキストサンプリング周波数(n=4の帯域 [12.5, 13.333] MHz内)

duration = 20e-6           # 20 マイクロ秒
t = np.arange(0, duration, 1 / fs)
x = np.sin(2 * np.pi * f1 * t) + 0.6 * np.sin(2 * np.pi * f2 * t)

# --- 理論的なエイリアス周波数(式(10)と同じ折り返し公式)---
def alias(f, fs):
    return abs(f - round(f / fs) * fs)

a1, a2 = alias(f1, fs), alias(f2, fs)
print(f"理論エイリアス: {f1/1e6} MHz -> {a1/1e6:.3f} MHz, {f2/1e6} MHz -> {a2/1e6:.3f} MHz")

# --- FFTで実際に観測される周波数を確認 ---
N = len(x)
X = np.fft.rfft(x)
freqs = np.fft.rfftfreq(N, 1 / fs)
mag = 2 / N * np.abs(X)

top_idx = np.argsort(mag)[::-1][:2]
for i in sorted(top_idx):
    print(f"観測ピーク: f={freqs[i]/1e6:.4f} MHz, 振幅={mag[i]:.4f}")

実行結果:

理論エイリアス: 22.0 MHz -> 4.000 MHz, 24.0 MHz -> 2.000 MHz
観測ピーク: f=2.0000 MHz, 振幅=0.6000
観測ピーク: f=4.0000 MHz, 振幅=1.0000

理論式 \((10)\) から予測したエイリアス周波数(22 MHz → 4.0 MHz、24 MHz → 2.0 MHz)と、実際にFFTで観測されたピーク周波数が完全に一致しています。振幅も元の係数(1.0 と 0.6)をそのまま保持しており、\(f_s = 13\) MHz という低いレートでの標本化でも情報が失われていないことが数値的に確認できます。

帯域通過標本化:fs=13MHzで[20,25]MHz帯の2トーンを標本化すると、ベースバンド上の2MHzと4MHzに一意に写像される

アップサンプリングとダウンサンプリング

デシメーション(ダウンサンプリング)

整数比 \(M\) のダウンサンプリング(デシメーション)は次の手順で行います。

  1. アンチエイリアシングLPF:カットオフ \(f_s/(2M)\) でフィルタリング
  2. \(M\) サンプルに1つを保持:\(y[n] = x[Mn]\)

手順1を省略すると、新しいナイキスト周波数 \(f_s/(2M)\) を超える成分がエイリアシングを起こします。scipy.signal.decimate はこの2段階を自動処理します。

\[y[n] = x_{\text{filtered}}[Mn] \tag{12}\]

補間(アップサンプリング)

整数比 \(L\) のアップサンプリング(補間)は次の手順で行います。

  1. ゼロ挿入:\(L-1\) 個のゼロを各サンプル間に挿入
  2. イメージ除去LPF:カットオフ \(f_s/(2L)\) (ゲイン \(L\) )でフィルタリング

ゼロ挿入だけではスペクトルに「イメージ」(スペクトルの複製)が残るため、LPFで除去します。

\[x_L[n] = \begin{cases} x[n/L] & n \text{ が } L \text{ の倍数} \\ 0 & \text{その他} \end{cases} \tag{13}\]

ポリフェーズフィルタリング

実用的なリサンプリング実装では、ポリフェーズフィルタが用いられます。ポリフェーズ分解により、フィルタを \(L\) 個のサブフィルタに分解し、ゼロ乗算を避けることで計算効率を大幅に改善します。計算量は素朴な実装の \(1/L\) になります。

scipy.signal.resample_poly はポリフェーズフィルタリングを実装した高効率なリサンプリング関数です。任意の有理数比 \(L/M\) のリサンプリングが可能です。

from scipy.signal import resample_poly

# 1000 Hz → 441 Hz(L=441, M=1000 のポリフェーズ)
x_poly = resample_poly(x_orig, 441, 1000)
print(f"ポリフェーズリサンプリング: {len(x_orig)}{len(x_poly)} サンプル")

実行結果:

ポリフェーズリサンプリング: 1000 → 441 サンプル

ポリフェーズ分解と resample_poly の内部実装の詳細(Type-1分解、Noble identitiesによる計算量削減など)は マルチレート信号処理の理論とPython実装 で深く扱っているので、そちらを参照してください。本記事ではサンプリング定理の帰結として、整数比・有理数比のレート変換にAAF/イメージ除去フィルタが必須である点までを扱います。

発展:線形再構成を超えるサンプリング理論

本記事で解説したナイキスト-シャノン定理は、線形(sinc)補間による帯域制限信号の復元を前提としています。近年はこの前提を緩めた研究が進んでいます。Najaf and Ongie (2024) “Towards a Sampling Theory for Implicit Neural Representations”( arXiv:2405.18410 )は、信号を単純なニューラルネットワーク(ReLU活性化・単一隠れ層の Implicit Neural Representation)で近似可能であるという仮定のもとで、低域フーリエ係数から連続領域の信号(画像)を完全復元するために必要なサンプル数を理論的に導出しています。線形補間を前提とする古典的サンプリング定理とは異なり、信号のクラスに関する非線形な事前知識(ニューラル表現で表現可能という制約)を使うことで、同じ復元精度をより少ないサンプル数で達成できる可能性を示す研究です。医療画像(MRI)の分野でも、Implicit Neural Representationを使ってサブナイキストレートで取得したk空間データから画像を復元する手法が実用化されつつあり、本記事で解説した古典的な線形サンプリング理論を補完する形で研究が進展しています。

まとめ

本記事で解説した内容を次の表にまとめます。

概念定義ポイント
サンプリング定理\(f_s > 2f_{\max}\)この条件が完全復元の必要十分条件
ナイキスト周波数\(f_N = f_s / 2\)システムが表現できる最大周波数
ナイキストレート\(f_{NR} = 2f_{\max}\)信号に要求される最小サンプリング周波数
エイリアス周波数\(f_\text{alias} = \|f - \text{round}(f/f_s)\cdot f_s\|\)折り返し後の観測周波数
アンチエイリアシングフィルタADC前のアナログLPFカットオフ \(\leq f_s/2\) 、ダウンサンプリング前に必須
デシメーションAAF → \(M\) サンプルに1つ保持scipy.signal.decimate
補間ゼロ挿入 → イメージ除去LPFscipy.signal.resample_poly
帯域通過標本化式 \((14)\) :\(2f_H/n \leq f_s \leq 2f_L/(n-1)\)\(f_s \ll 2f_H\) でも情報を保持できる(SDRのIFサンプリング)

サンプリング理論の理解は、DSPシステム設計の出発点です。 窓関数とPSDの理論ローパスフィルタの設計と比較 などと合わせて、信号処理パイプライン全体を俯瞰することを推奨します。

おすすめ書籍

はじめて学ぶディジタル・フィルタと高速フーリエ変換(三上直樹、CQ出版)

ディジタルフィルタと FFT の原理を基礎から丁寧に解説した定番の入門書です。本記事で扱った処理の理論的背景を体系的に学べます。

※ 上記は Amazon アソシエイトのリンクです。

関連記事

参考文献

  • Shannon, C. E. (1949). “Communication in the Presence of Noise”. Proceedings of the IRE, 37(1), 10-21.
  • Nyquist, H. (1928). “Certain Topics in Telegraph Transmission Theory”. Transactions of the AIEE, 47(2), 617-644.
  • Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
  • Proakis, J. G., & Manolakis, D. G. (2006). Digital Signal Processing (4th ed.). Prentice Hall.
  • Najaf, M., & Ongie, G. (2024). “Towards a Sampling Theory for Implicit Neural Representations”. 58th Asilomar Conference on Signals, Systems, and Computers. arXiv:2405.18410
  • SciPy Signal Processing documentation