为什么有斩波前端的 ADC 还是标有一个 1/f 的噪音值？（YUNSWJ 仿真版）

f = np.logspace(-3, 3, 1000)   # 1 mHz ~ 1 kHz
en_white = 10e-9               # 10 nV/√Hz 白噪声本底
fc = 1.0                        # 1 Hz 拐点

# 无 1/f：纯白噪声
en_flat = en_white * np.ones_like(f)

# 有 1/f：低频上升 ~ sqrt(fc/f)
en_1f = en_white * np.sqrt(1 + fc / f)

Fs = 1000.0       # 采样率 1 kS/s
T_total = 200.0   # 总时长 200 s
N = int(Fs * T_total)

np.random.seed(0)
x_white = np.random.normal(0, 1.0, N)   # 纯白噪声

Xw = np.fft.rfft(x_white)
freqs = np.fft.rfftfreq(N, d=1/Fs)
H_1f = np.sqrt(1 + fc / np.maximum(freqs, 1e-6))
H_1f[0] = 0.0  # DC 处避免无穷大
X_1f = Xw * H_1f
x_1f = np.fft.irfft(X_1f, n=N)

freqs_est, asd_white = estimate_asd(x_white, Fs)
_, asd_1f_est = estimate_asd(x_1f, Fs)

def rms_of_averaged(x, Fs, T_avg_list):
    rms_list = []
    for T in T_avg_list:
        L = int(T * Fs)                 # 每块样本数
        nseg = len(x) // L
        segs = x[:nseg*L].reshape(nseg, L)
        means = segs.mean(axis=1)       # 每块的平均值
        rms_list.append(np.sqrt(np.mean(means**2)))  # 这些平均值的 RMS
    return np.array(rms_list)

T_list = np.logspace(-1, 2, 20)  # 0.1s ~ 100s
rms_white = rms_of_averaged(x_white, Fs, T_list)
rms_1f = rms_of_averaged(x_1f, Fs, T_list)

import numpy as np
import matplotlib.pyplot as plt

# =========================
# 1. 频谱对比: 有 1/f vs 无 1/f
# =========================

# 频率轴: 1e-3 Hz ~ 10^3 Hz
f = np.logspace(-3, 3, 1000)
en_white = 10e-9   # 白噪声本底 10 nV/√Hz
fc = 1.0           # 1/f 拐点频率 1 Hz

# 无 1/f: 纯白噪声
en_flat = en_white * np.ones_like(f)

# 有 1/f: 低于 fc 时上升 ~ sqrt(fc/f)
en_1f = en_white * np.sqrt(1 + fc / f)

plt.figure()
plt.loglog(f, en_flat*1e9, label="无 1/f (斩波, 白噪声平坦)")
plt.loglog(f, en_1f*1e9, label="有 1/f (普通放大器)")
plt.axvline(fc, linestyle="--")
plt.xlabel("频率 / Hz")
plt.ylabel("电压噪声谱密度 / nV/√Hz")
plt.title("有 1/f 噪声 vs 无 1/f 噪声 的 ASD 对比")
plt.legend()
plt.grid(True, which="both")

# =========================
# 2. 时间域仿真: 平均 10s 时 RMS 噪声
# =========================

Fs = 1000.0          # 采样率 1 kS/s
T_total = 200.0      # 总时长 200 s
N = int(Fs * T_total)

# 生成一段白噪声
np.random.seed(0)
x_white = np.random.normal(0, 1.0, N)

# 通过频域塑形生成 "带 1/f" 噪声
Xw = np.fft.rfft(x_white)
freqs = np.fft.rfftfreq(N, d=1/Fs)
H_1f = np.sqrt(1 + fc / np.maximum(freqs, 1e-6))
H_1f[0] = 0.0  # DC 分量设 0，避免发散
X_1f = Xw * H_1f
x_1f = np.fft.irfft(X_1f, n=N)

# 验证频谱: 估算两种噪声的 ASD
def estimate_asd(x, Fs, nfft=16384):
    w = np.hanning(nfft)
    step = nfft // 2
    segs = (len(x) - nfft) // step
    psd = np.zeros(nfft//2 + 1)
    for i in range(segs):
        seg = x[i*step:i*step+nfft] * w
        X = np.fft.rfft(seg)
        psd += (np.abs(X)**2) / (Fs * np.sum(w**2))
    psd /= segs
    freqs = np.fft.rfftfreq(nfft, d=1/Fs)
    asd = np.sqrt(psd)
    return freqs, asd

freqs_est, asd_white = estimate_asd(x_white, Fs)
_, asd_1f_est = estimate_asd(x_1f, Fs)

plt.figure()
plt.loglog(freqs_est[1:], asd_white[1:], label="无 1/f (仿真)")
plt.loglog(freqs_est[1:], asd_1f_est[1:], label="有 1/f (仿真)")
plt.axvline(fc, linestyle="--")
plt.xlim(1e-1, 1e3)
plt.xlabel("频率 / Hz")
plt.ylabel("电压噪声谱密度 (任意单位)")
plt.title("仿真得到的噪声谱密度: 有 / 无 1/f")
plt.legend()
plt.grid(True, which="both")

# =========================
# 3. 统计平均时间 vs RMS 噪声
# =========================

def rms_of_averaged(x, Fs, T_avg_list):
    rms_list = []
    for T in T_avg_list:
        L = int(T * Fs)
        nseg = len(x) // L
        if nseg == 0:
            rms_list.append(np.nan)
            continue
        segs = x[:nseg*L].reshape(nseg, L)
        means = segs.mean(axis=1)
        rms_list.append(np.sqrt(np.mean(means**2)))
    return np.array(rms_list)

T_list = np.logspace(-1, 2, 20)  # 0.1 s ~ 100 s
rms_white = rms_of_averaged(x_white, Fs, T_list)
rms_1f = rms_of_averaged(x_1f, Fs, T_list)

plt.figure()
plt.loglog(T_list, rms_white, 'o-', label="无 1/f (白噪声)")
plt.loglog(T_list, rms_1f, 'o-', label="有 1/f")
plt.xlabel("平均时间 T / s")
plt.ylabel("平均值的 RMS 噪声 (任意单位)")
plt.title("平均时间 vs RMS 噪声：有 / 无 1/f")
plt.grid(True, which="both")
plt.legend()
plt.show()