アナログ信号 x(t)
  ↓ [アンチエイリアシングフィルタ] (f_max 以上をカット)
  ↓ [標本化] (f_s > 2*f_max)
  x[n] = x(n*T_s)
  ↓ [量子化] (B ビット一様量子化)
  x_q[n] = Q(x[n])
  ↓ [符号化] (自然二進符号 or グレイ符号)
  ビット列 b_0, b_1, ..., b_{B-1}
  ↓ [送信]

import numpy as np
import matplotlib.pyplot as plt

# ===== サンプリングとエイリアシングの可視化 =====
# 元の信号: 5 Hz の正弦波
f_signal = 5  # 信号周波数 [Hz]
t_cont = np.linspace(0, 1, 10000)  # 連続時間（近似）
x_cont = np.sin(2 * np.pi * f_signal * t_cont)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# (1) 適切なサンプリング (fs = 40 Hz > 2*5 = 10 Hz)
fs1 = 40
ts1 = np.arange(0, 1, 1/fs1)
xs1 = np.sin(2 * np.pi * f_signal * ts1)
ax = axes[0, 0]
ax.plot(t_cont, x_cont, 'b-', alpha=0.5, label=f'Original ({f_signal} Hz)')
ax.stem(ts1, xs1, linefmt='r-', markerfmt='ro', basefmt='k-', label=f'Sampled ($f_s$={fs1} Hz)')
ax.set_title(f'Proper Sampling ($f_s$={fs1} Hz > $2f_{{max}}$={2*f_signal} Hz)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Amplitude')

# (2) ナイキストぎりぎり (fs = 12 Hz)
fs2 = 12
ts2 = np.arange(0, 1, 1/fs2)
xs2 = np.sin(2 * np.pi * f_signal * ts2)
ax = axes[0, 1]
ax.plot(t_cont, x_cont, 'b-', alpha=0.5, label=f'Original ({f_signal} Hz)')
ax.stem(ts2, xs2, linefmt='r-', markerfmt='ro', basefmt='k-', label=f'Sampled ($f_s$={fs2} Hz)')
ax.set_title(f'Marginal Sampling ($f_s$={fs2} Hz)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Amplitude')

# (3) エイリアシング発生 (fs = 8 Hz < 2*5 = 10 Hz)
fs3 = 8
ts3 = np.arange(0, 1, 1/fs3)
xs3 = np.sin(2 * np.pi * f_signal * ts3)
# エイリアス周波数: f_alias = f_s - f_signal = 3 Hz
f_alias = fs3 - f_signal
x_alias = np.sin(2 * np.pi * f_alias * t_cont)
ax = axes[1, 0]
ax.plot(t_cont, x_cont, 'b-', alpha=0.5, label=f'Original ({f_signal} Hz)')
ax.plot(t_cont, x_alias, 'g--', alpha=0.7, label=f'Alias ({f_alias} Hz)')
ax.stem(ts3, xs3, linefmt='r-', markerfmt='ro', basefmt='k-', label=f'Sampled ($f_s$={fs3} Hz)')
ax.set_title(f'Aliasing ($f_s$={fs3} Hz < $2f_{{max}}$={2*f_signal} Hz)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Amplitude')

# (4) スペクトルの比較
ax = axes[1, 1]
# 各サンプリングレートでのスペクトル（レプリカ）を描画
for fs, color, label in [(40, 'blue', '$f_s$=40Hz'),
                          (12, 'green', '$f_s$=12Hz'),
                          (8, 'red', '$f_s$=8Hz')]:
    freqs = np.linspace(-60, 60, 1000)
    spectrum = np.zeros_like(freqs)
    for k in range(-3, 4):
        center = k * fs
        # 三角形で信号スペクトルを近似表現
        spectrum += np.maximum(0, 1 - np.abs(freqs - center) / f_signal)
        spectrum += np.maximum(0, 1 - np.abs(freqs + center) / f_signal)
    ax.plot(freqs, spectrum, color=color, alpha=0.6, label=label)

ax.set_title('Spectrum Replicas')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Amplitude')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim([-30, 30])

plt.tight_layout()
plt.savefig('sampling_aliasing.png', dpi=150, bbox_inches='tight')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# ===== sinc補間による信号復元 =====
f_signal = 5  # 信号周波数 [Hz]
fs = 20       # サンプリング周波数 [Hz]
Ts = 1 / fs

# 元の連続信号
t_cont = np.linspace(0, 1, 10000)
x_cont = np.sin(2 * np.pi * f_signal * t_cont)

# サンプリング
n_samples = int(1 * fs)
t_samples = np.arange(n_samples) * Ts
x_samples = np.sin(2 * np.pi * f_signal * t_samples)

# sinc補間による復元
def sinc_interpolation(t, t_samples, x_samples, Ts):
    """ホイッタカー-シャノンの補間公式"""
    x_reconstructed = np.zeros_like(t)
    for n in range(len(t_samples)):
        x_reconstructed += x_samples[n] * np.sinc((t - t_samples[n]) / Ts)
    return x_reconstructed

x_reconstructed = sinc_interpolation(t_cont, t_samples, x_samples, Ts)

# プロット
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# 上段: 復元結果
ax = axes[0]
ax.plot(t_cont, x_cont, 'b-', linewidth=2, label='Original signal', alpha=0.7)
ax.plot(t_cont, x_reconstructed, 'r--', linewidth=2, label='Reconstructed (sinc)')
ax.stem(t_samples, x_samples, linefmt='g-', markerfmt='go', basefmt='k-',
        label='Samples')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Amplitude')
ax.set_title(f'Sinc Interpolation ($f_s$={fs} Hz, $f_{{signal}}$={f_signal} Hz)')
ax.legend()
ax.grid(True, alpha=0.3)

# 下段: 各sinc関数の寄与
ax = axes[1]
for n in range(len(t_samples)):
    sinc_contribution = x_samples[n] * np.sinc((t_cont - t_samples[n]) / Ts)
    ax.plot(t_cont, sinc_contribution, alpha=0.4, linewidth=0.8)
ax.plot(t_cont, x_reconstructed, 'k-', linewidth=2, label='Sum of sinc functions')
ax.stem(t_samples, x_samples, linefmt='r-', markerfmt='ro', basefmt='k-')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Amplitude')
ax.set_title('Individual sinc contributions')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('sinc_interpolation.png', dpi=150, bbox_inches='tight')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# ===== 量子化シミュレーション =====
# 1kHz正弦波を量子化
f0 = 1000  # 信号周波数 [Hz]
fs = 44100  # サンプリング周波数 [Hz]
duration = 0.01  # 表示期間 [s]
t = np.arange(0, duration, 1/fs)
x = np.sin(2 * np.pi * f0 * t)

# 量子化関数
def uniform_quantize(x, B, V_max=1.0):
    """一様量子化"""
    L = 2**B
    Delta = 2 * V_max / L
    # 量子化
    x_clipped = np.clip(x, -V_max, V_max - Delta)
    x_q = Delta * np.floor(x_clipped / Delta + 0.5)
    return x_q, Delta

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

bit_depths = [2, 4, 8, 16]
for idx, B in enumerate(bit_depths):
    ax = axes[idx // 2, idx % 2]
    x_q, Delta = uniform_quantize(x, B)
    error = x - x_q

    ax.plot(t * 1000, x, 'b-', linewidth=1.5, label='Original', alpha=0.7)
    ax.step(t * 1000, x_q, 'r-', linewidth=1.5, label=f'Quantized ({B} bit)',
            where='mid')
    ax.set_xlabel('Time [ms]')
    ax.set_ylabel('Amplitude')
    ax.set_title(f'{B}-bit Quantization ($\\Delta$={Delta:.4f}, L={2**B})')
    ax.legend()
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('quantization_demo.png', dpi=150, bbox_inches='tight')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# ===== 量子化SNRの理論値と実測値の比較 =====
f0 = 1000
fs = 44100
duration = 1.0  # 1秒間
t = np.arange(0, duration, 1/fs)
x = np.sin(2 * np.pi * f0 * t)

def uniform_quantize(x, B, V_max=1.0):
    """一様量子化"""
    L = 2**B
    Delta = 2 * V_max / L
    x_clipped = np.clip(x, -V_max, V_max - Delta)
    x_q = Delta * np.floor(x_clipped / Delta + 0.5)
    return x_q, Delta

bit_range = np.arange(2, 25)
snr_theory = 6.02 * bit_range + 1.76  # 理論値 [dB]
snr_measured = []

for B in bit_range:
    x_q, _ = uniform_quantize(x, B)
    error = x - x_q
    P_signal = np.mean(x**2)
    P_noise = np.mean(error**2)
    if P_noise > 0:
        snr_measured.append(10 * np.log10(P_signal / P_noise))
    else:
        snr_measured.append(np.inf)

snr_measured = np.array(snr_measured)

plt.figure(figsize=(10, 6))
plt.plot(bit_range, snr_theory, 'b-', linewidth=2, label='Theory: $6.02B + 1.76$ dB')
plt.plot(bit_range, snr_measured, 'ro', markersize=6, label='Measured (sine wave)')
plt.xlabel('Number of Bits $B$', fontsize=12)
plt.ylabel('SNR [dB]', fontsize=12)
plt.title('Quantization SNR vs Number of Bits', fontsize=13)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('quantization_snr.png', dpi=150, bbox_inches='tight')
plt.show()

for B in [8, 12, 16, 24]:
    idx = B - 2
    print(f"B={B:2d}: Theory={snr_theory[idx]:.1f} dB, Measured={snr_measured[idx]:.1f} dB")

import numpy as np
import matplotlib.pyplot as plt

# ===== μ-law 圧伸の特性と効果 =====
def mu_law_compress(x, mu=255):
    """μ-law 圧縮"""
    return np.sign(x) * np.log(1 + mu * np.abs(x)) / np.log(1 + mu)

def mu_law_expand(y, mu=255):
    """μ-law 伸張（逆変換）"""
    return np.sign(y) * (1/mu) * ((1 + mu)**np.abs(y) - 1)

def a_law_compress(x, A=87.6):
    """A-law 圧縮"""
    y = np.zeros_like(x)
    mask_small = np.abs(x) <= 1/A
    mask_large = np.abs(x) > 1/A
    y[mask_small] = A * np.abs(x[mask_small]) / (1 + np.log(A)) * np.sign(x[mask_small])
    y[mask_large] = (1 + np.log(A * np.abs(x[mask_large]))) / (1 + np.log(A)) * np.sign(x[mask_large])
    return y

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# (1) μ-law 圧縮特性
ax = axes[0, 0]
x_range = np.linspace(-1, 1, 1000)
for mu_val in [0, 5, 50, 255]:
    if mu_val == 0:
        y = x_range  # 線形（一様量子化）
        label = '$\\mu$=0 (linear)'
    else:
        y = mu_law_compress(x_range, mu_val)
        label = f'$\\mu$={mu_val}'
    ax.plot(x_range, y, linewidth=2, label=label)
ax.set_xlabel('Input $x$')
ax.set_ylabel('Output $F_\\mu(x)$')
ax.set_title('$\\mu$-law Compression Characteristics')
ax.legend()
ax.grid(True, alpha=0.3)

# (2) A-law 圧縮特性
ax = axes[0, 1]
ax.plot(x_range, x_range, 'b-', linewidth=2, label='Linear')
ax.plot(x_range, a_law_compress(x_range), 'r-', linewidth=2, label='A-law (A=87.6)')
ax.plot(x_range, mu_law_compress(x_range), 'g--', linewidth=2, label='$\\mu$-law ($\\mu$=255)')
ax.set_xlabel('Input $x$')
ax.set_ylabel('Compressed Output')
ax.set_title('Comparison of Companding Laws')
ax.legend()
ax.grid(True, alpha=0.3)

# (3) 量子化SNRの比較（一様 vs μ-law）
ax = axes[1, 0]
B = 8  # 8ビット量子化
input_levels_dB = np.linspace(-50, 0, 50)  # 入力レベル [dB]

snr_uniform = []
snr_mulaw = []

for level_dB in input_levels_dB:
    amplitude = 10**(level_dB / 20)
    # テスト信号
    t = np.linspace(0, 0.1, 10000)
    x = amplitude * np.sin(2 * np.pi * 100 * t)

    # 一様量子化
    L = 2**B
    Delta = 2.0 / L
    x_q_uniform = Delta * np.floor(np.clip(x, -1, 1-Delta) / Delta + 0.5)
    e_uniform = x - x_q_uniform
    P_s = np.mean(x**2)
    P_e_uniform = np.mean(e_uniform**2)
    snr_uniform.append(10*np.log10(P_s / P_e_uniform) if P_e_uniform > 0 else 0)

    # μ-law 圧伸量子化
    x_compressed = mu_law_compress(x, mu=255)
    x_compressed_q = Delta * np.floor(np.clip(x_compressed, -1, 1-Delta) / Delta + 0.5)
    x_expanded = mu_law_expand(x_compressed_q, mu=255)
    e_mulaw = x - x_expanded
    P_e_mulaw = np.mean(e_mulaw**2)
    snr_mulaw.append(10*np.log10(P_s / P_e_mulaw) if P_e_mulaw > 0 else 0)

ax.plot(input_levels_dB, snr_uniform, 'b-', linewidth=2, label='Uniform')
ax.plot(input_levels_dB, snr_mulaw, 'r-', linewidth=2, label='$\\mu$-law ($\\mu$=255)')
ax.set_xlabel('Input Level [dB]')
ax.set_ylabel('SNR [dB]')
ax.set_title(f'Quantization SNR vs Input Level ({B}-bit)')
ax.legend()
ax.grid(True, alpha=0.3)

# (4) PCM符号化の例
ax = axes[1, 1]
B_demo = 3
t_demo = np.linspace(0, 1, 1000)
x_demo = 0.8 * np.sin(2 * np.pi * 2 * t_demo)

L_demo = 2**B_demo
Delta_demo = 2.0 / L_demo
x_q_demo = Delta_demo * np.floor(np.clip(x_demo, -1, 1-Delta_demo) / Delta_demo + 0.5)

# サンプリング
fs_demo = 20
t_s_demo = np.arange(0, 1, 1/fs_demo)
x_s_demo = 0.8 * np.sin(2 * np.pi * 2 * t_s_demo)
x_sq_demo = Delta_demo * np.floor(np.clip(x_s_demo, -1, 1-Delta_demo) / Delta_demo + 0.5)

ax.plot(t_demo, x_demo, 'b-', alpha=0.5, label='Original')
ax.step(t_s_demo, x_sq_demo, 'r-', where='mid', linewidth=2, label=f'PCM ({B_demo}-bit)')
ax.stem(t_s_demo, x_sq_demo, linefmt='r-', markerfmt='ro', basefmt='k-')
# 量子化レベルを点線で表示
levels = np.arange(-1 + Delta_demo/2, 1, Delta_demo)
for lv in levels:
    ax.axhline(y=lv, color='gray', linestyle=':', alpha=0.3)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Amplitude')
ax.set_title(f'PCM Encoding ({B_demo}-bit, $f_s$={fs_demo} Hz)')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('pcm_companding.png', dpi=150, bbox_inches='tight')
plt.show()