ΔΣノイズシェーピングの式解析

1. 基本式

ΔΣ変調器の出力は次式で表されます：

Y(z) = STF(z) * X(z) + NTF(z) * E(z)

• STF(z): 信号伝達関数 (Signal Transfer Function)
• NTF(z): 雑音伝達関数 (Noise Transfer Function)
• E(z): 量子化雑音（白色雑音仮定）

---

2. 一般化されたNTF

プログラム中の apply_ntf() で実装されているのは、

NTF(z) = (1 - z⁻¹)ⁿ

つまり n階差分フィルタ を量子化雑音にかけていることと等価です。
次数 n が増えるほど、低周波数域の雑音が強力に抑圧されます。

---

3. 周波数領域での近似

z = e^{jω} を代入すると、

NTF(e^{jω}) = (1 - e^{-jω})ⁿ

小さい周波数 ω に対して近似すると：

|NTF(e^{jω})| ≈ |jω|ⁿ

→ 雑音スペクトルは低周波域で ωⁿ に比例して減少します。

---

4. dBスケールでの傾き

パワースペクトル密度 (PSD) の傾きは：

20 * n [dB/decade] （振幅）

雑音電力を扱うため +10 dB が加わり、

slope_theory = 20 * n + 10 [dB/decade]

• 1次 (n=1): 30 dB/decade
• 2次 (n=2): 50 dB/decade
• 3次 (n=3): 70 dB/decade
• 4次 (n=4): 90 dB/decade

---

5. プログラム中の対応

• np.diff(y) → 差分処理（1 - z⁻¹）
• forループで order 回繰り返す → (1 - z⁻¹)ⁿ を適用
• welch() → PSD推定
• np.polyfit(log10(f), Pxx_dB, 1) → 傾き [dB/decade] を実測評価
• 理論値 theory_slope = 20 * order + 10 と比較

---

6. 式のまとめ

最終的な出力のノイズ成分は：

Y_noise(z) = (1 - z⁻¹)ⁿ * E(z)

周波数領域での雑音伝達特性は：

|Y_noise(e^{jω})|² ≈ ω^(2n)

→ 周波数が1桁上がるごとに雑音が 20n + 10 dB 増加（シェーピングの傾き）。

---

!
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch, lfilter

# ==============================
# Parameters
# ==============================
fs = 1.0          # normalized sampling freq
N = 2**16         # number of samples
N_order = 9       # <<< set here: order of ΔΣ modulator (例: 9)

np.random.seed(0)  # reproducibility
q = np.random.randn(N)  # white quantization noise

# ==============================
# Function: apply NTF = (1 - z^-1)^N_order
# ==============================
def apply_ntf_direct(noise, order):
    """Apply NTF(z) = (1 - z^-1)^order using FIR filter coefficients"""
    b = np.poly1d([1, -1])**order   # FIR polynomial expansion
    b = b.coeffs
    y = lfilter(b, [1], noise)
    return y

# ==============================
# Apply NTF
# ==============================
y = apply_ntf_direct(q, N_order)

# PSD estimation
f, Pxx = welch(y, fs=fs, nperseg=4096)
Pxx_dB = 10 * np.log10(Pxx + 1e-20)

# ==============================
# Theory line (slope = 20*N_order+10)
# ==============================
theory_slope = 20 * N_order + 10
f0 = 1e-3
y0 = np.interp(f0, f, Pxx_dB)
theory_line = y0 + theory_slope * np.log10(f/f0)

# ==============================
# Plot
# ==============================
plt.figure(figsize=(10,6))
plt.semilogx(f, Pxx_dB, label=f"{N_order}-th order PSD", color="blue")
plt.semilogx(f, theory_line, 'k--', label=f"Theory slope {theory_slope} dB/dec (n={N_order})")

plt.xlabel("Normalized Frequency")
plt.ylabel("PSD [dB]")
plt.title(f"Noise Shaping of {N_order}-th Order ΔΣ Modulator")
plt.grid(True, which="both")
plt.legend()
plt.show()