从ADALM2000看网络分析仪原理

def acquire_from_m2k(f_sig, fs, n_cycles, vin_amp):
    # 1. 配置 AWG：频率 f_sig、幅度 vin_amp
    # 2. 配置 Scope：采样率 fs，通道1测 Vin，通道2测 Vout
    # 3. 等待系统稳定（可以丢掉前几十个周期）
    # 4. 抓取 N 点数据 → vin, vout
    return t, vin, vout

# 原来：
t, vin = generate_sine(...)
vout = simulate_dut(t, vin, f_sig, fc=fc_dut)

# 替换为：
t, vin, vout = acquire_from_m2k(f_sig, fs, n_cycles, vin_amp)

import numpy as np
import matplotlib.pyplot as plt

def generate_sine(f_sig, fs, n_cycles, amplitude=1.0):
    """
    生成某个频率下的正弦波测试信号。
    f_sig: 被测信号频率 (Hz)
    fs: 采样率 (Hz)
    n_cycles: 采样的周期数
    amplitude: 正弦幅度
    """
    T = n_cycles / f_sig           # 总时长
    N = int(np.round(T * fs))      # 采样点数
    t = np.arange(N) / fs
    x = amplitude * np.sin(2 * np.pi * f_sig * t)
    return t, x


def rc_lowpass_H(f, fc):
    """
    一阶 RC 低通的理论传递函数 H(j2πf) = 1 / (1 + j f/fc)
    f: 频率 (Hz) 或 numpy 数组
    fc: 截止频率 (Hz)
    """
    s = 1j * f / fc
    return 1.0 / (1.0 + s)


def simulate_dut(t, x, f_sig, fc=1e3):
    """
    用“理想一阶 RC 低通”来模拟 DUT。
    这边不做真正的滤波器卷积，而是直接用理论幅度/相位来构造 y(t)：
        y(t) = |H(f)| * sin(2π f t + φ)
    这样能保持算法结构，又不依赖 SciPy 等库。
    将来你用 M2K 时，把这一函数替换成“采集真实 y(t)”即可。
    """
    H = rc_lowpass_H(f_sig, fc)
    mag = np.abs(H)
    phase = np.angle(H)
    y = mag * np.sin(2 * np.pi * f_sig * t + phase)
    return y


def measure_transfer_fft(x, y, fs, f_sig):
    """
    用 FFT 在 f_sig 对应的频率 bin 上估计 H(f) = Y/X。
    x, y: 两路同步采样信号 (numpy 数组)
    fs: 采样率 (Hz)
    f_sig: 当前扫频点 (Hz)
    返回: H_est (复数)，包含幅度和相位
    """
    # 统一长度
    N = min(len(x), len(y))
    x = x[:N]
    y = y[:N]

    # 加窗，减小泄漏（这里用 Hann 窗）
    window = np.hanning(N)
    xw = x * window
    yw = y * window

    # 只做 rFFT（实信号），只关注正频
    X = np.fft.rfft(xw)
    Y = np.fft.rfft(yw)
    freqs = np.fft.rfftfreq(N, d=1.0/fs)

    # 找到最接近 f_sig 的 bin
    k = np.argmin(np.abs(freqs - f_sig))

    # 估计传递函数
    # 注意：真实系统里要处理 X[k] ≈ 0 的情况，这里简单略过
    H_est = Y[k] / X[k]
    return H_est


def sweep_network_analyzer(
    f_start=10.0,
    f_stop=1e5,
    points_per_decade=20,
    fs=1e6,
    n_cycles=200,
    vin_amp=1.0,
    fc_dut=1e3
):
    """
    扫频测量 DUT 的幅频 + 相频（这里 DUT 用一阶 RC 低通模拟）。
    将来如果接 M2K，只要改“获取 y(t)”那一步即可。

    参数含义：
    - f_start, f_stop: 扫频起止频率 (Hz)
    - points_per_decade: 每十倍频点的采样数（对数扫频）
    - fs: 采样率 (Hz)
    - n_cycles: 每个频点采多少个周期
    - vin_amp: 输入正弦幅度
    - fc_dut: 虚拟 RC 低通的截止频率 (Hz)
    """
    # 生成对数扫频频点
    n_decades = np.log10(f_stop) - np.log10(f_start)
    n_points = int(n_decades * points_per_decade)
    freqs = np.logspace(np.log10(f_start), np.log10(f_stop), n_points)

    mag_db_list = []
    phase_deg_list = []

    for f_sig in freqs:
        # 1) 生成输入正弦
        t, vin = generate_sine(f_sig, fs, n_cycles, amplitude=vin_amp)

        # 2) 通过 DUT（此处为仿真；真实 M2K 时用采集）
        vout = simulate_dut(t, vin, f_sig, fc=fc_dut)

        # 3) 用 FFT 算传递函数
        H_est = measure_transfer_fft(vin, vout, fs, f_sig)

        mag_db = 20 * np.log10(np.abs(H_est))
        phase_deg = np.angle(H_est, deg=True)

        mag_db_list.append(mag_db)
        phase_deg_list.append(phase_deg)

    mag_db_array = np.array(mag_db_list)
    phase_deg_array = np.unwrap(np.deg2rad(phase_deg_list))  # 弧度展开
    phase_deg_array = np.rad2deg(phase_deg_array)

    return freqs, mag_db_array, phase_deg_array


def plot_bode(freqs, mag_db, phase_deg, title="Simulated Network Analyzer"):
    """
    画一个双子图 Bode 图：上幅度，下相位。
    """
    fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(8, 6))

    ax1.semilogx(freqs, mag_db)
    ax1.set_ylabel("Magnitude (dB)")
    ax1.grid(True, which="both", ls=":")

    ax2.semilogx(freqs, phase_deg)
    ax2.set_ylabel("Phase (deg)")
    ax2.set_xlabel("Frequency (Hz)")
    ax2.grid(True, which="both", ls=":")

    fig.suptitle(title)
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    # 示例：被测 DUT 是 fc = 1 kHz 的一阶 RC 低通
    fc = 1e3

    freqs, mag_db, phase_deg = sweep_network_analyzer(
        f_start=10,
        f_stop=1e5,
        points_per_decade=20,
        fs=1e6,
        n_cycles=200,
        vin_amp=1.0,
        fc_dut=fc
    )
    plot_bode(freqs, mag_db, phase_deg,
              title=f"Simulated RC Lowpass (fc = {fc:.0f} Hz)")