デジタル信号を正しく再生するには? ~サンプリング定理の意味

Twitterで以下のように発言をしたら、案の定、あまり理解されなかったので頑張って説明してみます。

サンプリング定理は「信号の最大周波数」の2倍より早い速度でサンプリングすれば元信号の情報は完全に再現できる(一意に決まる)、とは言ってるけども、サンプリングしたデータをそのまま再生したとき元波形が再現できるとは一言も言ってない。

https://twitter.com/nabe_abk/status/777874934424940544

なおこの記事では、フーリエ変換が存在する連続信号(関数)のみを考えることにします。

フーリエ変換とは

フーリエ変換とは、無限に続く任意の信号を、無限の周波数までのsin波とcos波の重ねあわせとして表す変換のことです。

f(t) = 2\pi\int_{-\infty}^{\infty}F(2\pi f) (\cos(2\pi ft)+j\sin(2\pi ft) ) df

この式が意味していることは、「信号f(t)F(2\pi f)の周波数fについての積分」によって得られるということで、より工学的に解釈すればF(\omega)は周波数ごとの信号の強さ(位相差)を表していると言えます。

ですので、F(\omega)のことを周波数スペクトルと言うことがあります。

より詳しい解説は専門文献を当たってください。また「FFTとは? ~本当は正しくないFFTの周波数特性~」も参考に。

式の導出

※苦手な方は読み飛ばしてください。

f(t)のフーリエ変換F(\omega)を数式で書くと以下のようになります。

F(\omega) = \frac{1}{2\pi}\int_{-\infty}^{\infty}f(t) \mathrm{e}^{-j\omega t} dt

このフーリエ変換F(\omega)を用いて、元のf(t)を表すと次のようになります。

\begin{eqnarray} f(t) &=& \int_{-\infty}^{\infty}F(\omega) \mathrm{e}^{j\omega t} d\omega \nonumber \\ &=& \int_{-\infty}^{\infty}F(\omega) ( \cos(\omega t) + j\sin(\omega t) ) d\omega \nonumber \end{eqnarray}

ここで、\omega=2\pi fと置けば最初に書いた定義式になります。

サンプリング定理とは?

いくつか定義の仕方がありますが、ここでは工学的に少し分かりやすく定義することにします。「以上」「以下」「超えた」「未満」の表現に注意してください。*1

周波数f_mHz以下の成分しか含まない信号f(t)のフーリエ変換をF(\omega)としたとき、2f_mHzを超えた速度でサンプリングした信号g(t)のフーリエ変換G(\omega)は、f_mHz以下の範囲でF(\omega)と一致する。

この意味を理解するため、例を示して説明します。

sampling_01.png

青い信号がf(t)、緑の信号がそれをサンプリングした信号g(t)です。

f(t)g(t)はサンプリング定理の条件である「2倍を超えた周波数でサンプリングする」を満たしていますので、これらのフーリエ変換であるF(\omega)G(\omega)はサンプリング周波数の半分である2.25Hzより小さい周波数において一致します。

F(2\pi f) = G(2\pi f) \qquad (f < 2.25 = \frac{4.5}{2})

ちなみに、f(t)f_mHz(例なら2.25Hz)以上の信号を含む場合、エイリアシング(折り返し歪)と呼ばれる現象が起こり、フーリエ変換の結果は一致しなくなります。

*1 : 例えば「最大周波数の2倍以上でサンプリングすれば良い」は明確に誤りです。正しくは「最大周波数の2倍を超えた速度でサンプリングすれば良い」になります。同じく「サンプリング周波数の1/2以下の周波数しか含まなければ良い」は誤りで、「1/2未満周波数」が正解です。

サンプリング定理が示すもの

  • 元信号:f(t)
  • サンプリング信号:g(t)
  • サンプリング周波数:f_s
  • それぞれのフーリエ変換:F(\omega), G(\omega)
  • f(t)f_s/2未満の信号しか含まない

このとき、F(\omega)G(\omega)\displaystyle \frac{f_s}{2}Hz未満の領域について一致する。

さらに次のことが言えます。

  • f_s/2を超えた部分で、F(\omega)G(\omega)は一致しない。
    • f(t)f_s/2を超えた周波数成分を含まない。
    • g(t)f_s/2を超えた周波数成分を含む。
  • サンプリング定理が成り立つのは、無限に続く連続信号においてである。

サンプリングされた信号とその再生

工学書の多くで次に来るのは「サンプリング定理より信号の最大周波数の2倍より高い(速い)周波数でサンプリングすれば良い」という解説であり、「その条件を満たす信号を再生し理想LPF(ローパスフィルタ)に通せば元の信号に戻る」という解説です。

サンプリング定理の前提条件を満たす限り実際にそれは正しいのですが、それと同時にこれを読んだ人の多くが誤解をする内容にもなっています。それを理解するために、サンプリングされたデータを実際にアナログ信号として再生することを考えます。

一般的なデータ再生

sampling_02.png

サンプリングされたデータをそのまま再生する場合は*2、通常赤い線のようになります。これにLPFをかけるとどうなるでしょうか?

sampling_03.png

手書きで申し訳ないですが、緑の線のようになります。

フィルタの次数を上げて理想LPFに近づけろと言う方が出てくると思うので、サンプリング周波数を1kHzにアップサンプリングした上で1023点のFIR型LPFで処理した例を次に示します。

sampling_fir_1023.png

まだ足りないという方のために4095点のFIR型LPFの結果も示します。

sampling_fir_4095.png

形は近くなりましたが遅延が酷いですね。*3

なぜこのようなことになるかと言うと、理想LPFなんて作れないからです。それも「現実だと誤差あるし、しょうがないね」というレベルではなく、近似としても理想LPFを実現することはかなり困難です。

*2 : このような再生方法を0次ホールドと言います。

*3 : 44.1kHzのオーディオデータだとしたら、4095点FIRでは信号が50ms程度遅延する計算になります。現実的ではないと言える遅延でしょう。

サンプリング定理が成り立つ条件

  1. 理想LPFが存在すること。
  2. ある一点のデータを再生するとき、その前後のデータが無限に存在すること。
  3. インパルスを再生できること。

(1)だけで理解している気になっている人がとても多い。理想LPFってどうやって作るか知ってますか?

「ある周波数f_cHzでスパンと切れる周波数特性を逆フーリエ変換したフィルタを作ればいい」と思った人、正しいです。でもその関数がどんな形をしているかちゃんと知っていますか?

その関数はsinc関数と呼ばれる形になります。

sampling_sinc.png

この関数、見て分かるように正負方向に無限に続いています。理想LPFは、サンプリングした各点のデータとsinc関数を畳み込み演算することで実現できます。

つまりここから導き出される条件が(2)です。過去と未来の無限のデータを持ってきて処理するなんて不可能です。*4

また「データを再生して理想LPFに通す」という教科書の解説をよく考えずに鵜呑みしていると、(3)サンプリングしたデータの再生はインパルスで行うことを見落とします。サンプリング定理の項目で教科書に書かれるデータの再生とは「インパルスでのデータ再生」を示します。インパルスとは幅0で面積1のパルスのことで、これを生成することも不可能です。

インパルスで再生してからLPFで処理すると

そのまま再生(0次ホールド)ではなく、疑似インパルスで再生し、4095点FIR LPF処理した場合の図を掲載しておきます。

sampling_fir_4095_impulse.png

これでもまだ誤差がありますね。

*4 : 有限長の音声データ、音楽1曲とかなら頑張れば実現できるだろうとツッコミを入れる方がおられるかもしれません。現実の再生装置はリアルタイム処理をしなければならないので、様々な制約により、どんなに頑張っても1000点ぐらいのデータ処理が限界です。力技でデータ処理を可能にしたとしても、ある有限区間にも無限個のデータが存在しない限り理想LPFにはなりません。

デジタル信号を正しく再生する方法

前節で述べたとおり、すべての点をインパルスで再生し、それらについてsinc関数を畳み込み(フィルタする)必要があるので、まとめると以下のようになります。

すべてのデータ1点1点に対して、その記録してある値がピークとなるようなsinc関数を生成し、無限の時間と分解能でそのすべての点のsinc関数を加算する。

正確を期すため数式でも示しておきますが、苦手な人は読み飛ばしてください。

  • 元の信号 : f(t)
  • サンプリングデータ点数 : N
  • サンプリングデータ列 : g(n) \quad (n=0 \sim N-1)
  • サンプリング周波数 : f_s Hz
  • サンプリング周期 : T = 1/f_s sec
f(t) = \sum_{k=0}^{N-1} g(k)\frac{\sin(\pi(t/T-k))}{\pi(t-kT)}

sinc関数\mathrm{sinc}(x)=\sin(\pi x)/\pi xを使って書くとこんな感じ。

f(t) = T \sum_{k=0}^{N-1} g(k)\mathrm{sinc}\frac{t-kT}{T}

サンプリング周期Tずつずらしながら、すべでの点についてsinc関数を重ね合わせるイメージです。

例。T=0.25, N=4, g(n)=1, 3, 0, -2のとき。

f(t) = 0.25 \left( \mathrm{sinc}\frac{t}{0.25} + 3 \mathrm{sinc}\frac{t-0.25}{0.25} - 2 \mathrm{sinc}\frac{t-0.75}{0.25} \right)

オーディオ再生ではどうやって解決しているのか?

最初に「サンプリングしたデータをそのまま再生したときに元波形が再現できるとは一言も言ってない」と書きました。そして、そのまま再生した後でLPFを通し元の波形を再現することは困難であることを解説しました。では「そのまま再生しなければ」何か手立てがあるのでしょうか?

通常のオーディオ用DAC(DA変換)ではデータをそのまま再生するのではなく「オーバーサンプリング」という技術を使用してデータを再生します。

まず元々の信号を8倍ないしは16倍にアップサンプリングします。例えば、48kHzの信号の8倍オーバーサンプリングならば「384kHz」の信号に変換します。そしてsinc関数を使った畳み込み処理でデータを補完し再生します。*5

現実のオーディオ再生におけるオーバーサンプリングの実例は、「お気楽オーディオさんのNOSDACの記事」の下部にあるオシロ波形が参考になります。左側がオーバーサンプリングなし、右側が8倍オーバーサンプリングです。6つあるオシロ画像の右下の(疑似)インパルス応答を見ると、sinc関数でデータを再生している様子がよくわかります。

ただ、オーバーサンプリングでもやってることは一般的なデータ再生の項目で述べたことと原理的に近く、この方法も近似に過ぎません。

オーバーサンプリングについて

オーバーサンプリングの目的は元の波形を再現することだけではありません。主目的はノイズシェーピングと呼ばれる処理にあります。説明は割愛しますので気になる方は調べてください。

またオーバーサンプリング(やその前提であるアップサンプリング)にはデジタルフィルタが必須で、これもまた完璧な処理方法は存在しません。よってオーバーサンプリングの処理方法によって信号は変化してしまいますし、そこに問題を見出す人も居ます。

ただ現実問題として、今市販されてるほぼすべてのオーディオ用の再生装置はオーバーサンプリング技術を使用していて、オーバーサンプリングでないオーディオ再生装置を探すことは困難です。

しかし不思議なことに、オーディオ以外の用途ではオーバーサンプリングはあまり使われていません。

*5 : 畳み込みの点数(FIRのタップ数に相当)はDACチップによって異なりますし特に公表もされていませんが、遅延や回路規模の問題もありせいぜい100点ぐらいだと思われます。

まとめ

  • サンプリング定理より、サンプリング周波数の半分未満の信号しか含まない信号をサンプリングしたとき、理論上は元の信号を一意に再現できる。
  • しかし実際には「理想LPFが存在せず」「無限のデータ列を扱うことができない」ので、再現は困難である。
  • 現実のオーディオ再生装置では、オーバーサンプリングという技術で元波形に近づけようとしている。