ADC 位数在频谱分析中对高频的影响是巨大的（Σ-Δ 上大分）

# 设定输入噪声密度 50 nV/√Hz（随便举例）
e_n = 50e-9

# 4. 窗函数：Hann
window = np.hanning(N)

fs = 100_000.0      # 采样率 100 kS/s
N = 16384           # FFT 点数

B = 12              # 12 bit ADC

# 将噪声强度降低或者直接去除噪声
e_n = 10e-9  # 降低噪声密度
# 或者直接去掉噪声
# noise = np.zeros(N)
# x_analog += noise

# 尝试其他窗函数，比如 Hamming 或 Blackman
window = np.hamming(N)  # Hamming窗
# 或者使用Blackman窗
# window = np.blackman(N)

# 增加 FFT 点数 N 或者提高采样率 fs
fs = 200_000.0  # 提高采样率
N = 32768  # 增加 FFT 点数

# 通过增加 ADC 位数，减少量化误差
B = 16  # 使用 16 位 ADC

import numpy as np
import matplotlib.pyplot as plt

# 1. ADC / FFT 参数
fs = 1000_000.0      # 采样率 100 kS/s
N = 32768*2           # FFT 点数
B = 48              # 12 bit ADC
V_FS = 1.0          # 满量程 ±1 V，示例

# 2. 构造模拟输入信号：1 kHz 正弦 + 白噪声
f_sig = 1000.0
A_sig = 0.2         # 峰值 0.2 V → 约 0.141 Vrms

t = np.arange(N) / fs
x_analog = A_sig * np.sin(2 * np.pi * f_sig * t)

# 设定输入噪声密度 50 nV/√Hz（随便举例）
e_n = 10e-9
# 对应到每个采样点的时域噪声标准差：sigma = e_n * sqrt(fs/2)
# （粗略关系，准确推导更复杂，这里当示范）
sigma_noise = e_n * np.sqrt(fs/2)
noise = np.random.normal(scale=sigma_noise, size=N)

x_analog += noise

# 3. ADC 量化建模：12 bit，范围 ±V_FS
# 3.1 限幅
x_clip = np.clip(x_analog, -V_FS, V_FS)

# 3.2 量化步长
Delta = 2 * V_FS / (2**B)
codes = np.round(x_clip / Delta)   # 量化到整数码
codes = np.clip(codes, -(2**(B-1)), 2**(B-1)-1)

# 3.3 再转换回电压（这一步和 M2K 驱动做的类似）
x_digital = codes * Delta

# 4. 窗函数：Hann
window = np.hamming(N)
CG = np.sum(window) / N           # 相干增益
ENBW_factor = np.sum(window**2) * N / (np.sum(window)**2)
df = fs / N
ENBW = df * ENBW_factor

xw = x_digital * window

# 5. FFT（只看正频）
X = np.fft.rfft(xw)
freqs = np.fft.rfftfreq(N, d=1/fs)

# 6. 从 FFT 幅值恢复单频正弦的峰值、RMS、电压谱等

# 6.1 幅度谱（峰值，针对单频信号）
#     A_est[k] = 2 * |X[k]| / (N * CG)
A_est = 2 * np.abs(X) / (N * CG)
Vrms_line = A_est / np.sqrt(2)

# 6.2 把正弦的那个 bin 标出来看看
k_sig = np.argmin(np.abs(freqs - f_sig))
print("信号频率 bin:", k_sig, "对应频率:", freqs[k_sig], "Hz")
print("估计 Vrms:", Vrms_line[k_sig], "V")

# 7. 噪声谱：计算近似“电压噪声密度”（V/√Hz）
#    Vrms_line 是每个 bin 的等效 RMS（包含信号+噪声）。
#    对于噪声，可以除以 sqrt(ENBW) 得到近似噪声密度。
Vn_density = Vrms_line / np.sqrt(ENBW)

# 8. 转换为 dBV 与 dBV/√Hz
spec_dBV = 20 * np.log10(Vrms_line / 1.0 + 1e-20)
spec_dBV_Hz = 20 * np.log10(Vn_density / 1.0 + 1e-30)

# 9. 作图
plt.figure(figsize=(10, 6))
plt.semilogx(freqs, spec_dBV, label="Line spectrum (dBV)")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Amplitude (dBV)")
plt.grid(which="both", ls=":")
plt.legend()
plt.title("FFT Spectrum (windowed, Hann)")
plt.tight_layout()
plt.show()

plt.figure(figsize=(10, 6))
plt.semilogx(freqs, spec_dBV_Hz, label="Noise density (dBV/√Hz)")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Noise density (dBV/√Hz)")
plt.grid(which="both", ls=":")
plt.legend()
plt.title("Estimated Noise Spectral Density")
plt.tight_layout()
plt.show()