ナイキスト周波数とは?エイリアシングが起きる仕組みを導出から理解する

映画やテレビで走っている車を見ていると、ホイールが止まって見えたり、ときには進行方向と逆にゆっくり回って見えたりします。実際の車輪は猛烈な速さで前向きに回っているのに、カメラが 1 秒に 24 枚しか撮っていないせいで、私たちの目には「ゆっくり逆回転する車輪」という存在しない現象が見えてしまう。これがエイリアシング(aliasing、折り返し)です。

同じことが、音声のA/D変換にも、オシロスコープの波形取得にも、レーダーの距離計測にも、CMOSイメージセンサの画像にも起こります。しかも厄介なことに、いったん折り返ってしまった成分は、あとからディジタル処理でどれだけ頑張っても原理的に取り除けません。エイリアシングは「あとで直せるノイズ」ではなく「情報が壊れる事故」なのです。

この事故が起きるか起きないかを分ける境界線が ナイキスト周波数 $f_N = f_s/2$ です。サンプリング周波数 $f_s$ の半分、それだけの話に見えますが、「なぜ半分なのか」「半分を超えると具体的に何 Hz に化けるのか」「化けないようにするにはどんなフィルタが要るのか」を説明できる人は意外と多くありません。

この記事では、ナイキスト周波数を丸暗記の公式ではなく、周波数領域でスペクトルが複製されるという一枚の絵から導出します。その絵さえ手に入れば、折り返し周波数の計算式も、アンチエイリアシングフィルタに必要な阻止量も、オーバーサンプリングが効く理由も、すべて同じ絵の帰結として説明できます。

まず、これから扱う現象の全体像を一枚の絵で押さえておきましょう。

粗いストロボで速い波を観測すると遅い波に化けるというエイリアシングの概念図

灰色の細かい波が本物の信号(1300 Hz)、青い丸が 1000 Hz のストロボ(サンプラー)で拾った値、赤い破線が観測側に見える偽物(300 Hz)です。注目すべきは、青い丸が灰色の波の上にも赤い破線の上にも同時に乗っているという一点に尽きます。標本値だけを渡された受信側には、元が速い波だったのか遅い波だったのかを区別する材料が一切残っていない — これがエイリアシングという事故の正体で、この記事はこの絵を数式で裏づけ、境界線がどこにあるかを求めていく作業です。

具体的な応用としては、少なくとも次の 2 つを念頭に置いてください。

  • ディジタルオーディオのA/D変換設計 — 「なぜ CD は 44.1 kHz なのか」「なぜ現代のADCは 192 kHz でオーバーサンプリングしてから間引くのか」は、この記事で導く必要次数の式でそのまま説明できます。
  • 計測器・レーダー・SDRのフロントエンド設計 — サンプリング周波数を決めるとき、同時に「アンチエイリアシングフィルタの遷移帯域と阻止量」を決めないと設計は完結しません。逆に、意図的に折り返しを使う帯域通過サンプリング(アンダーサンプリング)も、同じ式から設計できます。

本記事の内容

  • ナイキスト周波数の直感的な意味と、「ナイキストレート」との区別
  • インパルス列との積としてサンプリングを書き、スペクトルが $f_s$ 間隔で複製されることの導出
  • 複製が重ならない条件から $f_s > 2f_{\max}$、すなわち $f_{\max} < f_N$ を導く
  • 折り返し後の見かけ周波数 $f_a = |f – \mathrm{round}(f/f_s)\,f_s|$ の導出と、位相が反転する理由
  • 周波数を掃引したときに観測周波数が三角波状に折り返す様子
  • アンチエイリアシングフィルタに必要な阻止量と次数の見積り、オーバーサンプリングの効果
  • Python による可視化(低周波への化け方・スペクトログラムでの折り返し・LPFの必要次数)

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

ナイキスト周波数とは — まず直感から

サンプリングとは、連続時間の波形をストロボで照らして、そのときの値だけを記録する作業です。ストロボの点滅が速ければ波形の形は忠実に残りますが、点滅が遅いと「点と点のあいだで何が起きたか」がわからなくなります。ここで自然に浮かぶ問いは、「1 周期あたり最低何点あれば、その正弦波を正弦波だと言い切れるのか」 です。

紙の上でやってみるとすぐわかります。正弦波の 1 周期に点が 1 点しかなければ、その点を通る正弦波は無数に描けます。2 点あると、山と谷を 1 回ずつ捉えられるので、ようやく「上がって下がる」という振動の存在がわかる。3 点、4 点と増やせば、波形の形はどんどん確からしくなります。逆に言えば、1 周期あたり 2 点が、振動の存在を主張できるぎりぎりの下限です。

1周期あたり1点・2点・4点で標本化したときに波形がどこまでわかるかの比較

左のパネルでは 1 周期に 1 点しかなく、標本値が毎回同じ高さになるため、直流(破線)でも 2 Hz の正弦波(点線)でも同じ標本列が得られてしまいます。中央の 1 周期 2 点になると、標本値が必ず正負に振れるので「何かが振動している」ことだけは主張できます。右の 1 周期 4 点まで来ると、山と谷の位置と傾きの向きまで見え、元の波形の形がほぼ再現されます。2 点というのは「形がわかる」下限ではなく「振動の存在がわかる」下限だ、という点が読み取りどころです。

この「1 周期あたり 2 点」を周波数の言葉に翻訳します。サンプリング周波数 $f_s$ で 1 秒間に $f_s$ 点取れるので、周波数 $f$ の正弦波 1 周期あたりに取れる点数は $f_s / f$ 点です。これが 2 点以上であればよいので

$$ \frac{f_s}{f} > 2 \quad \Longleftrightarrow \quad f < \frac{f_s}{2} $$

となります。この右辺 $f_s/2$ が ナイキスト周波数 です。

$$ \begin{equation} f_N \equiv \frac{f_s}{2} \end{equation} $$

ナイキスト周波数は「サンプリング周波数 $f_s$ で記録したときに、正しく記録できる信号周波数の上限」を表します。$f_s = 48\,\mathrm{kHz}$ なら $f_N = 24\,\mathrm{kHz}$、$f_s = 44.1\,\mathrm{kHz}$ なら $f_N = 22.05\,\mathrm{kHz}$ です。単位は「Hz」であり、サンプリング周波数と同じ次元の量である点を押さえてください。

ナイキストレートと混同しない

紛らわしいのが ナイキストレート(Nyquist rate)です。こちらは信号側の性質から決まる量で、「信号の最高周波数 $f_{\max}$ の 2 倍」、つまり

$$ f_{\text{Nyquist rate}} = 2 f_{\max} $$

と定義されます。必要なサンプリング周波数の下限を表す量です。

  • ナイキスト周波数 $f_N = f_s/2$ … サンプラー側の都合。「$f_s$ で取ると、何 Hz まで記録できるか」
  • ナイキストレート $2f_{\max}$ … 信号側の都合。「この信号を記録するには、何 Hz 以上で取るべきか」

同じ人名がついていますが、片方は $f_s$ を 2 で割り、もう片方は $f_{\max}$ を 2 倍します。設計の場面では「信号帯域 $f_{\max} = 20\,\mathrm{kHz}$ だからナイキストレートは 40 kHz、余裕を見て $f_s = 48\,\mathrm{kHz}$ にすると $f_N = 24\,\mathrm{kHz}$ でガードバンドが 4 kHz 取れる」という使い方をします。

ここまでは「1 周期に 2 点」という時間領域の直感で $f_s/2$ を出しました。しかしこの説明では、上限を超えたときに具体的に何が起きるのかが見えません。「情報が失われる」だけでなく「別の周波数の信号として出現する」という、もっと具体的で厄介な現象が起きます。それを説明できるのは周波数領域の描像だけです。次の節で、サンプリングという操作を数式に落として、周波数領域で何が起きているかを見ていきましょう。

サンプリングを数式で書く — インパルス列との積

サンプリングを数式で扱うとき、いちばん見通しがよいのは「連続信号に、等間隔に立った針を掛け算する」という描像です。針の位置でだけ値が残り、それ以外はゼロになる。この針の列を インパルス列(くし関数、Dirac comb)と呼びます。

$$ \begin{equation} \delta_{T_s}(t) = \sum_{n=-\infty}^{\infty} \delta(t – n T_s), \qquad T_s = \frac{1}{f_s} \end{equation} $$

サンプリングされた信号 $x_s(t)$ は、元の連続信号 $x(t)$ とこのインパルス列の積として書けます。

$$ \begin{equation} x_s(t) = x(t)\,\delta_{T_s}(t) = \sum_{n=-\infty}^{\infty} x(nT_s)\,\delta(t – nT_s) \end{equation} $$

最後の等号では、$\delta(t – nT_s)$ が $t = nT_s$ でしか値を持たないので $x(t)$ をその点の値 $x(nT_s)$ に置き換えられる、という $\delta$ 関数の性質(篩い出し)を使いました。つまり $x_s(t)$ は「高さが標本値 $x(nT_s)$ の針が $T_s$ 間隔で並んだもの」です。

連続信号とインパルス列の積としてサンプリングを表した3段の図

上段の連続信号 $x(t)$ に、中段の等間隔($T_s$)に立った針の列を掛けると、下段のように「元の波形の高さを持つ針の列」になります。針以外の場所は値がゼロになるので、下段には元の波形の情報が標本点でしか残っていません。ここで大事なのは、サンプリングという操作が「間引き」ではなく掛け算として書けたという一点です。掛け算にできたおかげで、次に見るように周波数領域では畳み込みという扱いやすい形に変わります。

なぜわざわざ $\delta$ 関数を持ち出すのでしょうか。理由は単純で、時間領域の積は周波数領域の畳み込みになるというフーリエ変換の性質を使いたいからです。掛け算で書けてさえいれば、周波数領域で何が起きるかは自動的に決まります。

インパルス列のフーリエ変換

まずインパルス列そのもののスペクトルを求めます。$\delta_{T_s}(t)$ は周期 $T_s$ の周期関数なので、複素フーリエ級数に展開できます。

$$ \delta_{T_s}(t) = \sum_{k=-\infty}^{\infty} c_k \, e^{j 2\pi k f_s t} $$

係数 $c_k$ は 1 周期分の積分で計算します。区間を $[-T_s/2,\, T_s/2]$ に取ると、この中にあるインパルスは $t=0$ のもの 1 本だけです。

$$ c_k = \frac{1}{T_s}\int_{-T_s/2}^{T_s/2} \delta_{T_s}(t)\, e^{-j2\pi k f_s t}\, dt = \frac{1}{T_s}\int_{-T_s/2}^{T_s/2} \delta(t)\, e^{-j2\pi k f_s t}\, dt $$

ここで $\delta$ 関数の篩い出し(積分は被積分関数の $t=0$ での値を拾う)を使うと $e^{-j2\pi k f_s \cdot 0} = 1$ なので

$$ c_k = \frac{1}{T_s} = f_s $$

となります。すべての $k$ について同じ値という点が重要です。したがって

$$ \delta_{T_s}(t) = f_s \sum_{k=-\infty}^{\infty} e^{j2\pi k f_s t} $$

次に、この両辺をフーリエ変換します。複素指数のフーリエ変換は $\mathcal{F}\{e^{j2\pi f_0 t}\} = \delta(f – f_0)$ ですから、$f_0 = k f_s$ を代入して和を取ると

$$ \begin{equation} \Delta(f) \equiv \mathcal{F}\{\delta_{T_s}(t)\} = f_s \sum_{k=-\infty}^{\infty} \delta(f – k f_s) \end{equation} $$

つまり、時間領域で $T_s$ 間隔に並んだインパルス列は、周波数領域でも $f_s$ 間隔に並んだインパルス列になる、というきれいな双対性が成り立ちます。時間間隔を詰める($T_s$ を小さくする)ほど、周波数側の間隔 $f_s$ は広がります。この「時間で詰めると周波数で広がる」感覚が、以降の議論の全体を支配します。

インパルス列の時間間隔を半分にすると周波数間隔が2倍に広がる双対性の図

上段($T_s = 0.1$ s)と下段($T_s = 0.05$ s)を見比べてください。時間軸で針の間隔を半分に詰めると、周波数軸では針の間隔が 10 Hz から 20 Hz へと 2 倍に広がっています。同時に周波数側の針の高さも $f_s$ に比例して伸びており、式 (4) の係数 $f_s$ がそのまま見えている形です。この「詰めると広がる」関係が、次に導くスペクトルの複製間隔をコントロールするつまみになります。

サンプリング後のスペクトル — 周期化

準備ができたので、$x_s(t) = x(t)\,\delta_{T_s}(t)$ をフーリエ変換します。積のフーリエ変換は畳み込みになるので

$$ X_s(f) = X(f) * \Delta(f) $$

右辺に式 (4) を代入します。畳み込みは線形なので和の外に出せて

$$ X_s(f) = X(f) * \left[ f_s \sum_{k=-\infty}^{\infty} \delta(f – k f_s) \right] = f_s \sum_{k=-\infty}^{\infty} \left[ X(f) * \delta(f – k f_s) \right] $$

ここで「$\delta(f – a)$ との畳み込みは、関数を $a$ だけ平行移動する操作」という性質を使います。すなわち $X(f) * \delta(f – kf_s) = X(f – kf_s)$ です。これを代入すると、最終的に

$$ \begin{equation} X_s(f) = f_s \sum_{k=-\infty}^{\infty} X(f – k f_s) \end{equation} $$

を得ます。これが本記事でいちばん大事な式です。日本語に訳すと、サンプリングすると、元のスペクトルの形が $f_s$ 間隔でコピー&ペーストされて無限に並ぶ(スペクトルの周期化)ということです。原本($k=0$)はそのまま残り、その左右に $\pm f_s$、$\pm 2f_s$、… ずらした複製がずらりと並びます。

サンプリングによって元のスペクトルがfs間隔で複製されて並ぶ様子

上段が元の信号のスペクトル $X(f)$($f_{\max} = 300$ Hz で帯域制限、幅は $2f_{\max}$)、下段がサンプリング後の $X_s(f)$ です。濃い色の $k=0$ が原本で、$\pm 1000$ Hz、$\pm 2000$ Hz、… に薄い色の複製が等間隔で並んでいるのが見えます。ここでは複製どうしのあいだに緑の矢印で示した隙間があるので、ナイキスト周波数 $f_N = 500$ Hz でスパッと切り出せば原本だけを取り戻せます。逆に言えば、この隙間が消えた瞬間に復元が不可能になる — 次節はその境界を式にする作業です。

別ルートでの確認 — DTFT が $f_s$ 周期になる理由

式 (5) は、標本列そのもののフーリエ変換からも同じ形で確認できます。式 (3) の右辺を直接フーリエ変換すると

$$ X_s(f) = \int_{-\infty}^{\infty} \left[\sum_n x(nT_s)\delta(t – nT_s)\right] e^{-j2\pi f t}\, dt = \sum_{n=-\infty}^{\infty} x(nT_s)\, e^{-j2\pi f n T_s} $$

これは標本列 $x[n] = x(nT_s)$ の離散時間フーリエ変換(DTFT)そのものです。ここで $f$ を $f + f_s$ に置き換えてみると、指数部の追加分は $e^{-j2\pi f_s n T_s} = e^{-j2\pi n} = 1$ なので値は変わりません。したがって

$$ X_s(f + f_s) = X_s(f) $$

が任意の $f$ で成り立ちます。離散信号のスペクトルは必ず $f_s$ 周期の周期関数である — これは式 (5) の言い換えであり、$\delta$ 関数を使わずに導ける事実でもあります。「サンプリングしたら周波数軸が $f_s$ で巻き取られる」と覚えておくと、以降の議論がすべて腑に落ちます。

ここまでで、サンプリングの周波数領域での正体が「複製の無限列」だとわかりました。となると、次に気になるのはその複製たちが互いにぶつからない条件です。ぶつからなければ元のスペクトルを切り出して復元でき、ぶつかれば混ざって復元できない。この条件を書き下すことが、そのままナイキスト周波数の導出になります。

複製が重ならない条件からナイキスト周波数を導く

元の信号 $x(t)$ が 帯域制限 されているとします。すなわち、ある $f_{\max}$ が存在して

$$ X(f) = 0 \quad \text{for} \quad |f| > f_{\max} $$

が成り立つとします。実信号のスペクトルは正負両側に広がるので、原本 $k=0$ の複製は周波数軸上で区間

$$ [-f_{\max},\; +f_{\max}] $$

を占めます。幅は $2f_{\max}$ です。

$k=1$ の複製は、この区間をそのまま $+f_s$ だけ平行移動したものなので、区間

$$ [\,f_s – f_{\max},\; f_s + f_{\max}\,] $$

を占めます。同様に $k=-1$ の複製は $[-f_s – f_{\max},\, -f_s + f_{\max}]$ です。

隣り合う複製が重ならない条件は、「原本の右端」が「$k=1$ の複製の左端」より小さいことです。式にすると

$$ f_{\max} < f_s - f_{\max} $$

両辺に $f_{\max}$ を足して整理すると

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

すなわち

$$ \begin{equation} f_s > 2 f_{\max} \end{equation} $$

が得られます。これが サンプリング定理(標本化定理)の条件です。$k=-1$ 側についても、原本の左端 $-f_{\max}$ が $k=-1$ の複製の右端 $-f_s + f_{\max}$ より大きいこと、つまり $-f_{\max} > -f_s + f_{\max}$ を要求すると、まったく同じ不等式に帰着します。左右対称なので条件は一つで済みます。

式 (6) の両辺を 2 で割ると

$$ \begin{equation} f_{\max} < \frac{f_s}{2} = f_N \end{equation} $$

信号の最高周波数がナイキスト周波数より小さいこと — これが、時間領域の「1 周期あたり 2 点」の直感に対応する、周波数領域での厳密な条件です。二つの導出が同じ答えにたどり着いたことで、$f_s/2$ という数字の意味がだいぶ立体的になったはずです。

帯域制限が足りる場合と足りない場合でスペクトルの複製が重なるかどうかを比べた図

上段は $f_{\max} = 400\,\mathrm{Hz} < f_N$ の場合で、複製どうしのあいだに隙間があり、$[-f_N, f_N]$ を切り出せば原本がそっくり手に入ります。下段は $f_{\max} = 700\,\mathrm{Hz} > f_N$ の場合で、$k = \pm 1$ の複製の裾が原本に食い込み、赤い太線で示した「実際に観測されるスペクトル」は原本と複製のになってしまっています。和になった時点で、どこまでが原本でどこからが複製かを分ける情報は失われており、どんなフィルタを掛けても元の三角形は取り出せません。式 $f_{\max} < f_s - f_{\max}$ が要求しているのは、まさにこの「裾が届かないだけの距離」です。

条件が満たされているとき何ができるか

条件 (7) が満たされていれば、$X_s(f)$ の中で原本 $k=0$ の複製は $[-f_N, f_N]$ の中に孤立して存在し、他の複製と一切混ざりません。そこで、遮断周波数 $f_N$ の理想ローパスフィルタ

$$ H(f) = \begin{cases} 1/f_s & |f| < f_N \\ 0 & |f| \ge f_N \end{cases} $$

を掛ければ、$X(f)$ をそっくり取り出せます。これを時間領域に戻すと、理想ローパスフィルタのインパルス応答が $\mathrm{sinc}$ 関数であることから、有名な sinc 補間 の式

$$ x(t) = \sum_{n=-\infty}^{\infty} x(nT_s)\,\mathrm{sinc}\!\left(\frac{t – nT_s}{T_s}\right) $$

が得られます。「標本値さえあれば、間の値も含めて元の連続波形が完全に復元できる」という、サンプリング定理のもう一つの顔です(sinc 補間の詳細はサンプリング定理の記事を参照してください)。

「ちょうど 2 倍」ではなぜ足りないのか

式 (6) が $\ge$ ではなく厳密な $>$ である点は、試験でもよく問われます。$f_s = 2f_{\max}$、つまり $f_{\max} = f_N$ ちょうどの場合を考えてみましょう。このとき原本の右端 $f_{\max}$ と $k=1$ の複製の左端 $f_s – f_{\max}$ がぴったり一致し、両者が $f = f_N$ という一点で接触します。接触点でスペクトルが足し合わされてしまうため、理想フィルタで切り出しても元の値には戻りません。

時間領域で見るともっと生々しくわかります。周波数 $f_N$ の正弦波を $f_s = 2f_N$ でサンプリングすると

$$ x[n] = \sin\!\left(2\pi f_N \frac{n}{f_s} + \varphi\right) = \sin(\pi n + \varphi) = (-1)^n \sin\varphi $$

となります。位相 $\varphi$ が 0 なら標本値はすべてゼロ、$\varphi = \pi/2$ なら $+1, -1, +1, \dots$ と最大振幅で振れます。つまり、同じ周波数・同じ振幅の信号でも、位相しだいで観測される振幅が 0 から最大まで変わってしまう。振幅情報が復元できない以上、これは「正しくサンプリングできている」とは言えません。だから条件は厳密な不等号なのです。

ナイキスト周波数ちょうどの正弦波を位相を変えて標本化したときの観測振幅の違い

3 つのパネルはいずれも振幅 1 の 500 Hz 正弦波を $f_s = 1000$ Hz で標本化したもので、違うのは位相だけです。左(位相 0)では標本点がすべてゼロ交差に当たり、標本列からは「無音」としか読めません。中央(位相 $\pi/4$)では振幅 0.71、右(位相 $\pi/2$)では振幅 1.0 に見えます。同じ信号なのに観測振幅が 0 から 1 まで変わるという事実が、$f_s = 2f_{\max}$ の等号が使えない理由をいちばん生々しく示しています。

実務では、この境界ぴったりを狙うことはまずありません。$f_{\max}$ と $f_N$ の間にガードバンドを設け、その帯域をアンチエイリアシングフィルタの遷移帯域に充てます。CD の 44.1 kHz は $f_N = 22.05\,\mathrm{kHz}$ で、可聴上限 20 kHz に対して 2.05 kHz のガードバンドを確保している、という設計です。

さて、条件が満たされていれば安全だとわかりました。では満たされなかったとき、具体的に何 Hz の偽物が現れるのか。ここが実務でいちばん役に立つ部分です。次の節で、折り返し後の周波数を与える公式を導きます。

折り返し周波数の導出 — 何 Hz に化けるのか

ナイキスト周波数を超えた成分は「消える」のではありません。別の周波数の本物そっくりの信号として、帯域内に化けて出てきます。この化けた先の周波数を エイリアス周波数(折り返し周波数)$f_a$ と呼びます。

導出はきわめて簡単で、しかも「1 行の三角関数の変形」で終わります。周波数 $f$、位相 $\varphi$ の実正弦波をサンプリングした標本列を書き下します。

$$ x[n] = \cos\!\left(2\pi f \frac{n}{f_s} + \varphi\right) $$

ここで、$f$ を「$f_s$ の整数倍」と「残り」に分解します。

$$ f = m f_s + r, \qquad m = \mathrm{round}\!\left(\frac{f}{f_s}\right) $$

$\mathrm{round}$ は四捨五入(最近接整数)なので、残り $r$ は必ず

$$ -\frac{f_s}{2} \le r \le \frac{f_s}{2}, \qquad \text{つまり} \quad -f_N \le r \le f_N $$

の範囲に収まります。この $f = mf_s + r$ を標本列に代入すると、位相の中身は

$$ 2\pi (m f_s + r)\frac{n}{f_s} + \varphi = 2\pi m n + 2\pi r \frac{n}{f_s} + \varphi $$

となります。ここが山場です。$m$ と $n$ はどちらも整数なので $2\pi m n$ は $2\pi$ の整数倍であり、余弦関数の周期性から、この項はまるごと捨てられます。残るのは

$$ \begin{equation} x[n] = \cos\!\left(2\pi r \frac{n}{f_s} + \varphi\right) \end{equation} $$

つまり、標本列だけを見ると、周波数 $f$ の信号は周波数 $r$ の信号と完全に区別がつきません。これがエイリアシングの正体です。

最後に $r$ の符号を処理します。$r \ge 0$ ならそのまま周波数 $r$、位相 $\varphi$ の正弦波です。$r < 0$ の場合は、余弦が偶関数($\cos(-\theta) = \cos\theta$)であることから

$$ \cos\!\left(-2\pi |r| \frac{n}{f_s} + \varphi\right) = \cos\!\left(2\pi |r| \frac{n}{f_s} – \varphi\right) $$

と書き直せます。周波数は $|r|$、ただし位相の符号が反転します。いずれの場合も観測される周波数は $|r|$ なので、まとめると

$$ \begin{equation} f_a = \left| f – \mathrm{round}\!\left(\frac{f}{f_s}\right) f_s \right| \end{equation} $$

が折り返し周波数の公式です。定義から $0 \le f_a \le f_N$ が常に成り立ちます。どんなに高い周波数を入力しても、観測される周波数は必ず 0 〜 $f_N$ の範囲に押し込められるわけです。

位相反転が意味すること — 車輪が逆に回る理由

$r < 0$ のときに位相が反転する、という副産物は、実は肉眼で確認できる現象です。複素指数 $e^{j2\pi f t}$ で同じ計算をすると

$$ x[n] = e^{j2\pi f n/f_s} = e^{j2\pi m n} e^{j2\pi r n /f_s} = e^{j2\pi r n/f_s} $$

となり、$r$ の符号がそのまま残ります。複素平面での回転方向は $r$ の符号で決まるので、$r < 0$ なら回転が逆向きに見えるのです。

数値を入れてみましょう。フレームレート $f_s = 24\,\mathrm{fps}$ のカメラで、スポーク 12 本のホイールを撮影します。ホイールが毎秒 $\nu$ 回転すると、スポークが同じ位置を通過する周波数は $f = 12\nu$ Hz です。

  • $\nu = 1.9\,\mathrm{rps}$ のとき $f = 22.8\,\mathrm{Hz}$。$f/f_s = 0.95$ なので $m = 1$、$r = 22.8 – 24 = -1.2\,\mathrm{Hz}$。負なので逆回転に見え、見かけの速さは $1.2\,\mathrm{Hz} \div 12 = 0.1\,\mathrm{rps}$ のゆっくりした回転になります。
  • $\nu = 2.0\,\mathrm{rps}$ のとき $f = 24\,\mathrm{Hz}$。$r = 0$ なので完全に静止して見えます。
  • $\nu = 2.1\,\mathrm{rps}$ のとき $f = 25.2\,\mathrm{Hz}$。$r = +1.2\,\mathrm{Hz}$ なので、今度は正方向にゆっくり回って見えます。

24fpsのカメラで車輪を撮ったときフレームごとのスポーク位置が逆回転・静止・正回転になる図

矢印はフレームごとのスポークの見かけ位置で、番号がフレーム番号です。左($r = -1.2$ Hz)ではフレームが進むほど矢印が時計回りに動いており、実際には前向きに回っている車輪が逆回転して見えます。中央($r = 0$)は 8 フレームすべてが同じ位置に重なり、車輪が止まって見える状態です。右($r = +1.2$ Hz)では反時計回りにゆっくり進みます。真ん中を境に見かけの回転方向が反転するのは、$r$ の符号が $+$ から $-$ に変わるからで、式 (8) の位相反転がそのまま目に見えているわけです。

車が加速していくと、車輪が「前回り → 停止 → 逆回り → 停止 → 前回り」を繰り返すのは、$r$ が $-f_N \to 0 \to +f_N$ を周期的に往復するからです。導いた公式が、日常の見え方をそのまま説明してしまうのは気持ちがいいところです。

実信号のスペクトルで見る折り返し

もう一つ、周波数領域の絵でも同じ結論を確認しておきます。実信号のスペクトルは $X(-f) = X^*(f)$ という共役対称性を持つため、周波数 $f$ の正弦波は $+f$ と $-f$ の 2 本のスペクトル線を持ちます。サンプリングによって、これらが $f_s$ 間隔で複製されます。

$+f$ 由来の線は $f – k f_s$ の位置に、$-f$ 由来の線は $-f + k f_s$ の位置に現れます。このうち $[0, f_N]$ に入るものだけが「観測される周波数」です。たとえば $f_s = 1000\,\mathrm{Hz}$、$f = 700\,\mathrm{Hz}$ のとき、$+700$ 由来の複製は $700, -300, 1700, \dots$、$-700$ 由来の複製は $-700, 300, 1300, \dots$ に並びます。$[0, 500]$ に入るのは $300\,\mathrm{Hz}$ だけ。式 (9) でも $\mathrm{round}(0.7) = 1$ より $f_a = |700 – 1000| = 300\,\mathrm{Hz}$ で、確かに一致します。

700Hzの正弦波を1000Hzで標本化したときのスペクトル線が300Hzに折り返す様子

青い線が $+700$ Hz 由来の複製($700, -300, 1700, \dots$)、赤い線が $-700$ Hz 由来の複製($-700, 300, 1300, \dots$)です。緑で塗った $[0, f_N]$ が実際に観測できる帯域で、その中に入っているのは赤い $300\,\mathrm{Hz}$ の線ただ 1 本だけであることが読み取れます。つまり「負の周波数側の複製がひとつ持ち上がって正の側に着地した」結果が、私たちの見る 300 Hz のピークの正体です。

「負の周波数側の複製が正の側に飛び込んでくる」というこの描像から、折り返し(folding)という呼び名が来ています。周波数軸を $f_N$ のところで鏡のように折り曲げると、高い周波数がぱたんと低い側に倒れ込む — その絵をイメージすると忘れません。

公式が手に入ったので、次は入力周波数を連続的に上げていったときに観測周波数がどう動くのかを追いかけてみましょう。単調に増えるわけでも、単に消えるわけでもない、独特の振る舞いが見えてきます。

観測周波数は三角波を描く

式 (9) を、入力周波数 $f$ の関数として眺めてみます。$f$ を 0 からじわじわ上げていくと、$f_a$ は次のように動きます。

  • $0 \le f \le f_N$:$m = 0$ なので $f_a = f$。観測周波数は入力どおりに増えます(正常領域)。
  • $f_N \le f \le f_s$:$m = 1$ なので $f_a = f_s – f$。入力を上げるほど観測周波数は下がります
  • $f_s \le f \le \tfrac{3}{2} f_s$:$m = 1$ なので $f_a = f – f_s$。再び増加に転じます。
  • $\tfrac{3}{2} f_s \le f \le 2f_s$:$m = 2$ なので $f_a = 2f_s – f$。また減少。

一般には、整数 $k$ に対して

$$ f_a(f) = \begin{cases} f – k f_s & \left(k f_s \le f \le \left(k + \tfrac{1}{2}\right) f_s\right) \\[4pt] (k+1) f_s – f & \left(\left(k + \tfrac{1}{2}\right) f_s \le f \le (k+1) f_s\right) \end{cases} $$

とまとめられます。つまり $f_a(f)$ は、周期 $f_s$、振幅 $0$ 〜 $f_N$ の三角波です。折り返しの「折り返し」という語感そのままのグラフになります。

観測周波数が入力周波数に対して三角波を描き同じ観測値を与える入力が複数存在することを示す図

赤い折れ線が式 (9) の $f_a(f)$ で、$f_N = 500$ Hz と 0 Hz のあいだを $f_s = 1000$ Hz 周期で往復しています。緑で重ねた区間($f < f_N$)だけが「入力どおりに記録できる」領域で、そこから右はすべて折り返し領域です。青い水平線(観測値 300 Hz)と折れ線の交点に注目すると、$300, 700, 1300, 1700, 2300, 2700\,\mathrm{Hz}$ の 6 本(さらに上も無限に)が同じ 300 Hz として観測されることがわかります。この多対一の構造こそが「情報が原理的に失われる」という表現の中身です。

この三角波の形から、実務上ただちに 2 つの帰結が読み取れます。

帰結 1:折り返し先は一意に決まるが、逆はそうでない。 ある入力周波数がどこに化けるかは式 (9) で一意に決まります。しかし観測された $f_a$ から元の $f$ を復元することはできません。$f_a = 300\,\mathrm{Hz}$($f_s = 1000\,\mathrm{Hz}$)を観測したとき、元の信号は $300, 700, 1300, 1700, 2300, \dots$ のどれでもありえます。情報が原理的に失われているというのは、この多対一の写像のことを指しています。

帰結 2:ナイキスト周波数のすぐ上が最も危険。 $f_N$ をわずかに超えた成分は $f_N$ のすぐ下に折り返してきます。信号帯域の上端付近に化けるので、「なんだか高域が汚い」という形で現れ、原因に気づきにくい。一方、$f_s$ に近い周波数の成分は DC 付近に折り返します。こちらは「なぜかオフセットが揺らぐ」という症状になります。どちらも、時間波形を眺めているだけでは絶対に見抜けません。

さらに三角波の描像は、どの帯域をフィルタで殺せばよいかも教えてくれます。信号帯域を $[0, f_p]$ とすると、$[0, f_p]$ に折り返してくるのは $f_a \le f_p$ となる $f$、すなわち三角波が低い値を取る区間です。$f_N$ より上でいちばん近いそのような区間は $[f_s – f_p,\ f_s + f_p]$ です。逆に言えば、$f_N$ から $f_s – f_p$ までの成分は、折り返しても $(f_p, f_N]$ というガードバンドの中にしか落ちないので、あとからディジタルフィルタで除去できます。この事実が、次節でアンチエイリアシングフィルタの遷移帯域を決める根拠になります。

具体例で計算してみる

抽象的な式だけでは身につかないので、数値を入れて確かめます。$f_s = 1000\,\mathrm{Hz}$($f_N = 500\,\mathrm{Hz}$)のサンプリングを考えます。

入力周波数 $f$ [Hz] $f/f_s$ $m = \mathrm{round}(f/f_s)$ $r = f – mf_s$ [Hz] 観測周波数 $f_a$ [Hz] 位相
300 0.30 0 +300 300 そのまま
500 0.50 0 or 1 ±500 500 振幅不定(境界)
700 0.70 1 −300 300 反転
1000 1.00 1 0 0(直流)
1300 1.30 1 +300 300 そのまま
1700 1.70 2 −300 300 反転
2200 2.20 2 +200 200 そのまま

$300, 700, 1300, 1700$ がすべて $300\,\mathrm{Hz}$ に化ける点に注目してください。標本列だけを見せられたら、この 4 つを区別する手立てはありません。$f = 1000\,\mathrm{Hz}$ が直流になるのも面白い例で、これは「1 周期にちょうど 1 点」しか取らないため、毎回まったく同じ位相の点を拾ってしまい、値が変化しないからです。

身近な数値例:オーディオ

$f_s = 44.1\,\mathrm{kHz}$(CD)で $f_N = 22.05\,\mathrm{kHz}$ です。仮にアンチエイリアシングフィルタが不十分で、$30\,\mathrm{kHz}$ の超音波成分が漏れ込んだとします。

$$ \frac{30}{44.1} = 0.680 \;\Rightarrow\; m = 1, \quad f_a = |30 – 44.1| = 14.1\,\mathrm{kHz} $$

$30\,\mathrm{kHz}$ は本来聞こえない超音波なのに、サンプリングを通ると $14.1\,\mathrm{kHz}$ という、はっきり聞こえる高音として現れます。しかも元の音楽とは何の関係もない周波数なので、金属的な異音として知覚されます。「聞こえない成分だから放っておいてよい」が通用しないのは、このためです。

身近な数値例:回転機械の振動計測

$f_s = 2\,\mathrm{kHz}$ で軸受の振動を測っているとします。$f_N = 1\,\mathrm{kHz}$ です。軸受の欠陥に起因する $2.4\,\mathrm{kHz}$ の高周波成分があると

$$ \frac{2.4}{2.0} = 1.2 \;\Rightarrow\; m = 1, \quad f_a = |2.4 – 2.0| = 0.4\,\mathrm{kHz} = 400\,\mathrm{Hz} $$

$400\,\mathrm{Hz}$ に現れます。もし回転数由来の成分がたまたま $400\,\mathrm{Hz}$ 付近にあれば、存在しない異常を検出したり、逆に本物の異常を見逃したりします。診断系では、サンプリング周波数を変えてスペクトルを取り直し、「ピークが動くかどうか」でエイリアスかどうかを判定する手法がよく使われます。式 (9) より、真の成分は $f_s$ を変えても位置が動かず、エイリアスは動くからです。

意図的にエイリアシングを使う — 帯域通過サンプリング

ここまでエイリアシングを「事故」として扱ってきましたが、条件を整えれば積極的に利用できます。信号が $[f_L, f_H]$ という帯域に限られていて、その中身だけに興味がある場合を考えます。帯域幅は $B = f_H – f_L$ です。

このとき、$f_s > 2f_H$ である必要はありません。$f_s > 2B$ を満たしつつ、複製同士が重ならないように $f_s$ を選べば、目的の帯域はきれいに $[0, f_N]$ の中へ折り返してきます。これを 帯域通過サンプリング(アンダーサンプリング、bandpass sampling)と呼びます。

たとえば $[100\,\mathrm{MHz},\ 102\,\mathrm{MHz}]$ の信号($B = 2\,\mathrm{MHz}$)を扱うのに、$204\,\mathrm{MHz}$ でサンプリングする必要はありません。$f_s = 20\,\mathrm{MHz}$ を選ぶと、$m = \mathrm{round}(100/20) = 5$ で $100\,\mathrm{MHz} \to 0\,\mathrm{MHz}$、$102\,\mathrm{MHz} \to 2\,\mathrm{MHz}$ に落ちてきます。ADC の速度を 1/10 にできるわけで、SDR やレーダー受信機で実際に使われている手法です。

ただし条件は厳しくなります。帯域外にわずかでも成分があれば、それも同じように折り返してくるからです。帯域通過サンプリングでは、ローパスではなくバンドパスのアンチエイリアシングフィルタが必須で、その阻止特性の要求は通常のローパスより厳しくなります。「エイリアシングを味方につける代わりに、フィルタ設計の負担が増える」というトレードオフです。

こうして見てくると、結局のところ勝負を決めているのはアンチエイリアシングフィルタだとわかります。では、そのフィルタにはどれだけの性能が必要なのでしょうか。次節で、必要な阻止量と次数を定量的に見積もります。

アンチエイリアシングフィルタの設計 — 必要な阻止量と次数

なぜ「あとから」では駄目なのか

まず最重要の事実を確認します。エイリアシングは A/D 変換の瞬間に起きる ので、ディジタル化された後で除去することは原理的に不可能です。折り返してきた $300\,\mathrm{Hz}$ の成分は、本物の $300\,\mathrm{Hz}$ の成分と数学的に完全に同一であり、両者を分離する情報がどこにも残っていません。

したがって、対策は A/D 変換の前に、アナログのローパスフィルタで $f_N$ 以上の成分を落としておく 以外にありません。これが アンチエイリアシングフィルタ(AAF)です。ADC のビット数を増やしても、サンプリング後にどんな高性能フィルタをかけても、この事故は防げません。「AAF は ADC の前段に必ず要る」というのが、計測系設計の鉄則です。

必要な阻止量は ADC の分解能で決まる

では、どれだけ落とせばよいのでしょうか。基準になるのは ADC の量子化雑音レベル です。$N$ ビットの理想 ADC のフルスケール正弦波に対する SN 比は、よく知られた式

$$ \begin{equation} \mathrm{SNR}_{\mathrm{dB}} = 6.02 N + 1.76 \;[\mathrm{dB}] \end{equation} $$

で与えられます。折り返してくる成分が、この量子化雑音より十分小さければ、実質的に「見えない」ことになります。したがって AAF に要求される阻止量 $A_{\mathrm{stop}}$ は

$$ A_{\mathrm{stop}} \ge 6.02 N + 1.76 \;[\mathrm{dB}] $$

が目安です。具体的には

  • 12 ビット ADC:$6.02 \times 12 + 1.76 = 74.0\,\mathrm{dB}$
  • 16 ビット ADC:$6.02 \times 16 + 1.76 = 98.1\,\mathrm{dB}$
  • 24 ビット ADC:$6.02 \times 24 + 1.76 = 146.2\,\mathrm{dB}$

となります。16 ビットで約 98 dB、つまり振幅比で約 8 万分の 1 まで落とせというかなり過酷な要求です。

遷移帯域はどこからどこまでか

前節の三角波の議論から、阻止帯域の開始周波数は $f_N$ ではなく

$$ \begin{equation} f_{\mathrm{stop}} = f_s – f_p \end{equation} $$

です($f_p$ は信号帯域の上端 = 通過端)。$f_N < f < f_s - f_p$ の成分は折り返しても $(f_p, f_N]$ というガードバンドの中にしか落ちないので、A/D 変換後にディジタルフィルタで安全に除去できるからです。この「$f_N$ ではなく $f_s - f_p$ でよい」という一歩が、実際の設計をかなり楽にします。

したがって AAF の仕様は次の 3 点にまとまります。

  1. 通過端 $f_p$ までは平坦(リプル許容値以内)
  2. 阻止端 $f_{\mathrm{stop}} = f_s – f_p$ 以上で $A_{\mathrm{stop}}$ dB 以上の減衰
  3. その間の遷移帯域幅は $f_{\mathrm{stop}} – f_p = f_s – 2f_p$

3 番目の式が本質的です。遷移帯域幅は $f_s – 2f_p$。$f_s$ を $2f_p$ ぎりぎりに設定すると遷移帯域幅がゼロに近づき、フィルタの要求が発散します。逆に $f_s$ を大きく取れば、遷移帯域はいくらでも広くできます。

アンチエイリアシングフィルタの通過帯域・ガードバンド・折り返し許容域・阻止帯域の関係図

$f_s = 48\,\mathrm{kHz}$、$f_p = 20\,\mathrm{kHz}$ の場合の帯域の割り付けです。緑が信号を通す通過帯域、橙が $f_p$ と $f_N = 24$ kHz のあいだのガードバンド、紫が $f_N$ から $f_{\mathrm{stop}} = f_s – f_p = 28$ kHz までの「折り返してもガードバンドにしか落ちない」領域、赤が本当に落とさなければならない阻止帯域です。阻止端が $f_N$ ではなく 28 kHz でよいおかげで、遷移帯域幅が 4 kHz ではなく $f_s – 2f_p = 8\,\mathrm{kHz}$ と倍に取れている点が読み取りどころです。青い曲線は必要阻止量 98 dB を阻止端でちょうど満たす 34 次バターワースで、この条件でもまだかなり急峻なフィルタが要ることがわかります。

バターワースフィルタの必要次数

具体的な次数を見積もるために、$n$ 次バターワースローパスフィルタ(遮断周波数 $f_c$)の振幅特性

$$ |H(f)|^2 = \frac{1}{1 + (f/f_c)^{2n}} $$

を使います。減衰量(dB)は

$$ A(f) = -10\log_{10}|H(f)|^2 = 10\log_{10}\!\left[1 + \left(\frac{f}{f_c}\right)^{2n}\right] $$

阻止帯域では $(f/f_c)^{2n} \gg 1$ なので 1 は無視でき、対数の中身が $\left(f/f_c\right)^{2n}$ になります。対数の指数を前に出すと

$$ A(f) \approx 10 \cdot 2n \log_{10}\frac{f}{f_c} = 20 n \log_{10}\frac{f}{f_c} \;[\mathrm{dB}] $$

つまり、1 次あたり $20\log_{10}(f/f_c)$ dB ずつ稼げるわけです。これを $f = f_{\mathrm{stop}}$、$f_c = f_p$ として $A_{\mathrm{stop}}$ 以上にする条件を $n$ について解くと

$$ \begin{equation} n \ge \frac{A_{\mathrm{stop}}}{20 \log_{10}\!\left(\dfrac{f_s – f_p}{f_p}\right)} \end{equation} $$

分母の $\log$ の中身が、遷移帯域の「比」です。この比が 1 に近い(= $f_s$ が $2f_p$ に近い)と分母が 0 に近づき、必要次数が爆発します。

数値で殴られてみる — CD の 44.1 kHz

$f_p = 20\,\mathrm{kHz}$、16 ビット相当で $A_{\mathrm{stop}} = 98\,\mathrm{dB}$ として、$f_s$ を変えながら式 (12) を計算します。

$f_s$ [kHz] $f_{\mathrm{stop}} = f_s – f_p$ [kHz] 比 $f_{\mathrm{stop}}/f_p$ 1 次あたり [dB] 必要次数 $n$
44.1 24.1 1.205 1.62 61
48 28 1.40 2.92 34
96 76 3.80 11.60 9
192 172 8.60 18.69 6
384 364 18.20 25.20 4

数字の落差に驚いてください。$f_s = 44.1\,\mathrm{kHz}$ で 16 ビット相当の阻止量を素直に得ようとすると 61 次のアナログフィルタが必要になります。アナログの世界で 61 次は、部品数・温度ドリフト・位相特性のどれを取っても現実的ではありません。初期の CD プレーヤーが「ブリックウォールフィルタ」と呼ばれる急峻なアナログフィルタを積み、その位相歪みが音質論争を呼んだのは、まさにこの表の 1 行目に立たされていたからです。

一方 $f_s = 192\,\mathrm{kHz}$ なら 6 次で足ります。6 次バターワースなら 2 次のセクションを 3 段重ねるだけで、オペアンプ数個で作れます。これが オーバーサンプリング の威力です。

オーバーサンプリングとデシメーション

現代の ADC がほぼ例外なく採っている戦略は次の通りです。

  1. 目標よりずっと高い $f_s$ でサンプリングする(例:目標 48 kHz に対して 6.144 MHz、128 倍オーバーサンプリング)
  2. アナログ AAF は緩い次数で済ませる(遷移帯域が $f_s – 2f_p$ と広大なので、2〜3 次で十分)
  3. ディジタル領域で急峻な FIR ローパスをかける(ディジタルなら 100 次以上でも平気で、しかも正確な線形位相にできる)
  4. デシメーション(間引き)で目標のサンプリングレートに落とす

緩いアナログフィルタと高速ADC・急峻なディジタルFIR・デシメーションを組み合わせた信号処理の流れ図

流れ図の左半分(アナログ側)に置かれているのは「緩い AAF」だけで、急峻なフィルタは右半分(ディジタル側)に移されています。6.144 MHz でサンプリングすれば $f_N$ は 3.072 MHz なので、20 kHz の信号帯域に対して遷移帯域が桁違いに広く、2〜3 次のアナログフィルタで足ります。難所を「部品点数と温度ドリフトに苦しむアナログ」から「係数を並べるだけのディジタル」へ付け替えた、というのがこの図の要点です。

このとき、ステップ 3 のディジタルフィルタは「ステップ 4 で間引くときのアンチエイリアシングフィルタ」として働きます。間引きは実質的に低いレートでのサンプリングなので、間引き後の $f_N$ を超える成分をあらかじめ落としておく必要があるわけです。アナログの難しい仕事を、ディジタルの得意な仕事に付け替える — これがオーバーサンプリングの本質です。

$\Delta\Sigma$(デルタシグマ)ADC はこの考えをさらに推し進め、量子化雑音を高域へ追いやる ノイズシェーピング と組み合わせることで、1 ビットの粗い量子化器から 24 ビット級の分解能を作り出します。ディジタルフィルタの設計そのものはディジタルフィルタの基礎(FIR/IIR)を解説で扱っていますので、あわせて読むと設計の全体像がつながります。

理論はここまでです。あとは実際に手を動かして、折り返しが起きる瞬間を自分の目で見るのがいちばん早い。次の節で、Python で 3 つの実験をします。

Python での実装と可視化

実験 1:高い正弦波が低い正弦波に化ける

まず、式 (9) を最も直接的に確かめます。$f_s = 1000\,\mathrm{Hz}$ で $f = 1300\,\mathrm{Hz}$ の正弦波をサンプリングすると、$f_a = 300\,\mathrm{Hz}$ に化けるはずです。標本点が「本物の 300 Hz 正弦波」の上にぴたりと乗るかどうかを見ます。

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

# 日本語フォント設定
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

fs = 1000.0        # サンプリング周波数 [Hz]
f_in = 1300.0      # 入力信号の周波数 [Hz]

# 折り返し周波数の公式 f_a = |f - round(f/fs)*fs|
def alias_freq(f, fs):
    return np.abs(f - np.round(f / fs) * fs)

f_a = alias_freq(f_in, fs)
print(f"入力 {f_in} Hz -> 観測 {f_a} Hz (fs={fs} Hz, fN={fs/2} Hz)")

# 連続波形の近似(十分高いレートで描く)
t_cont = np.linspace(0, 0.01, 5000)
x_cont = np.cos(2 * np.pi * f_in * t_cont)
x_alias = np.cos(2 * np.pi * f_a * t_cont)

# 実際の標本点
n = np.arange(0, int(0.01 * fs) + 1)
t_samp = n / fs
x_samp = np.cos(2 * np.pi * f_in * t_samp)

plt.figure(figsize=(11, 4.5))
plt.plot(t_cont * 1000, x_cont, color="gray", lw=1.0,
         label=f"本物の入力 {f_in:.0f} Hz(連続波形)")
plt.plot(t_cont * 1000, x_alias, color="crimson", lw=2.0, ls="--",
         label=f"折り返し先 {f_a:.0f} Hz")
plt.plot(t_samp * 1000, x_samp, "o", color="navy", ms=9,
         label=f"標本点({fs:.0f} Hz でサンプリング)")
plt.xlabel("時間 [ms]")
plt.ylabel("振幅")
plt.title("ナイキスト周波数を超えた正弦波は低い周波数に化ける")
plt.legend(loc="upper right", fontsize=9)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

1300Hzの正弦波を1000Hzで標本化した点が300Hzの正弦波の上にも乗ることを示すグラフ

このグラフから、エイリアシングの本質が一目でわかります。標本点(青丸)は、灰色の $1300\,\mathrm{Hz}$ の波形の上にあると同時に、赤い破線の $300\,\mathrm{Hz}$ の波形の上にも完全に乗っています。標本点だけを渡された受信側には、元が 1300 Hz だったのか 300 Hz だったのかを判断する材料が一切ありません。「情報が失われる」という抽象的な表現が、この図では「2 本の波が同じ点を通る」という具体的な絵になっています。

実験 2:観測周波数の三角波を実測で確かめる

次に、入力周波数を 0 から $3f_s$ まで掃引しながら、実際に FFT でピーク周波数を測り、式 (9) の予測と一致するかを確かめます。式を信じるのではなく、実測で検証するのが大事なところです。

import numpy as np
import matplotlib.pyplot as plt

fs = 1000.0
N = 1000                      # 1 秒分 → 周波数分解能 1 Hz
n = np.arange(N)

f_sweep = np.arange(0, 3001, 5, dtype=float)   # 入力周波数を掃引
f_measured = np.zeros_like(f_sweep)

for i, f in enumerate(f_sweep):
    x = np.cos(2 * np.pi * f * n / fs)         # サンプリングされた標本列
    X = np.abs(np.fft.rfft(x))                 # 片側スペクトル
    freqs = np.fft.rfftfreq(N, d=1 / fs)
    f_measured[i] = freqs[np.argmax(X)]        # ピーク周波数 = 観測周波数

f_theory = np.abs(f_sweep - np.round(f_sweep / fs) * fs)
err = np.max(np.abs(f_measured - f_theory))
print(f"理論式と実測ピークの最大誤差: {err:.3f} Hz")
plt.figure(figsize=(11, 4.5))
plt.plot(f_sweep, f_theory, color="crimson", lw=2.2,
         label=r"理論値 $f_a=|f-\mathrm{round}(f/f_s)f_s|$")
plt.plot(f_sweep, f_measured, "o", color="navy", ms=3.5, alpha=0.7,
         label="FFT で実測したピーク周波数")
plt.plot([0, fs / 2], [0, fs / 2], color="green", lw=3, alpha=0.35,
         label="正しく記録できる領域")
for k in range(1, 4):
    plt.axvline(k * fs, color="gray", ls=":", lw=1)
plt.axhline(fs / 2, color="orange", ls="--", lw=1.2,
            label=f"ナイキスト周波数 {fs/2:.0f} Hz")
plt.xlabel("入力信号の周波数 $f$ [Hz]")
plt.ylabel("観測される周波数 $f_a$ [Hz]")
plt.title("入力周波数を上げていくと観測周波数は三角波状に折り返す")
plt.legend(fontsize=9, loc="upper right")
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

FFTで実測したピーク周波数が理論式の三角波にぴたりと重なることを示すグラフ

実測点(青)が理論曲線(赤)にぴたりと重なり、最大誤差はゼロ(0.000 Hz)になります。読み取りどころは 3 つあります。第一に、$f < f_N = 500\,\mathrm{Hz}$ の領域(緑の直線)でだけ、観測周波数が入力周波数に一致します。第二に、$f_N$ を超えると観測周波数は下がり始め、$f = f_s = 1000\,\mathrm{Hz}$ で 0 Hz(直流)に達します。第三に、そこからまた上昇して $1500\,\mathrm{Hz}$ で $f_N$ に達し、以後 $f_s$ 周期で同じ三角波を繰り返します。同じ観測周波数を与える入力周波数が無数にあることが、グラフの水平線を引いてみればすぐ確認できます。

実験 3:チャープ信号のスペクトログラムで折り返しを見る

単一周波数だけでは折り返しの「動き」が見えません。周波数が時間とともに直線的に上がっていくチャープ信号を使うと、スペクトログラム上で折り返しが起きる瞬間が動画のように観察できます。真値(高レートでサンプリング)と、低レートでサンプリングしたものを並べます。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

T = 2.0            # 信号長 [s]
f0, f1 = 0.0, 1200.0   # チャープの開始・終了周波数 [Hz]
fs_hi = 8000.0     # 参照用の高いサンプリング周波数
fs_lo = 1000.0     # 折り返しを起こすサンプリング周波数(fN = 500 Hz)

t_hi = np.arange(0, T, 1 / fs_hi)
t_lo = np.arange(0, T, 1 / fs_lo)
x_hi = signal.chirp(t_hi, f0=f0, t1=T, f1=f1, method="linear")
x_lo = signal.chirp(t_lo, f0=f0, t1=T, f1=f1, method="linear")

fig, axes = plt.subplots(1, 2, figsize=(13, 4.8))
for ax, (x, fsr, ttl) in zip(axes, [
        (x_hi, fs_hi, f"参照:$f_s$={fs_hi:.0f} Hz(折り返しなし)"),
        (x_lo, fs_lo, f"問題:$f_s$={fs_lo:.0f} Hz($f_N$=500 Hz)")]):
    f_sp, t_sp, S = signal.spectrogram(x, fs=fsr, nperseg=256,
                                       noverlap=224, scaling="spectrum")
    ax.pcolormesh(t_sp, f_sp, 10 * np.log10(S + 1e-12),
                  shading="gouraud", cmap="magma", vmin=-70, vmax=0)
    ax.set_xlabel("時間 [s]")
    ax.set_ylabel("周波数 [Hz]")
    ax.set_title(ttl)
    ax.set_ylim(0, 1300 if fsr > 2000 else fsr / 2)
axes[1].axhline(fs_lo / 2, color="cyan", ls="--", lw=1.2)
axes[1].text(0.05, fs_lo / 2 - 45, "ナイキスト周波数", color="cyan", fontsize=9)
plt.tight_layout()
plt.show()

高いサンプリング周波数と低いサンプリング周波数でチャープ信号のスペクトログラムを比べた図

左(参照)のスペクトログラムでは、明るい線が左下から右上へまっすぐ伸びています。0 Hz から 1200 Hz まで直線的に上がるチャープの、素直な姿です。ところが右($f_s = 1000\,\mathrm{Hz}$)では、線が $t \approx 0.83\,\mathrm{s}$(瞬時周波数が 500 Hz に達する時刻)でナイキスト周波数にぶつかり、そこで折れて下向きに転じます。さらに $t \approx 1.67\,\mathrm{s}$(瞬時周波数 1000 Hz)で 0 Hz に達し、再び上向きに折れ返って、$t = 2\,\mathrm{s}$ では 200 Hz に到達します。式 (9) の予測($f = 1200 \to f_a = |1200-1000| = 200\,\mathrm{Hz}$)とぴったり一致しています。

この図が示す実務上の教訓は重い。もし右の図しか手元になければ、「途中で周波数が下がる不思議な信号」というまったく間違った解釈をしてしまいます。スペクトログラムで折れ線が $f_N$ や 0 Hz で反射しているように見えたら、まずエイリアシングを疑うべきです。

実験 4:アンチエイリアシングフィルタの必要次数

最後に、式 (12) が主張する「オーバーサンプリングでフィルタが劇的に楽になる」を、バターワース特性の実カーブで確かめます。

import numpy as np
import matplotlib.pyplot as plt

f_p = 20e3            # 信号帯域の上端(通過端)[Hz]
A_stop = 6.02 * 16 + 1.76   # 16bit ADC 相当の必要阻止量 [dB]
print(f"必要阻止量 A_stop = {A_stop:.2f} dB")

fs_list = [44.1e3, 48e3, 96e3, 192e3, 384e3]
rows = []
for fs in fs_list:
    f_stop = fs - f_p                      # 阻止帯域の開始周波数
    db_per_order = 20 * np.log10(f_stop / f_p)
    n_req = int(np.ceil(A_stop / db_per_order))
    rows.append((fs, f_stop, f_stop / f_p, db_per_order, n_req))
    print(f"fs={fs/1e3:6.1f} kHz | f_stop={f_stop/1e3:6.1f} kHz | "
          f"比={f_stop/f_p:5.2f} | {db_per_order:5.2f} dB/次 | n>={n_req}")
# バターワース特性の実カーブで確認
f = np.logspace(np.log10(1e3), np.log10(1e6), 2000)

plt.figure(figsize=(11, 5.2))
colors = ["#d62728", "#ff7f0e", "#2ca02c", "#1f77b4", "#9467bd"]
for (fs, f_stop, ratio, dbo, n_req), c in zip(rows, colors):
    A = 10 * np.log10(1 + (f / f_p) ** (2 * n_req))     # n_req 次バターワース
    plt.semilogx(f, -A, color=c, lw=1.8,
                 label=f"$f_s$={fs/1e3:.1f} kHz → {n_req} 次")
    plt.plot([f_stop], [-A_stop], "o", color=c, ms=8)
plt.axhline(-A_stop, color="k", ls="--", lw=1.2)
plt.text(1.2e3, -A_stop + 4, f"必要阻止量 {A_stop:.0f} dB(16bit相当)", fontsize=9)
plt.axvline(f_p, color="gray", ls=":", lw=1.2)
plt.text(f_p * 1.05, -140, "通過端 20 kHz", fontsize=9, color="gray")
plt.xlabel("周波数 [Hz]")
plt.ylabel("振幅応答 [dB]")
plt.title("必要次数はオーバーサンプリングで劇的に下がる(丸=阻止端で必要量を満たす点)")
plt.ylim(-160, 5)
plt.legend(fontsize=9, loc="lower left")
plt.grid(alpha=0.3, which="both")
plt.tight_layout()
plt.show()

サンプリング周波数を上げるとバターワースフィルタの必要次数が61次から4次まで下がることを示すグラフ

出力される表と図を突き合わせると、設計の勘所が見えてきます。$f_s = 44.1\,\mathrm{kHz}$ では阻止端が $24.1\,\mathrm{kHz}$ と通過端 $20\,\mathrm{kHz}$ のすぐ上にあり、1 次あたりわずか 1.6 dB しか稼げないため 61 次が必要です。$f_s$ を 192 kHz に上げると阻止端は 172 kHz まで遠ざかり、1 次あたり 18.7 dB も稼げるので 6 次で足ります。必要次数がほぼ 1 桁下がるわけです。

図の各曲線が、対応する丸印(阻止端)でちょうど破線(必要阻止量)に触れているのも確認してください。次数を増やすと曲線の傾きが急になり、$f_s$ を上げると丸印が右へ移動する。この 2 つのつまみのどちらを回すかがトレードオフで、アナログ部品を減らしたいならオーバーサンプリング(右へ)、サンプリングレートを上げられないなら高次フィルタ(急峻に)を選ぶ、という構図です。現代の設計がほぼ例外なく前者を選ぶのは、ディジタル側の処理コストが年々安くなっているからにほかなりません。

よくある誤解を整理する

最後に、現場で繰り返し見かける誤解を、これまで導いた式に照らして片付けておきます。

誤解 1:「$f_s = 2f_{\max}$ でちょうど足りる」 式 (6) は厳密な不等号です。等号では複製が 1 点で接触し、しかも $f_N$ ちょうどの成分は位相しだいで振幅が 0 になり得ます。実務では $f_s \ge 2.2\,f_{\max}$ 程度、オーディオのように急峻なフィルタを避けたい用途では $f_s \ge 2.5\,f_{\max}$ 以上を取ります。

誤解 2:「ビット数を増やせば折り返しも減る」 量子化雑音(振幅方向の誤差)とエイリアシング(周波数軸の巻き取り)は、まったく別の現象です。むしろビット数を増やすほど量子化雑音の床が下がるので、式 (10) が示すとおりAAF に要求される阻止量は増えます。24 ビット ADC を採用したなら、AAF も 146 dB 級の要求に耐えるものが必要になります。

誤解 3:「サンプリング後にディジタルフィルタをかければ直せる」 直せません。折り返してきた成分は帯域内で本物と完全に重なっており、両者を分離する情報が失われています。ディジタルフィルタで除去できるのは、$(f_p, f_N]$ のガードバンドに落ちた成分だけです。

誤解 4:「エイリアシングは音や振動だけの話」 画像でも起きます。2 次元の空間サンプリングでは、細かい縞模様(高い空間周波数)がセンサの画素ピッチのナイキスト空間周波数を超えると、モアレという低い空間周波数の縞に化けます。多くのカメラが撮像素子の前に光学ローパスフィルタ(OLPF)を置いているのは、まさにこの光学版アンチエイリアシングフィルタです。CG のレンダリングで「アンチエイリアシング」と呼ばれる処理も、原理は同じで、サンプリング前にぼかす(帯域制限する)かスーパーサンプリングする($f_s$ を上げる)かの二択です。

誤解 5:「観測されたスペクトルのピーク位置から、元の周波数がわかる」 式 (9) は多対一の写像なので、逆は一意に定まりません。$f_a$ を観測しただけでは、元が $f_a$ なのか $f_s – f_a$ なのか $f_s + f_a$ なのか判別できません。判別したければ、サンプリング周波数を変えて測り直し、ピークが動くかどうかを見るしかありません。

まとめ

本記事では、ナイキスト周波数とエイリアシングを、周波数領域の一枚の絵から導出しました。

  • ナイキスト周波数 $f_N = f_s/2$ は、サンプリング周波数 $f_s$ で正しく記録できる信号周波数の上限。信号側の量であるナイキストレート $2f_{\max}$ とは別物
  • サンプリングは インパルス列との積 として書け、インパルス列のフーリエ変換が $f_s\sum_k \delta(f-kf_s)$ になることから、$X_s(f) = f_s\sum_k X(f – kf_s)$ というスペクトルの周期化が導かれる
  • 複製が重ならない条件 $f_{\max} < f_s - f_{\max}$ から $f_s > 2f_{\max}$、すなわち $f_{\max} < f_N$ が出る。等号ではスペクトルが接触し、$f_N$ ちょうどの成分は位相しだいで振幅が消える
  • 折り返し周波数は $f_a = |f – \mathrm{round}(f/f_s) f_s|$。$f = mf_s + r$ と分解したとき $2\pi mn$ が $2\pi$ の整数倍で消えることが理由で、$r<0$ のときは位相が反転する(車輪が逆回転して見える正体)
  • 入力周波数を上げていくと観測周波数は 周期 $f_s$、振幅 $0$〜$f_N$ の三角波 を描く。写像は多対一なので、折り返した成分から元の周波数は復元できない
  • エイリアシングは A/D 変換の瞬間に起きるため、ディジタル処理では除去不能。対策は前段のアナログ アンチエイリアシングフィルタ 一択
  • 必要阻止量は ADC のビット数から $A_{\mathrm{stop}} \approx 6.02N + 1.76$ dB、阻止端は $f_s – f_p$。バターワースの必要次数は $n \ge A_{\mathrm{stop}} / \left[20\log_{10}\left((f_s-f_p)/f_p\right)\right]$ で見積もれる
  • $f_s$ を $2f_p$ ぎりぎりに取ると必要次数が爆発する(20 kHz / 44.1 kHz で 61 次)。オーバーサンプリング+ディジタルフィルタ+デシメーション が現代の標準解

ナイキスト周波数は「$f_s$ を 2 で割るだけ」の公式ですが、その背後には「スペクトルが $f_s$ 間隔で複製される」という一枚の絵があります。この絵さえ思い浮かべられれば、折り返し先の計算も、フィルタの遷移帯域の決め方も、帯域通過サンプリングの設計も、その場で導き直せます。公式ではなく絵で覚えてください。

次のステップとして、以下の記事も参考にしてください。