静かなピアノ曲の減衰していく残響を16ビットで録音して再生すると、音が消える寸前に「ジリジリ」「ザラザラ」という、明らかに元の音とは違う質感のノイズが乗ることがあります。単に雑音が乗っているのなら「サーッ」という無害なヒスノイズになるはずなのに、実際には音程がついた金属的な響きに聞こえる。これは量子化誤差が「ランダムな雑音」ではなく、信号そのものから決まってしまう歪みだからです。
もっと極端な例を作れます。振幅が量子化ステップ（1 LSB）の半分しかない正弦波を量子化すると、出力は全サンプルがゼロになります。信号は完全に消滅します。ところが、量子化する直前にわざと雑音を足してから量子化すると、出力の平均をとったときに元の正弦波がちゃんと復活します。1 LSB の 1/20 の振幅の信号すら、量子化器を通り抜けて生き残ります。この一見して魔法のような操作がディザ (dither) です。
ディザとノイズシェーピングは、次のような場面で日常的に使われています。
- オーディオのビット深度変換: 24ビットで制作したマスターを16ビットのCDフォーマットへ落とすとき、単純に切り捨てると小信号部に歪みが出ます。TPDFディザ＋ノイズシェーピングを使うと、可聴帯域の実効ダイナミックレンジを16ビットの理論値より十数dB稼げます。
- ΔΣ型AD/DA変換器: 1ビットという極端に粗い量子化器で20ビット相当の分解能を出せるのは、オーバーサンプリングとノイズシェーピングの組み合わせがあるからです。さらに1ビット変調器特有のリミットサイクル（アイドルトーン）を潰すためにディザが併用されます。
- 画像処理・ディスプレイ: 8ビットパネルに10ビットの階調を表示するときのテンポラルディザリング、印刷の誤差拡散法（Floyd–Steinberg）も、本記事で扱う誤差フィードバック型ノイズシェーピングと数学的には同じ構造です。
- 計測・レーダー: 低振幅信号を粗いADCで捉えるとき、意図的にノイズを注入して平均化することで、ADCの分解能以下の信号を検出します。
本題に入る前に、これから作る装置の全体像を1枚の絵にしておきます。
信号は左から入り、まず過去の量子化誤差を差し引かれ（緑の経路）、次にディザ雑音を足され（青の経路）、最後に粗い量子化器を通ります。誤差は「量子化器の出力」から「ディザを足す前の節点 $u[n]$」を引いて作られるので、ディザ雑音まで含んだ量になります。下段の3段構えが本記事の物語です。①そのまま量子化すると誤差が信号に張り付いて高調波（歪み）になる、②ディザを足すと誤差が信号から独立な白色雑音に変わる、③ノイズシェーピングでその雑音を帯域外へ動かす。量子化の粗さそのものは最後まで消えませんが、誤差の「性質」と「置き場所」は設計者が選べる、というのがこの図の主張です。
本記事の内容
- 量子化誤差がなぜ「雑音」ではなく「歪み」になるのか — 誤差の鋸歯波表現とフーリエ級数
- 正弦波を粗く量子化したときの高調波振幅をベッセル関数で導出する
- 非減算ディザの一般理論 — 全誤差の条件付き特性関数を導出する
- RPDF／TPDF の判定条件：$\Phi_\nu^{(k)}(\ell\Omega) = 0$ というモーメント条件
- RPDFに残る「雑音変調」を閉じた形で求める
- 誤差フィードバック型ノイズシェーパの伝達関数 $E(z)(1-z^{-1})^n$ の導出
- 帯域内雑音電力とオーバーサンプリング比の関係、Gerzon–Cravenの対数面積保存則
- Pythonで無ディザ／RPDF／TPDF／TPDF＋1次シェーピングの4条件を比較する
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
量子化誤差は本当に「雑音」なのか
量子化雑音の教科書的な扱いでは、誤差 $e$ は $[-\Delta/2, \Delta/2]$ に一様分布し、信号とは独立で、白色である、と仮定します。この仮定のもとで誤差電力は $\sigma_e^2 = \Delta^2/12$ となり、フルスケール正弦波に対する SNR は $6.02B + 1.76$ dB という有名な式が出てきます。
しかしこの仮定は、よく考えるとかなり怪しいものです。量子化器は決定論的な関数です。入力 $x$ を決めれば出力 $Q(x)$ は一意に決まり、誤差 $e = Q(x) – x$ も一意に決まります。そこにランダムさはどこにもありません。「一様分布の独立な雑音」という描像は、入力が量子化ステップに比べて十分大きく、かつ十分に激しく動き回るときに、結果的にそう見えるという近似でしかないのです。
入力が小さいとき、この近似は音を立てて崩れます。極端な例として、振幅 $A = 0.5\Delta$ の正弦波を考えます。四捨五入型（ミッドトレッド）の量子化器 $Q(x) = \Delta\,\mathrm{round}(x/\Delta)$ は、$|x| < \Delta/2$ の入力をすべて $0$ に潰します。したがって出力は恒等的にゼロ、誤差は $e = -x$、つまり誤差が信号そのものになります。これを「独立な雑音」と呼ぶのは無理があります。
もう少し大きい振幅、たとえば $A = 2.5\Delta$ ではどうでしょうか。出力は $-2\Delta, -\Delta, 0, \Delta, 2\Delta$ の5段の階段状になり、誤差は入力の周期に同期した周期波形になります。周期波形をフーリエ級数に展開すれば、そこには信号の高調波が現れます。実際にこの後の実験で見るように、3次高調波は基本波に対して $-21.9$ dBc という、耳では明らかに知覚できるレベルで立ちます。しかも高調波は信号周波数に追随して動くので、聴感上は「歪み」として、白色雑音よりずっと目立ちます。
この「決定論的である」という事実を、実際の波形で見ておきましょう。
左は誤差 $q(w) = Q(w)-w$ を入力 $w$ の関数として描いたものです。値が $\pm\Delta/2$ の間を規則正しく往復する鋸歯波で、ランダムな要素はどこにもありません。右は振幅 2.5 LSB の正弦波を量子化した様子で、出力（橙）は5段の階段になり、誤差（赤）は入力の周期にぴったり同期した規則的な波形になっています。雑音であればこんな整った形にはならないという一点だけで、「独立な白色雑音」というモデルが破綻していることが分かります。
さらに厄介なことに、誤差には信号と同じ周波数の成分も含まれます。後の実測では、無ディザの場合に基本波の振幅が $-0.853$ dB だけ縮んでいました。つまり量子化器は小信号に対して利得誤差まで持つ非線形素子として振る舞います。
まとめると、量子化誤差が問題を起こす理由は一つに集約されます。誤差が信号と相関していることです。相関があるから高調波が立ち、相関があるから利得が変わり、相関があるから信号が消えたときに誤差も一緒に消えて「音が途切れる」ような不自然な挙動が生まれます。ではこの相関を、どうすれば断ち切れるのでしょうか。まずは誤差そのものを数式で正確に書き下すところから始めます。
量子化誤差の鋸歯波表現とフーリエ級数
量子化器を $Q(w) = \Delta\,\mathrm{round}(w/\Delta)$ と定義します。ここで $\Delta$ は量子化ステップ（1 LSB）です。量子化誤差を
$$ \begin{equation} q(w) = Q(w) – w \end{equation} $$
と定義します。この $q(w)$ の形を考えてみましょう。$w$ が $-\Delta/2 < w < \Delta/2$ にあるとき $Q(w) = 0$ なので $q(w) = -w$。$w$ が $\Delta/2 < w < 3\Delta/2$ にあるとき $Q(w) = \Delta$ なので $q(w) = \Delta - w$。つまり $q$ は $-\Delta/2$ から $+\Delta/2$ の間を往復する周期 $\Delta$ の鋸歯波です。この「量子化誤差＝入力の鋸歯波関数」という見方が、以降のすべての議論の出発点になります。
周期関数ならフーリエ級数に展開できます。$\Omega = 2\pi/\Delta$ と置きます（$\Omega$ は「入力振幅の空間」における角周波数で、単位は 1/振幅であることに注意してください。時間周波数ではありません）。$q$ は奇関数なので正弦項だけが残ります。
$$ q(w) = \sum_{\ell=1}^{\infty} b_\ell \sin(\ell \Omega w) $$
係数 $b_\ell$ を求めます。基本区間 $(-\Delta/2, \Delta/2)$ では $q(w) = -w$ なので、
$$ b_\ell = \frac{2}{\Delta}\int_{-\Delta/2}^{\Delta/2} (-w)\sin(\ell\Omega w)\, dw $$
部分積分すると $\int w \sin(aw)\,dw = \sin(aw)/a^2 – w\cos(aw)/a$ です。$a = \ell\Omega$ と置くと $a\Delta/2 = \pi\ell$ なので $\sin(a\Delta/2) = 0$、$\cos(a\Delta/2) = (-1)^\ell$ となります。被積分関数は偶関数なので区間を半分にして2倍すると、
$$ \int_{-\Delta/2}^{\Delta/2} w\sin(aw)\,dw = 2\left[\frac{\sin(aw)}{a^2} – \frac{w\cos(aw)}{a}\right]_0^{\Delta/2} = -\frac{\Delta(-1)^\ell}{a} = -\frac{(-1)^\ell \Delta^2}{2\pi\ell} $$
最後の等号では $a = 2\pi\ell/\Delta$ を代入しました。符号を戻して $2/\Delta$ を掛けると
$$ \begin{equation} q(w) = \sum_{\ell=1}^{\infty} \frac{(-1)^\ell \Delta}{\pi \ell}\sin(\ell\Omega w), \qquad \Omega = \frac{2\pi}{\Delta} \end{equation} $$
を得ます。この式は「量子化器という強烈な非線形性が、入力に対する無限個の正弦波の重ね合わせで書ける」ことを意味します。$\ell = 1$ の項が主要項で、$\ell$ が大きい項は $1/\ell$ で減衰しますが、鋸歯波の不連続点付近では収束が遅く、ギブス現象が現れます。この級数がどう鋸歯波を作り上げるかを見てみます。
$\ell \le 1$（青）はただの正弦波ですが、項を足していくと黒い鋸歯波にみるみる近づきます。$\ell \le 60$（赤）ではほぼ完全に重なっており、式(2)が正しいことが目で確認できます。一方で不連続点（$w/\Delta = \pm0.5, \pm1.5$）の直前では、項数をいくら増やしても約9%のオーバーシュートが消えません。これがギブス現象で、高次項がいつまでも効き続ける＝高調波が高い次数まで残ることの図的な表現になっています。
正弦波入力に対する高調波振幅
いよいよ本題です。入力が正弦波 $w(t) = A\sin\theta$（$\theta = \omega_0 t$）のとき、誤差は
$$ q(t) = \sum_{\ell=1}^{\infty} \frac{(-1)^\ell \Delta}{\pi\ell}\sin\!\left(\ell\Omega A \sin\theta\right) $$
となります。ここで「正弦波の正弦波」という入れ子構造が出てきました。これを展開するのがヤコビ–アンガー展開です。
$$ \sin(z\sin\theta) = 2\sum_{m:\text{奇数}} J_m(z)\sin(m\theta) $$
$J_m$ は第1種ベッセル関数です。$z = \ell\Omega A$ を代入して $\ell$ の和と $m$ の和を入れ替えると、
$$ \begin{equation} q(t) = \sum_{m:\text{奇数}} \left[\sum_{\ell=1}^{\infty} \frac{2(-1)^\ell \Delta}{\pi\ell} J_m\!\left(\frac{2\pi \ell A}{\Delta}\right)\right] \sin(m\theta) \end{equation} $$
角括弧の中が $m$ 次高調波の振幅です。ここから重要なことが3つ読み取れます。第一に、奇数次高調波しか出ないこと（正弦波と量子化器が両方とも原点対称だからです）。第二に、高調波振幅がベッセル関数を通じて $A/\Delta$ に依存し、しかも $J_m$ が振動する関数なので、振幅を少し変えるだけで高調波の大きさが激しく変動すること。第三に、$m = 1$ の項も存在すること、つまり基本波と同じ周波数の誤差成分があり、これが実効的な利得誤差を生むことです。
$A = 2.5\Delta$ の場合、この式を数値的に評価すると3次高調波の振幅は $0.1996\,\Delta$（$\ell$ を40万項まで足した値）で、後述する数値実験でFFTから測った $0.1999\,\Delta$ とよく一致します。$\ell$ に関する和は $J_m(z) \sim \sqrt{2/(\pi z)}$ の漸近形から項が $\ell^{-3/2}$ でしか減衰せず、収束が非常に遅いので、数値確認の際は項数に注意が必要です。
この式が予測する「高調波の暴れ方」を、入力振幅の関数として描いてみます。
横軸が入力正弦波の振幅 $A/\Delta$、縦軸が誤差に含まれる各次成分の振幅です。$A < 0.5$ LSB では出力が恒等的にゼロなので $m=1$ の成分が $A$ そのもの（誤差 $=-x$）として直線的に立ち上がり、他の次数はゼロです。$A$ が量子化レベルをまたぐたびに曲線が不連続に跳び、その後もベッセル関数の振動を反映して激しく上下します。振幅を1%変えただけで3次高調波が数dB動くという、線形システムではありえない振る舞いです。破線が本記事の実験条件 $A=2.5$ LSB で、3次成分（赤）が約 $0.20$ LSB を読み取れます。
ここまでで、量子化誤差が信号と完全に結びついていること、そして正弦波に対しては明確な高調波スペクトルを生むことを数式で確認しました。次は、この「決定論的な結びつき」を切断する道具、ディザを導入します。
ディザとは — わざと雑音を足すという逆転の発想
古い機械式の計算機や爆撃照準器では、部品の静止摩擦で針が動かなくなるのを防ぐために、意図的に微小な振動を与えていました。振動していれば摩擦は動摩擦になり、装置は滑らかに応答します。この振動が dither の語源です。量子化におけるディザもまったく同じ発想で、量子化器という「階段」に微小な揺さぶりをかけて、階段を滑らかな傾斜に変えてしまうのです。
具体的には、量子化する直前に乱数 $\nu$ を加えます。
$$ y = Q(x + \nu) $$
この一手で何が起きるのかを、量子化器の入出力特性で見るのが一番分かりやすいでしょう。
赤がディザなしの $Q(x)$ で、1 LSB刻みの階段です。この階段こそが非線形性であり、高調波の発生源でした。ところがディザを足して出力の期待値 $E[Q(x+\nu)]$ をとると、RPDF（橙の太破線）でもTPDF（青）でも、点線で描いた理想の直線 $y=x$ にぴったり重なります。確率的な揺さぶりを加えて平均をとるだけで、階段が傾斜に化けるわけです。ここで「平均」は時間平均でも実現できるので、後段でローパスを掛ければ実際にこの直線的な応答が得られます。
$\nu$ は $x$ とは独立な確率変数で、これをディザと呼びます。加えたディザを受信側で引き算できる場合を減算ディザ (subtractive dither)、引き算しない場合を非減算ディザ (non-subtractive dither, NSD) と呼びます。実際のAD変換やビット深度変換ではディザの実現値を後段に伝える手段がないので、ほぼ常に非減算ディザです。本記事もNSDを主役にします。
非減算ディザにおける全誤差を
$$ \begin{equation} \varepsilon = y – x = Q(x + \nu) – x = q(x + \nu) + \nu \end{equation} $$
と定義します。最後の変形では $Q(w) = w + q(w)$ に $w = x + \nu$ を代入しました。ここが重要な点で、非減算ディザの誤差は「量子化器そのものの誤差 $q(x+\nu)$」と「加えたディザ $\nu$」の両方を含みます。だからディザを足せば必ず雑音は増えます。ディザは無料ではなく、SNRを対価として歪みを買い取る取引なのです。
ではどんな確率分布のディザを選べばよいのでしょうか。「とりあえず一様乱数」で十分なのか、それとも分布に条件があるのか。ここを厳密に決着させるのが、次に導出する条件付き特性関数の理論です。
非減算ディザの一般理論 — 全誤差の条件付き特性関数
目標
信号 $x$ を固定したときの全誤差 $\varepsilon$ の分布、すなわち条件付き分布 $p(\varepsilon \mid x)$ を調べます。もしこの分布が $x$ に依存しないなら、誤差は信号と完全に独立で、文字通りの「雑音」になります。ここでは分布そのものではなく、扱いやすい特性関数 $\Phi_{\varepsilon|x}(u) = E[e^{ju\varepsilon} \mid x]$ を計算します。特性関数を使う理由は、モーメントが原点での微分で得られるので、「何次のモーメントまで信号非依存にできるか」という問いに直接答えられるからです。
導出
$\varepsilon = Q(x+\nu) – x$ を代入し、$x$ は定数なので指数の外に出します。
$$ \Phi_{\varepsilon|x}(u) = E_\nu\!\left[e^{ju(Q(x+\nu)-x)}\right] = e^{-jux}\, E_\nu\!\left[e^{juQ(x+\nu)}\right] $$
次に $Q(w) = w + q(w)$ を使って指数を2つに分けます。
$$ e^{juQ(w)} = e^{juw}\,e^{juq(w)} $$
ここで $q(w)$ は周期 $\Delta$ の周期関数でしたから、$e^{juq(w)}$ もまた $w$ について周期 $\Delta$ の周期関数です。周期関数はフーリエ級数に展開できる、というのがこの導出の核心です。
$$ e^{juq(w)} = \sum_{\ell=-\infty}^{\infty} a_\ell(u)\, e^{j\ell\Omega w} $$
係数 $a_\ell(u)$ を求めます。基本区間 $(-\Delta/2, \Delta/2)$ では $q(w) = -w$ なので、
$$ a_\ell(u) = \frac{1}{\Delta}\int_{-\Delta/2}^{\Delta/2} e^{juq(w)} e^{-j\ell\Omega w}\,dw = \frac{1}{\Delta}\int_{-\Delta/2}^{\Delta/2} e^{-j(u+\ell\Omega)w}\,dw $$
この積分は矩形関数のフーリエ変換そのものです。$\beta = u + \ell\Omega$ と置くと $\frac{1}{\Delta}\cdot\frac{2\sin(\beta\Delta/2)}{\beta} = \operatorname{sinc}(\beta\Delta/2)$ となります。ここで $\operatorname{sinc}(z) = \sin z / z$（正規化しない定義）です。したがって
$$ a_\ell(u) = \operatorname{sinc}\!\left(\frac{(u+\ell\Omega)\Delta}{2}\right) $$
これを元の式に戻します。$e^{juw}$ と $e^{j\ell\Omega w}$ をまとめると
$$ e^{juQ(w)} = \sum_{\ell} \operatorname{sinc}\!\left(\frac{(u+\ell\Omega)\Delta}{2}\right) e^{j(u+\ell\Omega)w} $$
$w = x+\nu$ を代入して $\nu$ について期待値をとります。$E_\nu[e^{j\beta(x+\nu)}] = e^{j\beta x}\Phi_\nu(\beta)$（$\Phi_\nu$ はディザの特性関数）を使うと、
$$ E_\nu\!\left[e^{juQ(x+\nu)}\right] = \sum_\ell \operatorname{sinc}\!\left(\frac{(u+\ell\Omega)\Delta}{2}\right)\Phi_\nu(u+\ell\Omega)\, e^{j(u+\ell\Omega)x} $$
最後に前に付いていた $e^{-jux}$ を掛けます。$e^{-jux}e^{j(u+\ell\Omega)x} = e^{j\ell\Omega x}$ と簡約されて、
$$ \begin{equation} \Phi_{\varepsilon|x}(u) = \sum_{\ell=-\infty}^{\infty} \underbrace{\operatorname{sinc}\!\left(\frac{(u+\ell\Omega)\Delta}{2}\right)\Phi_\nu(u+\ell\Omega)}_{g_\ell(u)}\, e^{j\ell\Omega x} \end{equation} $$
という美しい結果を得ます。これがWannamaker・Lipshitz・Vanderkooyらによる非減算ディザ理論の中心式です。
式の読み方
この式は「信号非依存な部分」と「信号依存な部分」にきれいに分離しています。
$\ell = 0$ の項は $\operatorname{sinc}(u\Delta/2)\Phi_\nu(u)$ で、$x$ をまったく含みません。しかもこの形は、$[-\Delta/2,\Delta/2]$ 上の一様分布（その特性関数が $\operatorname{sinc}(u\Delta/2)$）とディザ $\nu$ を独立に足した確率変数の特性関数に一致します。つまり $\ell=0$ 項こそが、教科書的な「量子化雑音は一様分布の独立雑音」というモデル（PQN: pseudo quantization noise モデル）そのものです。
$\ell \neq 0$ の項はすべて $e^{j\ell\Omega x}$ という因子を持ち、$x$ の関数として周期 $\Delta$ で振動します。これらが歪みの正体です。ディザなしのとき $\Phi_\nu \equiv 1$ なので、$\ell\neq0$ の項は消えません。だから歪みが出るのです。
したがってディザ設計の目標は明確です。$\ell \neq 0$ の項の寄与をできるだけ消すような $\Phi_\nu$ を選ぶこと。ここで $\Phi_\nu$ をどう選んでも $\ell\neq0$ 項を完全にゼロにはできない（それには $\Phi_\nu$ が $|u| \ge \Omega/2$ で恒等的にゼロという、有限分散の分布では両立しない条件が要る）ことが知られています。非減算ディザでは「分布そのものを信号非依存にする」ことは原理的に不可能なのです。しかしモーメントを信号非依存にすることならできます。次にその条件を導きます。
モーメント条件の導出
特性関数の $m$ 階微分から $m$ 次モーメントが得られます。
$$ E[\varepsilon^m \mid x] = (-j)^m \left.\frac{d^m \Phi_{\varepsilon|x}(u)}{du^m}\right|_{u=0} = (-j)^m \sum_\ell g_\ell^{(m)}(0)\, e^{j\ell\Omega x} $$
したがって「$m$ 次モーメントが $x$ に依存しない」ための条件は、すべての $\ell\neq0$ について $g_\ell^{(m)}(0) = 0$ となることです。$g_\ell(u) = S(u+\ell\Omega)\Phi_\nu(u+\ell\Omega)$（$S(u) \equiv \operatorname{sinc}(u\Delta/2)$）は積なので、ライプニッツ則で展開します。
$$ g_\ell^{(m)}(0) = \sum_{k=0}^{m}\binom{m}{k} S^{(m-k)}(\ell\Omega)\, \Phi_\nu^{(k)}(\ell\Omega) $$
ここで決定的な事実に気づきます。$S(\ell\Omega) = \operatorname{sinc}(\pi\ell) = 0$（$\ell \neq 0$）です。$\ell\Omega\Delta/2 = \pi\ell$ となり、$\sin(\pi\ell) = 0$ だからです。つまり和の中で $k = m$ の項（$S^{(0)}(\ell\Omega)\Phi_\nu^{(m)}(\ell\Omega)$）は自動的に消えます。残るのは $k = 0, 1, \dots, m-1$ の項だけです。よって次の十分条件が得られます。
$$ \begin{equation} \Phi_\nu^{(k)}(\ell\Omega) = 0 \quad \text{for all } \ell \neq 0,\ k = 0,1,\dots,m-1 \end{equation} $$
言葉にすると、ディザの特性関数が $\Omega$ の整数倍の点で $m$ 位の零点を持てば、全誤差の $m$ 次までのモーメントが信号から独立になる、ということです。$\Omega$ の整数倍という点が現れるのは、量子化器の階段が周期 $\Delta$ だからで、$\Omega = 2\pi/\Delta$ はその「サンプリング周波数」に相当します。
RPDFディザ — 1次モーメントだけが救われる
$\nu$ を $[-\Delta/2, \Delta/2]$ の一様分布（RPDF: Rectangular Probability Density Function）にとります。特性関数は
$$ \Phi_\nu(u) = \operatorname{sinc}\!\left(\frac{u\Delta}{2}\right) $$
$u = \ell\Omega$ では $\operatorname{sinc}(\pi\ell) = 0$（$\ell\neq0$）なので $m=1$ の条件を満たします。よって平均が信号非依存になり、$E[\varepsilon\mid x] = E[\nu] = 0$ が全ての $x$ で成り立ちます。歪み（DC的な偏り、そして1次の意味での信号相関）は消えます。
しかし $m=2$ の条件は満たしません。$\operatorname{sinc}$ の導関数は $\operatorname{sinc}'(z) = (z\cos z – \sin z)/z^2$ で、$z=\pi\ell$ では $\operatorname{sinc}'(\pi\ell) = (-1)^\ell/(\pi\ell) \neq 0$ です。零点が単純零点なので、1階微分は生き残ります。したがって分散は信号に依存します。
実際に計算してみましょう。$m=2$ のとき、$S(\ell\Omega)=\Phi_\nu(\ell\Omega)=0$ より $k=0$ と $k=2$ の項が消え、$k=1$ の項だけが残ります。
$$ g_\ell”(0) = 2 S'(\ell\Omega)\Phi_\nu'(\ell\Omega) $$
$S(u) = \operatorname{sinc}(u\Delta/2)$ より $S'(\ell\Omega) = (\Delta/2)\operatorname{sinc}'(\pi\ell) = \frac{\Delta(-1)^\ell}{2\pi\ell}$、RPDFでは $\Phi_\nu’ $ も同じ値です。掛けて2倍すると
$$ g_\ell”(0) = 2\left(\frac{\Delta(-1)^\ell}{2\pi\ell}\right)^2 = \frac{\Delta^2}{2\pi^2\ell^2} $$
$\ell=0$ の項は $g_0”(0) = -(\Delta^2/12 + \sigma_\nu^2)$ です（$\operatorname{sinc}(z) = 1 – z^2/6 + \cdots$ から $S”(0) = -\Delta^2/12$、$\Phi_\nu”(0) = -\sigma_\nu^2$ を使いました）。RPDFでは $\sigma_\nu^2 = \Delta^2/12$ です。$E[\varepsilon^2|x] = -\sum_\ell g_\ell”(0)e^{j\ell\Omega x}$ に代入し、$\ell$ と $-\ell$ を組にして余弦にまとめると
$$ E[\varepsilon^2 \mid x] = \frac{\Delta^2}{6} – \frac{\Delta^2}{\pi^2}\sum_{\ell=1}^{\infty}\frac{\cos(\ell\Omega x)}{\ell^2} $$
この級数は既知の公式 $\sum_{\ell\ge1}\cos(\ell\theta)/\ell^2 = \pi^2/6 – \pi\theta/2 + \theta^2/4$（$0\le\theta\le2\pi$）で閉じた形になります。$\theta = \Omega x = 2\pi s$（$s = x/\Delta$ の小数部）を代入して整理すると、$\Delta^2/6$ が打ち消し合って、
$$ \begin{equation} E[\varepsilon^2 \mid x] = \Delta^2 \, s(1-s), \qquad s = \left(\frac{x}{\Delta}\right) \bmod 1 \end{equation} $$
という驚くほど簡潔な結果を得ます。$s=0$（入力がちょうど量子化レベル上）で分散はゼロ、$s=1/2$（レベルの中間）で最大値 $\Delta^2/4$。つまり入力が1 LSB動くあいだに雑音の大きさがゼロから $\Delta^2/4$ まで振れるわけです。これが「雑音変調 (noise modulation)」で、音量に応じて雑音のレベルが息をするように変動するため、聴感上とても目立ちます。
$s=0$ で分散ゼロになることは直接確認できます。$x=0$ なら $x+\nu = \nu \in [-\Delta/2,\Delta/2]$ なので $Q$ は必ず $0$ を返し、誤差は常に $0$ です。逆に $s=1/2$、つまり $x=\Delta/2$ なら $x+\nu \in [0,\Delta]$ で出力は $0$ か $\Delta$ が半々、誤差は $\pm\Delta/2$ が半々なので分散は $\Delta^2/4$。式と完全に一致します。
TPDFディザ — 2次モーメントまで救う
$m=2$ の条件を満たすには、$\Phi_\nu$ が $\ell\Omega$ で2位の零点を持てばよい。最も簡単な作り方は、独立な2つのRPDFを足すことです。
$$ \nu = \nu_1 + \nu_2, \qquad \nu_1,\nu_2 \sim \mathcal{U}[-\Delta/2, \Delta/2] $$
独立な確率変数の和の特性関数は積なので、
$$ \Phi_\nu(u) = \operatorname{sinc}^2\!\left(\frac{u\Delta}{2}\right) $$
2乗になったことで $\ell\Omega$ での零点は2位になり、$\Phi_\nu(\ell\Omega) = 0$ かつ $\Phi_\nu'(\ell\Omega) = 0$ が同時に成り立ちます。一様分布同士の畳み込みは三角分布なので、$\nu$ の密度は $[-\Delta, +\Delta]$ 上の二等辺三角形になります。これがTPDF (Triangular PDF) ディザです。ピーク間振幅が $2\Delta$、つまり 2 LSB であることが重要で、これより小さくても大きくても2位零点の条件は崩れます。
2種類のディザの形と、そのときに全誤差がとる値の確率を並べてみます。
左はディザ自身の密度で、RPDF（橙）が幅1 LSBの一様分布、TPDF（青）が幅2 LSBの三角分布です。中央と右は「入力 $x$ を固定したときに全誤差 $\varepsilon$ がとる値とその確率」で、$\varepsilon$ は $k\Delta – x$ という飛び飛びの値しかとらないため、連続密度ではなく確率質量として描いています。中央のRPDFでは $x=0$ のとき誤差が確率1でゼロ（2乗平均 $0.000$）なのに、$x=0.5\Delta$ では $\pm0.5$ が半々になり2乗平均が $0.250$ へ跳ね上がります。右のTPDFでは分布の形自体は $x$ とともに変わるものの、凡例の2乗平均は3条件すべて $0.250$ で固定されています。分布は独立化できないが、モーメントは独立化できるという理論の主張が、そのまま数値として現れています。
TPDFのもとで平均と分散は次のようになります。$\ell\neq0$ の寄与が $m\le2$ で消えるので、$\ell=0$ 項だけを見ればよく、
$$ E[\varepsilon\mid x] = 0, \qquad E[\varepsilon^2 \mid x] = \frac{\Delta^2}{12} + \sigma_\nu^2 = \frac{\Delta^2}{12} + \frac{\Delta^2}{6} = \frac{\Delta^2}{4} $$
$\sigma_\nu^2 = 2\times\Delta^2/12 = \Delta^2/6$ を使いました。平均も分散も $x$ に一切依存しない定数です。歪みも雑音変調も、両方まとめて消えました。
代償は雑音電力です。$\Delta^2/4 = 3\times(\Delta^2/12)$ なので、無ディザ時の理論値に対して3倍、すなわち $10\log_{10}3 = 4.77$ dB のSNR劣化です。RPDFの場合は $\Delta^2/12 + \Delta^2/12 = \Delta^2/6$ で 3.01 dB の劣化にとどまりますが、その代わり雑音変調が残ります。4.77 dBを払って完全な信号独立性を買うか、3.01 dBで妥協して雑音変調を残すかという設計判断になります。オーディオでは前者が標準です。
この論法は一般化できます。$m$ 次モーメントまで信号非依存にするには $\ell\Omega$ に $m$ 位の零点が必要で、それには独立なRPDFを $m$ 個畳み込んだディザ（分散 $m\Delta^2/12$）を使えばよい。全誤差の分散は $(m+1)\Delta^2/12$、劣化は $10\log_{10}(m+1)$ dB です。$m=3$ にすると 6.02 dB の劣化になり、得られるもの（3次モーメントの独立性）に対して代償が大きすぎるため、実用上は $m=2$（TPDF）で止めるのが定石です。
補足：減算ディザなら代償はゼロ
もしディザの実現値を後段が知っていて引き算できるなら、話はまったく変わります。減算ディザの誤差は $\varepsilon_{\text{sub}} = Q(x+\nu) – \nu – x = q(x+\nu)$ だけで、$\nu$ が加わりません。RPDFディザを使えば $x+\nu$ を $\Delta$ で割った余りが完全に一様分布になるので、$q(x+\nu)$ は $x$ と独立に $[-\Delta/2,\Delta/2]$ の一様分布になります（Schuchmanの条件）。誤差電力はきっかり $\Delta^2/12$、劣化はゼロです。ただし送受でディザ系列を共有する必要があるため、AD変換やファイル配布では使えません。だからこそ非減算ディザの理論が必要になるわけです。
ここまでで、ディザによって歪みを「素性の良い白色雑音」に変換できることが分かりました。しかし雑音は増えています。では増えた雑音を、聞こえない場所へ押しやることはできないでしょうか。それがノイズシェーピングです。
ノイズシェーピングの理論
誤差フィードバックという発想
量子化誤差は本質的に減らせません。しかし「どの周波数に置くか」は選べます。もし量子化誤差を測ることができれば、次のサンプルを量子化するときにその誤差を先取りして差し引いておける。この直感を回路にしたのが誤差フィードバック型ノイズシェーパです。
ここで「誤差を測れる」というのが効く場面を確認しておきます。24ビットのデータを16ビットに落とす再量子化のような、入力も出力もディジタルな処理では、量子化器の入力と出力の差はその場で計算できます。誤差は完全に既知です。一方、アナログ入力のADCでは量子化器の入力（アナログ値）が分からないので、この構造は直接使えません。ADCではループフィルタを量子化器の前に置く形（ΔΣ変調器）で等価な効果を得ます。本記事の主対象は前者、ディジタル再量子化です。
構造を信号の流れとして描くと次のようになります。
注目すべきは緑の経路の起点です。誤差は「量子化器の出力 $y[n]$」から「量子化器の入力節点 $u[n]$」を引いて作られており、$u$ の分岐点がディザ加算の手前にあります。引かれるのがディザを足す前の値なので、$\epsilon[n]$ にはディザ雑音も丸ごと含まれ、後で導く雑音伝達関数によってディザ雑音まで整形されます。信号 $x[n]$ の経路には量子化器以外に何も掛からないので、信号伝達関数は恒等的に1のままです。信号の通り道と雑音の通り道が構造的に分離されていることが、この図から読み取れる最大のポイントです。
伝達関数の導出
信号 $x[n]$、量子化器の入力節点 $u[n]$、出力 $y[n]$、量子化器で生じた誤差 $\epsilon[n]$ とします。構造は次の2式で定義されます。
$$ u[n] = x[n] – \sum_{k=1}^{n_{\text{ord}}} c_k\, \epsilon[n-k], \qquad y[n] = Q\big(u[n] + d[n]\big) $$
$d[n]$ はディザです。そして誤差を量子化器の入力節点 $u[n]$ を基準に定義します。
$$ \epsilon[n] = y[n] – u[n] $$
この定義がポイントで、ディザ $d[n]$ の寄与も $\epsilon[n]$ に含まれます。こうしておくと、後で見るようにディザ雑音もシェーピングの恩恵を受けます。
$y[n] = u[n] + \epsilon[n]$ に $u[n]$ の定義を代入すると、
$$ y[n] = x[n] – \sum_{k=1}^{n_{\text{ord}}} c_k \epsilon[n-k] + \epsilon[n] $$
両辺をz変換します。$\epsilon[n-k] \to z^{-k}E(z)$（時間シフト性質）を使うと
$$ Y(z) = X(z) + E(z)\left(1 – \sum_{k=1}^{n_{\text{ord}}} c_k z^{-k}\right) $$
すなわち
$$ \begin{equation} Y(z) = \underbrace{1}_{\text{STF}}\cdot X(z) + \underbrace{H_e(z)}_{\text{NTF}}\cdot E(z), \qquad H_e(z) = 1 – \sum_{k=1}^{n_{\text{ord}}} c_k z^{-k} \end{equation} $$
ここに、この構造の素晴らしさが凝縮されています。信号伝達関数 (STF) はきっかり1、つまり信号は一切色付けされずにそのまま通ります。一方で雑音伝達関数 (NTF) $H_e(z)$ は係数 $c_k$ で自由に設計できます。信号には触れずに雑音のスペクトルだけを整形できる、というのが誤差フィードバックの本質です。
NTFを $(1-z^{-1})^n$ にする
低周波側の雑音を減らしたいので、$H_e(z)$ を $z=1$（直流）に零点を持つハイパスにします。最も素直な選択が
$$ H_e(z) = (1 – z^{-1})^{n} $$
です。$z=1$ に $n$ 重零点があるので直流付近の雑音は $\omega^n$ で押しつぶされます。係数 $c_k$ は二項展開から決まります。
$$ (1-z^{-1})^n = \sum_{k=0}^{n}\binom{n}{k}(-1)^k z^{-k} = 1 + \sum_{k=1}^{n}\binom{n}{k}(-1)^k z^{-k} $$
$H_e(z) = 1 – \sum_k c_k z^{-k}$ と係数を比較すると
$$ c_k = (-1)^{k+1}\binom{n}{k} $$
$n=1$ なら $c_1 = 1$（前サンプルの誤差をそのまま引く）。$n=2$ なら $c_1 = 2,\ c_2 = -1$。$n=3$ なら $c_1=3, c_2=-3, c_3=1$。実装は数行の遅延と加算だけで済みます。
単位円上での振幅特性を求めます。$z = e^{j\omega}$ を代入し、$1 – e^{-j\omega}$ を半角に整理します。
$$ 1 – e^{-j\omega} = e^{-j\omega/2}\left(e^{j\omega/2} – e^{-j\omega/2}\right) = e^{-j\omega/2}\cdot 2j\sin\frac{\omega}{2} $$
絶対値をとると $|1-e^{-j\omega}| = 2|\sin(\omega/2)|$ なので、
$$ \begin{equation} \left|H_e(e^{j\omega})\right| = \left(2\sin\frac{\omega}{2}\right)^{n} \end{equation} $$
$\omega \to 0$ で $\omega^n$ として消え、$\omega = \pi$（ナイキスト）で $2^n$ という最大値をとります。1次なら高域で6 dB持ち上がり、2次なら12 dB持ち上がります。低域で掘った分は必ず高域に積み上がるのです。
帯域内雑音電力とオーバーサンプリング
信号として使う帯域を $|\omega| < \omega_B$ とし、オーバーサンプリング比を $\mathrm{OSR} = \pi/\omega_B$ と定義します（時間周波数で言えば $\mathrm{OSR} = f_s/(2f_B)$）。分散 $\sigma_\epsilon^2$ の白色な誤差の両側パワースペクトル密度は $\sigma_\epsilon^2/(2\pi)$ なので、帯域内雑音電力は
$$ P_{\text{in}} = \frac{\sigma_\epsilon^2}{2\pi}\int_{-\omega_B}^{\omega_B}\left|H_e(e^{j\omega})\right|^2 d\omega = \frac{\sigma_\epsilon^2}{\pi}\int_0^{\omega_B}\left(2\sin\frac{\omega}{2}\right)^{2n} d\omega $$
$n=0$（シェーピングなし）なら $P_{\text{in}} = \sigma_\epsilon^2\,\omega_B/\pi = \sigma_\epsilon^2/\mathrm{OSR}$ で、これはオーバーサンプリング単独の効果（OSRを2倍にすると3 dB改善）です。
$n=1$ の場合は倍角公式 $4\sin^2(\omega/2) = 2(1-\cos\omega)$ を使うと厳密に積分できます。
$$ P_{\text{in}} = \frac{2\sigma_\epsilon^2}{\pi}\int_0^{\omega_B}(1-\cos\omega)\,d\omega = \frac{2\sigma_\epsilon^2}{\pi}\left(\omega_B – \sin\omega_B\right) $$
一般の $n$ では、$\omega_B$ が小さいとき $2\sin(\omega/2)\approx\omega$ と近似できて、
$$ \begin{equation} P_{\text{in}} \approx \frac{\sigma_\epsilon^2}{\pi}\cdot\frac{\omega_B^{2n+1}}{2n+1} = \sigma_\epsilon^2\,\frac{\pi^{2n}}{(2n+1)\,\mathrm{OSR}^{2n+1}} \end{equation} $$
シェーピングなしの $\sigma_\epsilon^2/\mathrm{OSR}$ との比をとると、改善率は $\pi^{2n}/\big((2n+1)\mathrm{OSR}^{2n}\big)$ です。$n=1,\ \mathrm{OSR}=8$ なら $\pi^2/(3\times64) = 0.0514$、つまり $-12.9$ dB。$n=2$ なら $\pi^4/(5\times4096) = 0.00476$ で $-23.2$ dB。厳密積分でも $-12.92$ dB、$-23.31$ dB とほぼ一致します。
$\mathrm{OSR}$ を2倍にしたときの改善は $2^{2n+1}$ 倍、すなわち $(2n+1)\times 3.01$ dB です。1次シェーピングならOSR倍増ごとに9 dB（分解能で1.5ビット）、2次なら15 dB（2.5ビット）稼げます。ΔΣ変調器の分解能がオーバーサンプリングで急速に伸びる理由がこの式です。
タダ飯はない — 対数面積保存則
ただし、帯域内で稼いだ分はどこかで払っています。全帯域の雑音電力の倍率は
$$ \frac{1}{\pi}\int_0^{\pi}\left(2\sin\frac{\omega}{2}\right)^{2n} d\omega = \binom{2n}{n} $$
で、$n=1$ なら2倍（+3.01 dB）、$n=2$ なら6倍（+7.78 dB）、$n=3$ なら20倍（+13.0 dB）です。次数を上げるほど帯域外の雑音は激しく増え、後段のアナログ回路の非線形性で折り返してくるなどの実害が出ます。この収支を1枚にまとめておきます。
左は帯域内雑音電力をOSRの関数として描いたもので、対数軸上で直線になり、その傾きが次数とともに急峻になります。$n=0$ はOSR倍増あたり3.0 dBですが、$n=1$ で9.0 dB、$n=2$ で15.0 dB、$n=3$ で21.1 dBです。右が「タダ飯はない」の可視化で、赤が全帯域の雑音電力の増加、青がOSR=8での帯域内雑音の変化です。$n=2$ では帯域内を $-23.3$ dB 掘る代わりに全体では $+7.8$ dB 積み増しており、掘った深さと積んだ高さは別々の量だが、必ずセットで動くことが分かります。次数を上げるほど右の赤棒が伸び続ける点が、実装上の歯止めになります。
より深い制約がGerzon–Cravenの定理です。$H_e(z)$ が因果的かつ最小位相で $h_e[0]=1$（モニック）であるとき、
$$ \frac{1}{2\pi}\int_{-\pi}^{\pi}\ln\left|H_e(e^{j\omega})\right| d\omega = \ln|h_e[0]| = 0 $$
が成り立ちます。対数振幅特性の面積は必ずゼロ、つまりdBスケールで見たとき、ある周波数で雑音を下げたら別の周波数で同じ面積だけ上げなければならない。これはノイズシェーピングの根本的な保存則で、$(1-z^{-1})^n$ もこの条件を満たします（数値積分すると確かに $0$ になります）。だから実用的なノイズシェーパの設計とは「どこを掘ってどこに積むか」を、聴覚の等ラウドネス曲線などの重み付けに基づいて最適化する問題になります。単純な $(1-z^{-1})^n$ は「直流付近を最大限掘る」という一つの極端な選択にすぎません。
ディザはループの内側に入れる
先ほど誤差を $\epsilon[n] = y[n] – u[n]$ と定義したことを思い出してください。$y[n] = Q(u[n]+d[n])$ なので
$$ \epsilon[n] = Q(u[n]+d[n]) – u[n] = q(u[n]+d[n]) + d[n] $$
これは非減算ディザの全誤差そのものです。つまり $\epsilon$ にはディザ雑音も含まれており、それがNTFで整形されて出力に現れます。TPDFディザで払った 4.77 dB の代償も、シェーピングによって帯域外へ追い出されるわけです。ディザをループの外（ノイズシェーパの後段）で足してしまうと、この恩恵は受けられません。
またディザ理論の適用条件も満たされています。$d[n]$ は独立同分布で、$u[n]$ は過去の誤差 $\epsilon[n-k]$（したがって過去の $d[n-k]$）にしか依存しないので、$d[n]$ と $u[n]$ は独立です。よって「$u[n]$ を条件付けたときのモーメントが $u$ に依存しない」という結論がそのまま成立します。
理論の準備はすべて整いました。ここからは実際に数値で確かめます。
Pythonでの実装
以下では、量子化ステップを $\Delta = 1$（1 LSBを単位にとる）として、サンプリング周波数 $f_s = 384$ kHz、信号帯域 $0$–$24$ kHz（$\mathrm{OSR}=8$）、信号は振幅 $A = 2.5$ LSB の約1 kHz正弦波という設定で統一します。FFTの窓漏れを避けるため、周期がFFT長にちょうど収まるコヒーレントサンプリング（$N = 65536$、$k_0 = 171$ 周期）を使います。
準備と誤差の鋸歯波
まず量子化器と誤差関数を定義し、誤差が入力の鋸歯波関数であること、正弦波入力に対して周期的な誤差波形になることを見ます。
import numpy as np
import matplotlib
import 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
D = 1.0 # 量子化ステップΔ（1 LSBを単位にとる）
quant = lambda w: D * np.round(w / D) # ミッドトレッド量子化器
# (1) 誤差 q(w)=Q(w)-w は周期Δの鋸歯波
w = np.linspace(-3*D, 3*D, 20001)
# (2) 正弦波入力に対する誤差波形（A=2.5 LSB）
A, N, k0 = 2.5*D, 1 << 16, 171
n = np.arange(N)
x = A * np.sin(2*np.pi*k0*n/N)
e_nodither = quant(x) - x
fig, ax = plt.subplots(1, 2, figsize=(13, 4))
ax[0].plot(w/D, (quant(w)-w)/D, lw=1.5)
ax[0].set_xlabel("入力 $w/\\Delta$ [LSB]"); ax[0].set_ylabel("誤差 $q(w)/\\Delta$")
ax[0].set_title("量子化誤差は入力の鋸歯波関数（周期 $\\Delta$）"); ax[0].grid(alpha=.3)
m = 600 # 1周期弱を表示
ax[1].plot(n[:m], x[:m], lw=1.2, label="入力 $x[n]$（振幅2.5 LSB）")
ax[1].plot(n[:m], quant(x)[:m], lw=1.2, label="量子化出力 $y[n]$（5段の階段）")
ax[1].plot(n[:m], e_nodither[:m], lw=1.0, label="誤差 $\\varepsilon[n]$（周期波形）")
ax[1].set_xlabel("サンプル番号 $n$"); ax[1].set_ylabel("振幅 [LSB]")
ax[1].set_title("小振幅正弦波の量子化：誤差は信号に同期した周期波形")
ax[1].legend(fontsize=9); ax[1].grid(alpha=.3)
plt.tight_layout(); plt.show()
print("誤差の全電力 / (Δ²/12) =", np.mean(e_nodither**2)/(D**2/12))
左のグラフは、量子化誤差が入力の関数として完全に決定論的な鋸歯波であることを示しています。ランダムさはどこにもありません。右のグラフでは、出力が5段の階段になり、誤差が入力の周期に同期した規則的な三角波状の波形になっているのが見えます。ランダムな雑音であればこんな規則正しい形にはなりません。最後の出力は約 $1.23$ で、誤差電力が理論値 $\Delta^2/12$ から2割以上ずれていることを示します。小信号では「一様分布モデル」が成立していない証拠です。
無ディザの高調波スペクトル
規則的な誤差波形は、周波数領域では高調波として現れます。式(3)のベッセル関数による予測と、FFTで測った値を突き合わせます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jv
fs = 384000.0
freq = np.fft.rfftfreq(N, 1/fs)
spec_e = np.abs(np.fft.rfft(e_nodither)) / N * 2 # 片側振幅スペクトル
# ベッセル関数による m 次高調波振幅の理論値（式(3)）
def harm_theory(m, A, L=400000):
l = np.arange(1, L+1)
return abs(np.sum((-1)**l * (2*D/(np.pi*l)) * jv(m, 2*np.pi*l*A/D)))
print("次数 FFT実測[LSB] 理論(ベッセル)[LSB]")
for m in [3, 5, 7, 9]:
print(f" {m} {spec_e[m*k0]:.5f} {harm_theory(m, A):.5f}")
plt.figure(figsize=(11, 4.5))
plt.semilogy(freq/1e3, np.maximum(spec_e, 1e-8), lw=.7, color="C3")
for m in [3, 5, 7, 9, 11, 13, 15]:
plt.axvline(m*k0*fs/N/1e3, color="k", ls=":", lw=.6)
plt.axvline(24, color="C0", ls="--", lw=1.5, label="信号帯域端 24 kHz")
plt.xlim(0, 60); plt.ylim(1e-5, 1)
plt.xlabel("周波数 [kHz]"); plt.ylabel("誤差の振幅 [LSB]")
plt.title("無ディザ量子化の誤差スペクトル：奇数次高調波が林立する（点線が高調波位置）")
plt.legend(); plt.grid(alpha=.3); plt.show()
黒丸を付けた位置が奇数次成分で、1 kHzの基本波から23 kHzの23次まで、帯域内（24 kHz以下）に12本（基本波1本＋高調波11本）も並んでいます。振幅は3次で $0.20$ LSB、15次でもまだ $0.17$ LSB あり、次数が上がっても単調に減らない点がベッセル関数由来の特徴です。偶数次の位置（2次・4次…）の振幅は $10^{-14}$ 台、つまり完全な数値誤差レベルで、理論どおり偶数次が出ていないことも確認できます。マークの間にびっしり並ぶ低いレベルの線は、ナイキスト（192 kHz）を超える非常に高次の奇数次高調波が折り返してきたものです。
出力される表では、3次高調波が実測 $0.19995$ に対して理論 $0.19962$、5次が $0.11782$ に対して $0.11750$ と、ベッセル関数による予測が3〜4桁で一致します。式(3)の導出が正しいことの確認です（$\ell$ の和は $\ell^{-3/2}$ でしか減衰しないため、項数を減らすと数%ずれます）。グラフでは、点線で示した奇数次高調波の位置に鋭いスパイクが立ち、しかもその多くが信号帯域（24 kHz以下）の内側にあることが分かります。誤差電力を平坦な雑音とみなす議論が成り立たないのは明白です。
条件付きモーメント — 理論の直接検証
理論の中心である「$E[\varepsilon\mid x]$ と $E[\varepsilon^2\mid x]$ が信号に依存するか」を、ディザ分布に対する数値積分で直接確かめます。
import numpy as np
import matplotlib.pyplot as plt
def cond_moments(xs, kind, ngrid=200001):
"""入力xを固定したときの全誤差の平均と2乗平均を数値積分で求める"""
if kind == "none":
e = quant(xs) - xs
return e, e**2
if kind == "rpdf": # 一様分布 [-Δ/2, Δ/2]
v = np.linspace(-D/2, D/2, ngrid); p = np.ones(ngrid)
else: # 三角分布 [-Δ, Δ]
v = np.linspace(-D, D, ngrid); p = np.maximum(0, 1 - np.abs(v)/D)
p = p / p.sum()
m1, m2 = [], []
for xx in np.atleast_1d(xs):
e = quant(xx + v) - xx
m1.append((e*p).sum()); m2.append((e**2*p).sum())
return np.array(m1), np.array(m2)
xs = np.linspace(-1.5*D, 1.5*D, 241)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.2))
for kind, lab, c in [("none", "ディザなし", "C3"), ("rpdf", "RPDF", "C1"), ("tpdf", "TPDF", "C0")]:
m1, m2 = cond_moments(xs, kind)
ax[0].plot(xs/D, m1/D, color=c, lw=1.8, label=lab)
ax[1].plot(xs/D, m2/D**2, color=c, lw=1.8, label=lab)
s = (xs/D) % 1.0
ax[1].plot(xs/D, s*(1-s), "k--", lw=1.0, label="RPDF理論 $s(1-s)$")
ax[1].axhline(1/4, color="C0", ls=":", lw=1.2, label="TPDF理論 $\\Delta^2/4$")
ax[0].set_title("条件付き平均 $E[\\varepsilon|x]$"); ax[1].set_title("条件付き2乗平均 $E[\\varepsilon^2|x]$")
for a, yl in zip(ax, ["$E[\\varepsilon|x]/\\Delta$", "$E[\\varepsilon^2|x]/\\Delta^2$"]):
a.set_xlabel("入力 $x/\\Delta$ [LSB]"); a.set_ylabel(yl); a.legend(fontsize=9); a.grid(alpha=.3)
plt.tight_layout(); plt.show()
左のグラフでは、ディザなしの平均が振幅 $\pm0.5$ LSB の鋸歯波として大きく振れているのに対し、RPDFとTPDFはどちらも完全にゼロの直線に張り付いています。これが $\Phi_\nu(\ell\Omega)=0$ の効果です。右のグラフが本題で、RPDFの2乗平均は入力が1 LSB進むあいだにゼロから $0.25\Delta^2$ まで大きく脈動し、黒破線で描いた理論式 $s(1-s)$ とぴったり重なります。一方TPDFは入力によらず $0.25\Delta^2$ の水平線で、分散まで信号から独立になっていることが確認できます。RPDFが「歪みは消せるが雑音変調は残る」中途半端な選択であることが、この1枚で分かります。
ディザの特性関数を見る
なぜTPDFで2次まで救われるのかを、特性関数の零点の位数として可視化します。
import numpy as np
import matplotlib.pyplot as plt
Omega = 2*np.pi/D
u = np.linspace(-3.2*Omega, 3.2*Omega, 4001)
sinc = lambda z: np.sinc(z/np.pi) # np.sinc は正規化版なのでπで割る
Phi_rpdf = sinc(u*D/2)
Phi_tpdf = sinc(u*D/2)**2
plt.figure(figsize=(11, 4.2))
plt.plot(u/Omega, Phi_rpdf, lw=1.8, label="RPDF: $\\Phi_\\nu(u)=\\mathrm{sinc}(u\\Delta/2)$（単純零点）")
plt.plot(u/Omega, Phi_tpdf, lw=1.8, label="TPDF: $\\Phi_\\nu(u)=\\mathrm{sinc}^2(u\\Delta/2)$（2位零点）")
plt.axhline(0, color="k", lw=.6)
for l in range(-3, 4):
if l != 0:
plt.axvline(l, color="C3", ls=":", lw=1.0)
plt.plot([], [], "C3:", label="判定点 $u=\\ell\\Omega$（$\\ell\\neq0$）")
plt.xlabel("$u/\\Omega$"); plt.ylabel("$\\Phi_\\nu(u)$")
plt.title("ディザの特性関数：$\\ell\\Omega$ での零点の位数が救えるモーメント次数を決める")
plt.legend(fontsize=9); plt.grid(alpha=.3); plt.show()
for name, Phi in [("RPDF", Phi_rpdf), ("TPDF", Phi_tpdf)]:
i = np.argmin(np.abs(u - Omega))
slope = (Phi[i+1]-Phi[i-1])/(u[i+1]-u[i-1])
print(f"{name}: Φ(Ω)={Phi[i]:+.2e}, Φ'(Ω)={slope:+.4f}")
赤い点線が判定点 $u=\ell\Omega$ です。RPDF（橙）はこの点で値がゼロになりますが、右の拡大図を見ると曲線は斜めに横切っており傾きが残っています。出力される $\Phi'(\Omega) = -0.1592$ がその値で、理論値 $\frac{\Delta(-1)^\ell}{2\pi\ell}\big|_{\ell=1} = -\frac{1}{2\pi} = -0.1592$ と一致します（黒破線がこの接線です）。この残った傾きが雑音変調の源です。TPDF（青）は同じ点で曲線が下から滑らかに接して値も傾きもゼロになり、$\Phi'(\Omega)$ の数値がほぼ $0$ と表示されます。零点が単純か2位かという一点だけで、雑音変調が残るか消えるかが決まるという理論の主張が、この図の形として直接見えています。
NTFの周波数特性
次にノイズシェーピング側です。次数ごとのNTF振幅と、帯域内雑音電力の理論値を確認します。
import numpy as np
import matplotlib.pyplot as plt
wgrid = np.linspace(1e-6, np.pi, 200000)
OSR = 8; wB = np.pi/OSR
plt.figure(figsize=(11, 4.4))
for order, c in zip([0, 1, 2, 3], ["C7", "C0", "C1", "C2"]):
H = np.abs(1 - np.exp(-1j*wgrid))**order
plt.plot(wgrid/np.pi, 20*np.log10(np.maximum(H, 1e-6)), color=c, lw=1.8,
label=f"$n={order}$" + ("（シェーピングなし）" if order == 0 else ""))
inband = np.trapz(H[wgrid <= wB]**2, wgrid[wgrid <= wB])/np.pi
total = np.trapz(H**2, wgrid)/np.pi
print(f"n={order}: 帯域内雑音倍率={inband:.6f}（対n=0で{10*np.log10(inband/(wB/np.pi)):+6.2f} dB）"
f", 全帯域電力倍率={total:.3f}, 対数面積={np.trapz(np.log(np.maximum(H,1e-300)), wgrid)/np.pi:+.4f}")
plt.axvspan(0, wB/np.pi, color="C0", alpha=.12, label="信号帯域（OSR=8）")
plt.xlabel("正規化角周波数 $\\omega/\\pi$"); plt.ylabel("$|H_e(e^{j\\omega})|$ [dB]")
plt.ylim(-80, 25); plt.title("雑音伝達関数 $H_e(z)=(1-z^{-1})^n$：低域を掘った分だけ高域に積み上がる")
plt.legend(fontsize=9); plt.grid(alpha=.3); plt.show()
出力の数値は、$n=1$ で帯域内雑音が $-12.92$ dB、$n=2$ で $-23.31$ dB となり、近似式 $\pi^{2n}/((2n+1)\mathrm{OSR}^{2n})$ による予測（$-12.89$ dB、$-23.22$ dB）とよく合います。同時に全帯域電力倍率が $2.000, 6.000, 20.00$ と二項係数 $\binom{2n}{n}$ に一致し、次数を上げるほど総雑音が急増することも確認できます。対数面積の欄はすべて $0.0000$ で、Gerzon–Cravenの保存則が数値的に成り立っています。グラフでは、影を付けた信号帯域の中で曲線が深く沈む代わりに、高域で $6n$ dB 持ち上がっている様子がはっきり見えます。
ノイズシェーパの実装と4条件の比較
いよいよ本命です。無ディザ／RPDF／TPDF／TPDF＋1次シェーピングの4条件でスペクトルと帯域内SNRを比較します。
import numpy as np
from math import comb
def requantize(x, dither="tpdf", order=0, seed=1):
"""誤差フィードバック型ノイズシェーパ付き再量子化器（ディザはループ内）"""
rng = np.random.default_rng(seed)
N = len(x)
if dither == "none":
d = np.zeros(N)
elif dither == "rpdf":
d = rng.uniform(-D/2, D/2, N)
else: # TPDF = 独立な2つのRPDFの和
d = rng.uniform(-D/2, D/2, N) + rng.uniform(-D/2, D/2, N)
c = np.array([(-1)**(k+1)*comb(order, k) for k in range(1, order+1)], float)
y = np.zeros(N); eps = np.zeros(N)
for i in range(N):
fb = sum(c[j-1]*eps[i-j] for j in range(1, order+1) if i-j >= 0)
u = x[i] - fb # 過去の誤差を先取りして差し引く
y[i] = quant(u + d[i]) # ディザを足してから量子化
eps[i] = y[i] - u # 誤差はディザ込みで定義（ループ内ディザ）
return y
実装は驚くほど短く、遅延つきの加算とループが1本あるだけです。eps[i] = y[i] - u としてディザを誤差の定義に含めている点が、前節で述べた「ディザをループの内側に置く」に対応します。ここに任意の係数 $c_k$ を入れれば任意のFIR型NTFが実現できます。
import numpy as np
band = N // 16 # 0〜fs/16 = 24 kHz を信号帯域とする
conds = [("無ディザ", "none", 0), ("RPDF", "rpdf", 0), ("TPDF", "tpdf", 0),
("TPDF+1次シェーピング", "tpdf", 1), ("TPDF+2次シェーピング", "tpdf", 2)]
results = {}
print(f"{'条件':<22}{'帯域内SNR':>10}{'誤差全電力':>14}{'最大スプリアス':>14}{'利得誤差':>10}")
for name, dith, order in conds:
y = requantize(x, dith, order)
err = y - x
P = np.abs(np.fft.rfft(err)/N)**2 * 2 # 片側パワースペクトル
Pn = P[1:band+1].sum() - P[k0] # 基本波ビンを除いた帯域内雑音
snr = 10*np.log10((A**2/2) / Pn)
sub = P[1:band+1].copy(); sub[k0-1] = 0
spur = 10*np.log10(sub.max()/(A**2/2))
gain = 20*np.log10(abs(np.fft.rfft(y)[k0]) / abs(np.fft.rfft(x)[k0]))
results[name] = (y, snr)
print(f"{name:<22}{snr:8.2f} dB{np.mean(err**2)/(D**2/12):11.2f}×Δ²/12"
f"{spur:11.1f} dBc{gain:9.3f} dB")
実行結果は次のようになります。
| 条件 | 帯域内SNR | 誤差全電力 | 帯域内最大スプリアス | 利得誤差 |
|---|---|---|---|---|
| 無ディザ | 17.33 dB | 1.23 × $\Delta^2/12$ | $-21.9$ dBc | $-0.853$ dB |
| RPDF | 21.30 dB | 2.16 × $\Delta^2/12$ | $-48.5$ dBc | $-0.002$ dB |
| TPDF | 19.91 dB | 3.01 × $\Delta^2/12$ | $-47.4$ dBc | $-0.002$ dB |
| TPDF＋1次シェーピング | 32.78 dB | 6.00 × $\Delta^2/12$ | $-56.2$ dBc | $-0.000$ dB |
| TPDF＋2次シェーピング | 43.35 dB | 17.96 × $\Delta^2/12$ | $-65.3$ dBc | $+0.000$ dB |
この表には理論の主張がすべて詰まっています。無ディザでは $-21.9$ dBc というとてつもなく大きなスプリアス（3次高調波）が立ち、しかも基本波の利得が $0.853$ dB も縮んでいます。ディザを足した瞬間にスプリアスは $-47$ dB 以下へ落ち、利得誤差も $0.002$ dB に消えます。誤差全電力は理論どおり RPDFで約2倍（実測2.16はRPDFの分散が信号依存であるため理論の2.00から少しずれます）、TPDFできっかり3.01倍。TPDFはRPDFよりSNRが1.4 dB悪いのですが、その対価として雑音変調が消えています。1次シェーピングを加えると帯域内SNRがTPDF単体から $12.86$ dB 改善し、理論値 $12.92$ dB と一致します。2次では $23.43$ dB 改善（理論 $23.31$ dB）ですが、誤差全電力は18倍にまで膨らんでおり、帯域外にどれだけ雑音を積み上げているかが分かります。
import numpy as np
import matplotlib.pyplot as plt
freq = np.fft.rfftfreq(N, 1/fs)
plt.figure(figsize=(12, 5.5))
for (name, _, _), c in zip(conds[:4], ["C3", "C1", "C0", "C2"]):
y = results[name][0]
P = np.abs(np.fft.rfft(y - x)/N)**2 * 2
# 対数ビンで平滑化して見やすくする（スプリアスは最大値で保持）
sm = np.convolve(P, np.ones(64)/64, mode="same")
plt.semilogx(freq[1:]/1e3, 10*np.log10(np.maximum(sm[1:], 1e-16)),
color=c, lw=1.3, label=f"{name}（帯域内SNR {results[name][1]:.1f} dB）")
plt.axvline(24, color="k", ls="--", lw=1.5)
plt.text(25, -55, "信号帯域端 24 kHz", fontsize=9)
plt.xlim(0.3, 192); plt.ylim(-130, -40)
plt.xlabel("周波数 [kHz]（対数軸）"); plt.ylabel("誤差パワースペクトル密度 [dB]")
plt.title("4条件の誤差スペクトル比較：ディザで歪みが白色雑音になり、シェーピングで帯域内が沈む")
plt.legend(fontsize=9, loc="lower right"); plt.grid(alpha=.3, which="both"); plt.show()
このグラフが記事全体の要約です。赤（無ディザ）は低域に高調波のスパイクが飛び出し、しかも床のレベルが平坦ではありません。橙（RPDF）と青（TPDF）はスパイクが消えて完全に平坦な白色スペクトルになり、青のほうが約1.4 dB高い床になっています（帯域内雑音電力の実測比 $3.19\times10^{-2}$ 対 $2.32\times10^{-2}$）。理論上の差は 4.77 dB 対 3.01 dB の 1.76 dB ですが、RPDFの誤差電力が信号依存で理論値より膨らむぶん、実測の差はそれより小さく出ます。緑（TPDF＋1次シェーピング）は、24 kHz以下で床が大きく沈み、そのぶん高域が右肩上がりに持ち上がっています。歪みを雑音に変え、その雑音を聞こえない場所へ動かすという2段構えの戦略が、1枚の図として現れています。
表の3つの指標を棒グラフに並べ替えると、支払いと見返りの関係がさらにはっきりします。
左のSNRだけを見るとRPDF（21.3 dB）がTPDF（19.9 dB）に勝っており、ディザは弱いほうが得に見えます。しかし中央のスプリアスでは無ディザの $-21.9$ dBc がディザ導入で $-47$ dBc 台へ一気に落ち、右の総電力ではTPDFがきっかり $3.01$ 倍（$=\Delta^2/4$）という理論値どおりの代償を払っていることが分かります。RPDFとTPDFの1.4 dB差の中身は「雑音変調が残るか消えるか」であり、SNRという単一指標には現れません。シェーピングを加えた右2本は、総電力を6倍・18倍に増やしながら帯域内SNRを32.8 dB・43.3 dBまで押し上げており、雑音を「減らす」のではなく「移動させる」技術であることが数字に出ています。
LSB以下の信号は本当に生き残るのか
最後に、冒頭で予告した「1 LSBより小さい信号がディザで復活する」ことを確かめます。
import numpy as np
rng = np.random.default_rng(3)
N2, k2 = 1 << 18, 684
n2 = np.arange(N2)
print("振幅[LSB] 無ディザ出力の基本波 TPDF出力の基本波")
for amp in [0.5, 0.2, 0.05]:
x2 = amp * np.sin(2*np.pi*k2*n2/N2)
y0 = quant(x2)
d = rng.uniform(-D/2, D/2, N2) + rng.uniform(-D/2, D/2, N2)
y1 = quant(x2 + d)
f = lambda s: abs(np.fft.rfft(s)[k2])/N2*2
print(f" {amp:<9.2f}{f(y0):>14.5f}{f(y1):>20.5f}")
左は振幅 0.05 LSB の入力（黒）と無ディザ出力（赤）で、出力は全サンプルがゼロの直線、つまり信号が完全に消滅しています。中央はTPDFディザを足した場合で、生の出力（灰色）は $0, \pm1$ を不規則に行き来する2値的な系列にしか見えませんが、信号周波数まわりだけを通す狭帯域フィルタを掛けると（青、見やすさのため8倍表示）、元の正弦波（黒破線）と位相まで一致した波形が現れます。右のFFT実測では、無ディザがすべて $0.00000$ なのに対しTPDFは $0.49836 / 0.19992 / 0.05043$ と目標値をほぼ正確に再現しています。情報は個々のサンプル値ではなく、跳び方の確率のなかに埋め込まれていることが3枚並べると腑に落ちます。
出力は、無ディザではどの振幅でも基本波が完全に $0.00000$（出力は全サンプルがゼロ）であるのに対し、TPDFディザでは $0.49836$、$0.19992$、$0.05043$ と、元の振幅がほぼ正確に再現されます。1 LSBの1/20という、量子化器が原理的に表現できないはずの振幅の信号が、統計的に符号化されて出力に残っているわけです。量子化器の出力は依然として整数値の飛び飛びの数列ですが、ディザが信号の情報を「どのタイミングでどちらのレベルへ跳ぶか」という確率へ埋め込んでいます。これがディザの最も直感に反する、そして最も重要な効能です。
ここまでで理論と数値実験は一通り揃いました。最後に、この理論を実際の再量子化器やコンバータに落とし込むときに引っかかりやすい点を確認しておきます。
実務上の注意点
理論を実装に落とすときに引っかかりやすい点を整理します。
ディザは無音でも切らない。 曲のフェードアウトや無音区間でディザを止めると、その瞬間から歪みが復活します。むしろ小信号ほどディザが必要なので、レベルに応じてディザをオン・オフする設計は逆効果になりがちです。ヒスノイズが気になる場合は、ディザを切るのではなくノイズシェーピングの重み付けで聴感上の目立ちを下げるのが正しい対処です。
TPDFの振幅は 2 LSB ピーク間でなければならない。 実装ミスで $\pm\Delta/2$ の三角分布（ピーク間 1 LSB）にしてしまうと、$\Phi_\nu(u) = \operatorname{sinc}^2(u\Delta/4)$ となって零点が $2\ell\Omega$ に移り、$\ell\Omega$ ではゼロになりません。歪みも雑音変調も残ります。独立な一様乱数を2回引いて足す、という手順を守るのが確実です。
ノイズシェーピングは再量子化器の話であってADCの話ではない。 誤差フィードバック構造は量子化誤差が観測できることを前提にします。アナログ入力のADCでは同じことができないので、ΔΣ変調器では積分器を量子化器の前に置き、フィードバックループ全体として同じNTFを実現します。詳しくはΔΣ変調とオーバーサンプリングを参照してください。
シェーピング次数はやみくもに上げない。 誤差フィードバック型はNTFがFIRなので線形モデルでは常に安定ですが、次数を上げると帯域外の雑音振幅が大きくなり、量子化器がオーバーロード（$\epsilon$ が $\pm\Delta/2$ を超える領域に入る）してループが暴れます。実用上は $(1-z^{-1})^n$ をそのまま使うのではなく、零点を可聴帯域内に分散させたり、NTFの帯域外ゲインを制限した設計にします。
$(1-z^{-1})^n$ は最適ではない。 対数面積保存則のもとで「どこを掘るか」は自由なので、実際の高音質ビット深度変換では、聴覚の最小可聴閾に沿って重み付けした最適NTFを設計します。可聴帯域の中でも3〜4 kHz付近を深く掘り、15 kHz以上に積み上げる形にすると、同じ総雑音電力でも聴感上のダイナミックレンジを大きく稼げます。この「重み付きの最適NTF設計」まで踏み込むと、話は最適フィルタ設計の領域に入っていきます。ここまでの流れを最後に整理しておきましょう。
まとめ
本記事では、ディザとノイズシェーピングの理論を導出から実装まで通して解説しました。
- 量子化誤差は雑音ではない: 誤差 $q(w) = Q(w)-w$ は入力の周期 $\Delta$ の鋸歯波であり、フーリエ級数 $\sum_\ell \frac{(-1)^\ell\Delta}{\pi\ell}\sin(\ell\Omega w)$ で書ける。正弦波入力ではヤコビ–アンガー展開により奇数次高調波が現れ、その振幅はベッセル関数 $J_m(2\pi\ell A/\Delta)$ の和で決まる
- 全誤差の条件付き特性関数: $\Phi_{\varepsilon|x}(u) = \sum_\ell \operatorname{sinc}\big(\frac{(u+\ell\Omega)\Delta}{2}\big)\Phi_\nu(u+\ell\Omega)e^{j\ell\Omega x}$。$\ell=0$ 項がPQNモデル、$\ell\neq0$ 項が歪みの正体
- モーメント条件: $\Phi_\nu^{(k)}(\ell\Omega)=0$（$k=0,\dots,m-1$、$\ell\neq0$）なら $m$ 次までのモーメントが信号非依存。$\operatorname{sinc}$ の零点のおかげで $k=m$ の項は自動的に消える
- RPDF は平均だけ、TPDF は分散まで: RPDF（$\Phi_\nu=\operatorname{sinc}$、単純零点）は歪みを消すが分散が $\Delta^2 s(1-s)$ と脈動する。TPDF（$\operatorname{sinc}^2$、2位零点）は分散も $\Delta^2/4$ の定数にする。代償は 3.01 dB と 4.77 dB
- ノイズシェーピング: 誤差フィードバックにより $Y = X + H_e(z)E$、STFは1のまま。$H_e(z)=(1-z^{-1})^n$ で $|H_e| = (2\sin(\omega/2))^n$、帯域内雑音は $\sigma^2\pi^{2n}/\big((2n+1)\mathrm{OSR}^{2n+1}\big)$、OSR倍増あたり $(2n+1)\times3$ dB改善
- タダ飯はない: 全帯域雑音は $\binom{2n}{n}$ 倍に増え、Gerzon–Cravenの定理により対数振幅の面積は常にゼロ。設計とは「どこを掘ってどこに積むか」の選択
- 実測: TPDF＋1次シェーピング（OSR=8）は無ディザ比で帯域内SNRが 15.5 dB改善し、スプリアスは $-21.9$ dBc から $-56.2$ dBc へ、利得誤差は $-0.853$ dB から $0.000$ dB へ改善した
次のステップとして、以下の記事も参考にしてください。