現(xiàn)代我們經(jīng)常使用的域就是時(shí)域和頻域，那就出現(xiàn)了兩個(gè)技術(shù)流派，它們都能讓 ADC 獲得“更高帶寬 / 采樣速率”，但原理完全相反； “時(shí)間交織 (Time Interleaving)” 和 “頻譜分片 (Spectral Slicing)” 這兩種完全不同的架構(gòu)理念。

傳統(tǒng)電子 ADC：時(shí)域拼接（Time-Domain Interleaving）

把時(shí)間軸切成許多小片，每個(gè)子 ADC 輪流采樣不同的時(shí)間點(diǎn)。

例如：4 個(gè) ADC，每個(gè)采樣率 2 GS/s，時(shí)序錯(cuò)開 90°，→ 總體等效采樣率 = 8 GS/s。

每個(gè)子通道在不同時(shí)間點(diǎn)采樣；最后再把各自的采樣點(diǎn)按時(shí)間順序拼接回完整時(shí)域波形；所以這是 “時(shí)域拼接”。主要問題是時(shí)間偏差（skew），增益、帶寬不匹配→ 高頻時(shí)會(huì)出現(xiàn)“交織雜散 spur”。

光電混合 ADC（例如那篇 Kerr 微梳 ADC）：頻域拼接（Frequency-Domain Slicing）

把頻譜分成多個(gè)頻段，每個(gè)頻段用一個(gè)光載波“并行采樣”，最后在頻域合成。

也就是：同一時(shí)間段內(nèi)，信號(hào)被“分頻”送入不同子系統(tǒng)采集，各通道采的不是時(shí)間片，而是頻率片。

輸入信號(hào) → 調(diào)制到光載波上；通過 光學(xué)濾波器 (ISP1) 分成多個(gè)相鄰頻帶；每個(gè)頻帶由獨(dú)立的 IQ 相干接收機(jī) 檢測；每個(gè)接收機(jī)輸出低速電子信號(hào)；最后經(jīng) 數(shù)字頻譜拼接 (Spectral Stitching) 合成完整 0–320 GHz 波形。其中所有通道同時(shí)采樣同一時(shí)間段；各通道分到不同頻率子帶；在數(shù)字域把頻譜再拼接回一條完整頻帶；所以這是 “頻域拼接”。

很明顯不受電子采樣時(shí)鐘抖動(dòng)限制；可同時(shí)覆蓋數(shù)百 GHz 帶寬；頻域拼接可以精確標(biāo)定重疊區(qū)域，相位鎖定精度高。

展示如何頻域拼接

S_in(f) ──? H_eo(f) ──? 分片窗 slice_i(f)
                     └─? 各片 ×(未知復(fù)系數(shù) c_i) + 噪聲
                             └─? 通過 overlap 區(qū)比對(duì)恢復(fù) c_i
                                     └─? Σ_i slice_i × ?_i / Σ_i window_i  → 拼接譜
                                             └─? ×H_eo?1(f) (數(shù)字補(bǔ)償)
                                                     └─? IFFT 得到 y(t)

用 Python 仿真該系統(tǒng)的頻率響應(yīng)與譜片拼接過程（spectral-sliced ADC 重建），展示 320 GHz 等效帶寬 ADC 的概念,代碼我也給出來。

Fs = 640.0 GSa/s, N = 32768, df = 19.531 MHz, f_Nyquist = 320.0 GHz

True per-slice complex factors vs. estimated (relative):
Slice 0: true= 1.000 ∠ 0.0° | est= 1.000 ∠ 0.0°
Slice 1: true= 0.941 ∠ 37.5° | est= 0.730 ∠ -83.1°
Slice 2: true= 0.923 ∠ 157.0° | est= 1.483 ∠ -170.7°
Slice 3: true= 1.142 ∠ -7.6° | est= 1.196 ∠ -5.6°

低頻 ~0 dB，100 GHz 起平滑滾降（3 dB@100 GHz 的量級(jí)），290 GHz 處有明顯“陷波”（模擬論文里提到的 dip）， → 這決定了未經(jīng)補(bǔ)償?shù)南到y(tǒng)在高頻與 290 GHz 附近會(huì)被抑制。

每個(gè)譜片寬約 2×FSR≈80 GHz，相鄰譜片用 ≈5 GHz 余量 做余弦過渡，出現(xiàn)小重疊區(qū)，用于后續(xù)幅相對(duì)齊與譜片拼接。

原始輸入頻譜 vs 經(jīng)過 MZM 的頻譜

構(gòu)造了三段帶通信號(hào)：24.4 GHz（位于 0–80 GHz 范圍）、233.4 GHz（位于 160–240 GHz）、264.4 GHz（位于 240–320 GHz）；經(jīng)過 MZM 后，高頻端被壓低，290 GHz 附近凹陷明顯。

觀測到的各譜片頻譜（含噪）

對(duì)每片引入了未知復(fù)數(shù)增益/相位（模擬 LO 漂移/光纖相位），并疊加了隨頻率上升的噪聲（模擬 OCNR/ASE）；這些“未知量”隨后通過重疊區(qū)比對(duì)（中值估計(jì)）自動(dòng)恢復(fù)。

頻響平坦度（未補(bǔ)償）：范圍 54.89 dB
頻響平坦度（補(bǔ)償后）：范圍 52.58 dB

SINAD / ENOB 估計(jì)（未補(bǔ)償）：
Tone @ 2.0 GHz: SINAD = -47.6 dB, ENOB = -8.20 bits
Tone @ 56.0 GHz: SINAD = -80.6 dB, ENOB = -13.68 bits
Tone @ 280.0 GHz: SINAD = -199.8 dB, ENOB = -33.48 bits
Tone @ 308.0 GHz: SINAD = -87.7 dB, ENOB = -14.86 bits

拼接 + 數(shù)字補(bǔ)償后頻譜

對(duì)拼接結(jié)果乘以 H_eo 的正則化逆（含小正則項(xiàng)防放大噪聲過頭）；觀感上整體更平，但陷波附近的噪聲底會(huì)被抬升（符合論文中“數(shù)字補(bǔ)償會(huì)把噪聲一起放大”的現(xiàn)象）。

時(shí)域片段（原始 vs 拼接（未補(bǔ)償） vs 拼接+補(bǔ)償）

未補(bǔ)償版本的包絡(luò)被高頻抑制；補(bǔ)償后回到更接近原始，這對(duì)應(yīng)“通過頻域補(bǔ)償拉平頻響 → 時(shí)域恢復(fù)波形細(xì)節(jié)”。

對(duì)比

左圖：時(shí)域交織采樣（Time-Interleaved Sampling）

黑虛線：原始模擬信號(hào)（一個(gè)高頻包絡(luò)脈沖）。

彩色圓點(diǎn)：4 個(gè)子通道在不同時(shí)間相位錯(cuò)開的采樣點(diǎn)。

各通道輪流采樣時(shí)間軸，再拼接成完整波形，表示時(shí)間分時(shí)采樣，能直接觀測瞬態(tài)信號(hào)變化。

右圖：頻域分片采樣（Spectral Slicing & Stitching）

黑線：信號(hào)的真實(shí)頻譜。

彩色帶：不同通道負(fù)責(zé)的頻率子帶（譜片 0–3）。

每個(gè)通道獨(dú)立測量該頻帶內(nèi)容，最后再數(shù)字拼接（Stitching），表示頻率并行采樣，主要獲取穩(wěn)態(tài)頻譜信息，時(shí)域波形需要通過 IFFT 后處理重建，不能實(shí)時(shí)觀察瞬態(tài)。

兩者互為傅里葉對(duì)偶，正好體現(xiàn)了采樣定理的時(shí)間–頻率對(duì)稱性。

代碼

# -*- coding: utf-8 -*-
"""
Spectrally-sliced photonic-electronic ADC concept demo
- 320 GHz acquisition bandwidth (effective Fs = 640 GSa/s)
- Frequency slicing, overlap stitching, and optional EO transfer compensation
Notes:
- Uses frequency-domain synthesis with random complex content in 3 bands
- Simulates per-slice unknown complex gains/phases and recovers them via overlap stitching
- Demonstrates MZM roll-off and a notch near 290 GHz, and optional digital compensation
"""
import numpy as np
import matplotlib.pyplot as plt

# ----------------------------
# 1) Global simulation params
# ----------------------------
Fs = 640e9                 # effective sampling rate (Hz) -> Nyquist 320 GHz
N  = 2**15                 # FFT size (~1e4 pts); adjust for speed vs resolution
df = Fs/N
f = np.fft.rfftfreq(N, d=1/Fs)       # one-sided frequency grid (0..Fs/2)
fmax = f[-1]

# Sanity
print("Fs = %.1f GSa/s, N = %d, df = %.3f MHz, f_Nyquist = %.1f GHz" % (Fs/1e9, N, df/1e6, fmax/1e9))

# ----------------------------
# 2) Define input spectrum: 3 bandlimited regions (like paper)
# ----------------------------
rng = np.random.default_rng(42)

S_in = np.zeros_like(f, dtype=complex)

def add_band(center_GHz, baud_GBd, span_GHz, amp=1.0):
    # Fill |f-center| <= span/2 with random complex values shaped by raised-cosine like window
    center = center_GHz*1e9
    span = span_GHz*1e9
    idx = np.where(np.abs(f-center) <= span/2)[0]
    if len(idx)==0:
        return
    # random complex payload
    x = rng.standard_normal(len(idx)) + 1j*rng.standard_normal(len(idx))
    # smooth edges with cosine roll-off
    xw = x.copy()
    edges = np.linspace(-1, 1, len(idx))
    win = 0.5*(1+np.cos(np.pi*edges))  # raised cosine across the band to avoid sharp edges
    xw *= win
    S_in[idx] += amp * xw

# Three example signals (roughly mirroring the paper’s bands)
add_band(center_GHz=24.4,  baud_GBd=30, span_GHz=30)    # around 0-80 GHz band
add_band(center_GHz=233.4, baud_GBd=40, span_GHz=40)    # in 160-240 GHz band
add_band(center_GHz=264.4, baud_GBd=10, span_GHz=15)    # in 240-320 GHz band

# ----------------------------
# 3) EO modulator transfer function (magnitude-only model)
#    - gentle LP roll-off (3 dB @ 100 GHz, ~10 dB @ 300 GHz)
#    - notch near 290 GHz
# ----------------------------
def H_eo_mag(f_Hz):
    fG = f_Hz/1e9
    # base low-pass envelope (smooth knee around 100 GHz)
    lp = 1.0/np.sqrt(1 + (fG/100.0)**2.2)   # phenomenological slope
    # notch near 290 GHz (Gaussian)
    notch = 1 - 0.7*np.exp(-0.5*((fG-290)/6.0)**2)  # ~ -10 dB dip at center
    return lp*notch

H_eo = H_eo_mag(f)

# apply EO response to the optical modulation (here equivalent in our RF spectral surrogate)
S_after_eo = S_in * H_eo

# ----------------------------
# 4) Define spectral slices (M=4), each ~80 GHz wide, slight overlaps
# ----------------------------
FSR = 40.025e9
slice_bw = 2*FSR                # ~80.05 GHz per slice
M = 4

# Slice edges with ~5 GHz cosine overlaps
overlap = 5e9
slice_edges = [(k*slice_bw, (k+1)*slice_bw) for k in range(M)]
slice_windows = []

def cosine_ramp(x, x0, x1):
    # returns 0..1..0 raised-cosine style window across [x0,x1], 0 outside
    w = np.zeros_like(x, dtype=float)
    band = (x>=x0) & (x<=x1)
    if not np.any(band):
        return w
    xx = (x[band]-x0)/(x1-x0)  # 0..1
    w[band] = 0.5*(1 - np.cos(2*np.pi*xx))
    return w

# Build overlapping windows with flat center and cosine tapers of 'overlap' on both sides
for (f0, f1) in slice_edges:
    # extended edges for overlap
    a = max(0.0, f0 - overlap)
    b = min(Fs/2, f1 + overlap)
    # construct trapezoid with cosine skirts: 0->cosine->1(flat)->cosine->0
    w = np.zeros_like(f, dtype=float)
    # rise edge
    rise_end = f0
    rise_start = a
    rise = (f>=rise_start) & (f<=rise_end)
    if np.any(rise) and rise_end>rise_start:
        xx = (f[rise]-rise_start)/(rise_end-rise_start)  # 0..1
        w[rise] = 0.5*(1 - np.cos(np.pi*xx))
    # flat
    flat = (f>f0) & (f<f1)
    w[flat] = 1.0
    # fall edge
    fall_start = f1
    fall_end = b
    fall = (f>=fall_start) & (f<=fall_end)
    if np.any(fall) and fall_end>fall_start:
        xx = (f[fall]-fall_start)/(fall_end-fall_start)  # 0..1
        w[fall] = 0.5*(1 + np.cos(np.pi*xx))
    slice_windows.append(w)

# ----------------------------
# 5) Simulate unknown per-slice complex gains/phases (e.g., LO drift, fiber phase)
# ----------------------------
true_slice_cplx = []
for k in range(M):
    gain = 0.8 + 0.4*rng.random()               # 0.8..1.2
    phase = rng.uniform(-np.pi, np.pi)
    true_slice_cplx.append(gain*np.exp(1j*phase))
true_slice_cplx = np.array(true_slice_cplx, dtype=complex)

# Per-slice observed spectra (with EO response and slice windows + unknown complex factor)
Y_slices = []
for k in range(M):
    Yk = S_after_eo * slice_windows[k] * true_slice_cplx[k]
    # add AWGN to emulate ASE/OCNR; larger at higher f (simple model)
    noise_mag = 1e-3 * (1 + (f/1e11))   # grows with frequency
    noise = (rng.standard_normal(len(f)) + 1j*rng.standard_normal(len(f))) * noise_mag / np.sqrt(2)
    Yk = Yk + noise
    Y_slices.append(Yk)

# ----------------------------
# 6) "Measure" OE transfer per slice (here we assume known slice windows; unknown complex scalar only)
#    Recover per-slice complex factors via overlap stitching (relative to slice 0)
# ----------------------------
# Find overlaps between consecutive slices and compute complex scale that best matches in LS sense
est_slice_cplx = np.ones(M, dtype=complex)
for k in range(1, M):
    # overlapping region between slice k-1 and k: where both windows > threshold
    ov = (slice_windows[k-1] > 0.2) & (slice_windows[k] > 0.2)
    # To avoid divide-by-zero, select bins with sufficient power
    mask = ov & (np.abs(Y_slices[k])>1e-6) & (np.abs(Y_slices[k-1])>1e-6)
    if np.sum(mask) < 10:
        est_slice_cplx[k] = 1.0 + 0j
    else:
        # Estimate ratio that makes Y_{k} align to Y_{k-1}
        r = Y_slices[k-1][mask] / Y_slices[k][mask]
        # robust estimate: median in complex plane (use mean of normalized ratios)
        est = np.median(r)
        est_slice_cplx[k] = est * est_slice_cplx[k-1]  # chain relative to slice 0

# Normalize so that slice 0 factor = 1
est_slice_cplx = est_slice_cplx / est_slice_cplx[0]

print("\nTrue per-slice complex factors vs. estimated (relative):")
for k in range(M):
    print("Slice %d: true= %.3f ∠ %.1f° | est= %.3f ∠ %.1f°" % (
        k, np.abs(true_slice_cplx[k])/np.abs(true_slice_cplx[0]), 
        (np.angle(true_slice_cplx[k])-np.angle(true_slice_cplx[0]))*180/np.pi,
        np.abs(est_slice_cplx[k]), np.angle(est_slice_cplx[k])*180/np.pi))

# ----------------------------
# 7) Stitch: apply estimated factors and overlap-add (sum of weighted slices)
# ----------------------------
Y_corr = []
for k in range(M):
    Y_corr.append(Y_slices[k] * est_slice_cplx[k])

# Weighted sum with window-power normalization to avoid gain bumps in overlaps
Wsum = np.sum(slice_windows, axis=0) + 1e-12
Y_stitched = np.sum(Y_corr, axis=0) / Wsum

# Optional: digital compensation for EO roll-off (regularized inverse)
reg = 1e-2
H_inv = np.conj(H_eo) / (H_eo**2 + reg)
Y_comp = Y_stitched * H_inv

# ----------------------------
# 8) Compare spectra (magnitude)
# ----------------------------
def plot_mag_spectrum(freq, X, title):
    plt.figure(figsize=(8,4))
    plt.plot(freq/1e9, 20*np.log10(np.abs(X)+1e-15))
    plt.xlabel("Frequency (GHz)")
    plt.ylabel("Magnitude (dB)")
    plt.title(title)
    plt.grid(True)
    plt.tight_layout()
    plt.show()

# Plot EO transfer
plt.figure(figsize=(8,4))
plt.plot(f/1e9, 20*np.log10(H_eo+1e-15))
plt.xlabel("Frequency (GHz)")
plt.ylabel("EO |H_eo| (dB)")
plt.title("MZM EO 傳遞函數(shù)幅度（含 290 GHz 陷波與高頻滾降）")
plt.grid(True)
plt.tight_layout()
plt.show()

# Plot slice windows
for k in range(M):
    plt.figure(figsize=(8,3))
    plt.plot(f/1e9, slice_windows[k])
    plt.xlabel("Frequency (GHz)")
    plt.ylabel("Slice window")
    plt.title(f"譜片窗口 Slice {k}  (寬≈{slice_bw/1e9:.1f} GHz)")
    plt.grid(True)
    plt.tight_layout()
    plt.show()

# Original vs after EO
plot_mag_spectrum(f, S_in, "原始輸入頻譜 |S_in(f)|")
plot_mag_spectrum(f, S_after_eo, "經(jīng) MZM 后的頻譜 |S_in(f)·H_eo(f)|")

# Per-slice observed (magnitude)
for k in range(M):
    plot_mag_spectrum(f, Y_slices[k], f"觀測譜片 Slice {k}（含隨機(jī)幅相與噪聲）")

# Stitched spectra
plot_mag_spectrum(f, Y_stitched, "拼接后頻譜 |Y_stitched(f)|（未做 EO 補(bǔ)償）")
plot_mag_spectrum(f, Y_comp, "拼接+數(shù)字補(bǔ)償后頻譜 |Y_comp(f)|")

# ----------------------------
# 9) Time-domain reconstruction and quick SINAD test with single tone sweep
# ----------------------------
# Build time-domain from stitched (without and with compensation)
# Complete Hermitian spectrum for IFFT: we used rfft frequencies, so use irfft
y_time = np.fft.irfft(Y_stitched, n=N)
y_comp_time = np.fft.irfft(Y_comp, n=N)
x_time = np.fft.irfft(S_in, n=N)  # ideal reference (pre-EO)

# Simple alignment (normalize power)
def rms(a): return np.sqrt(np.mean(np.abs(a)**2))
scale = rms(x_time)/max(rms(y_time),1e-12)
scale_c = rms(x_time)/max(rms(y_comp_time),1e-12)

# Plot short time segments
def plot_time(t_ns, sig, title):
    plt.figure(figsize=(8,3))
    plt.plot(t_ns, sig)
    plt.xlabel("Time (ns)")
    plt.ylabel("Amplitude (a.u.)")
    plt.title(title)
    plt.grid(True)
    plt.tight_layout()
    plt.show()

t = np.arange(N)/Fs
sel = slice(0, 4000)  # first 4000 samples for visibility
plot_time(t[sel]*1e9, x_time[sel], "時(shí)域：原始信號(hào) 片段")
plot_time(t[sel]*1e9, (scale*y_time)[sel], "時(shí)域：拼接后（未補(bǔ)償） 片段")
plot_time(t[sel]*1e9, (scale_c*y_comp_time)[sel], "時(shí)域：拼接+補(bǔ)償 片段")

# ----------------------------
# 10) Frequency response via single-tone sweep (measure amplitude error)
# ----------------------------
tones_GHz = np.linspace(5, 315, 40)  # 40 tones spanning band
meas_gain_no_comp = []
meas_gain_comp = []

for ft in tones_GHz*1e9:
    S_t = np.zeros_like(f, dtype=complex)
    # narrow tone represented by a small bin (nearest bin)
    kbin = np.argmin(np.abs(f-ft))
    S_t[kbin] = 1.0
    # EO + slicing + unknown slice factors + stitching + optional comp
    Yt_slices = []
    for k in range(M):
        Yk = S_t * H_eo * slice_windows[k] * true_slice_cplx[k]
        # add small noise
        noise = (rng.standard_normal(len(f)) + 1j*rng.standard_normal(len(f))) * 1e-4 / np.sqrt(2)
        Yk = Yk + noise
        Yt_slices.append(Yk)
    # estimate factors using same overlap estimator as above (recompute for fairness)
    est_c = np.ones(M, dtype=complex)
    for k in range(1, M):
        ov = (slice_windows[k-1] > 0.2) & (slice_windows[k] > 0.2)
        mask = ov & (np.abs(Yt_slices[k])>1e-9) & (np.abs(Yt_slices[k-1])>1e-9)
        if np.sum(mask) < 1:
            est = 1+0j
        else:
            r = Yt_slices[k-1][mask] / Yt_slices[k][mask]
            est = np.median(r)
        est_c[k] = est * est_c[k-1]
    est_c = est_c/est_c[0]

    Wsum = np.sum(slice_windows, axis=0) + 1e-12
    Yt_st = np.sum([Yt_slices[k]*est_c[k] for k in range(M)], axis=0)/Wsum
    mag_no = np.abs(Yt_st[kbin])
    # compensation
    Yt_comp = Yt_st * H_inv
    mag_co = np.abs(Yt_comp[kbin])
    meas_gain_no_comp.append(mag_no)
    meas_gain_comp.append(mag_co)

meas_gain_no_comp = np.array(meas_gain_no_comp)
meas_gain_comp = np.array(meas_gain_comp)

# Normalize to low-frequency median to see flatness
g0 = np.median(meas_gain_no_comp[0:5])
g0c = np.median(meas_gain_comp[0:5])
flat_no = 20*np.log10(meas_gain_no_comp/g0 + 1e-15)
flat_co = 20*np.log10(meas_gain_comp/g0c + 1e-15)

# Print quick flatness stats
print("\n頻響平坦度（未補(bǔ)償）：范圍 %.2f dB" % (np.max(flat_no)-np.min(flat_no)))
print("頻響平坦度（補(bǔ)償后）：范圍 %.2f dB" % (np.max(flat_co)-np.min(flat_co)))

# Plot frequency response curves
plt.figure(figsize=(8,4))
plt.plot(tones_GHz, flat_no, marker='o')
plt.xlabel("Frequency (GHz)")
plt.ylabel("Relative gain (dB)")
plt.title("頻響：拼接后（未補(bǔ)償）")
plt.grid(True)
plt.tight_layout()
plt.show()

plt.figure(figsize=(8,4))
plt.plot(tones_GHz, flat_co, marker='o')
plt.xlabel("Frequency (GHz)")
plt.ylabel("Relative gain (dB)")
plt.title("頻響：拼接+數(shù)字補(bǔ)償后")
plt.grid(True)
plt.tight_layout()
plt.show()

# ----------------------------
# 11) Quick SINAD/ENOB at a few tones (post-stitched, no-comp)
# ----------------------------
def measure_sinad(ftone_Hz, snr_noise=1e-4):
    # synthesize tone in freq domain, run through pipeline, back to time, compute SINAD
    S_t = np.zeros_like(f, dtype=complex)
    kbin = np.argmin(np.abs(f-ftone_Hz))
    S_t[kbin] = 1.0
    Yt_slices = []
    for k in range(M):
        Yk = S_t * H_eo * slice_windows[k] * true_slice_cplx[k]
        noise = (rng.standard_normal(len(f)) + 1j*rng.standard_normal(len(f))) * snr_noise/np.sqrt(2)
        Yk = Yk + noise
        Yt_slices.append(Yk)
    # stitch with estimated factors
    est_c = np.ones(M, dtype=complex)
    for k in range(1, M):
        ov = (slice_windows[k-1] > 0.2) & (slice_windows[k] > 0.2)
        mask = ov & (np.abs(Yt_slices[k])>1e-9) & (np.abs(Yt_slices[k-1])>1e-9)
        if np.sum(mask) < 1:
            est = 1+0j
        else:
            r = Yt_slices[k-1][mask] / Yt_slices[k][mask]
            est = np.median(r)
        est_c[k] = est * est_c[k-1]
    est_c = est_c/est_c[0]
    Wsum = np.sum(slice_windows, axis=0) + 1e-12
    Yst = np.sum([Yt_slices[k]*est_c[k] for k in range(M)], axis=0)/Wsum
    y = np.fft.irfft(Yst, n=N)
    # compute SINAD in time domain: tone bin power / rest power
    # estimate tone amplitude by projecting onto sin/cos at ftone
    t = np.arange(N)/Fs
    s = np.sin(2*np.pi*ftone_Hz*t); c = np.cos(2*np.pi*ftone_Hz*t)
    A_s = 2*np.mean(y*s); A_c = 2*np.mean(y*c)
    tone_power = 0.5*(A_s**2 + A_c**2)
    total_power = np.mean(y**2)
    noise_dist_power = max(total_power - tone_power, 1e-20)
    sinad = tone_power/noise_dist_power
    sinad_db = 10*np.log10(sinad+1e-20)
    enob = (sinad_db - 1.76)/6.02
    return sinad_db, enob

test_tones = [2e9, 56e9, 280e9, 308e9]
print("\nSINAD / ENOB 估計(jì)（未補(bǔ)償）：")
for ft in test_tones:
    sdb, eb = measure_sinad(ft)
    print("Tone @ %.1f GHz: SINAD = %.1f dB, ENOB = %.2f bits" % (ft/1e9, sdb, eb))