如何在Python中获取含噪正弦波的幅度、相位及直流偏移

原代码核心问题梳理 #

• numpy.fft.fftfreq参数错误：第一个参数需传入采样总点数，而非points_num-1，否则频率轴与FFT结果无法对应。
• 手动清零负频率的方式低效且易引入误差：直接通过索引截取正频率段即可，无需修改FFT结果数组。
• 未实现直流偏移计算：直流偏移可通过信号均值或FFT的0频分量直接获取。
• 相位计算未匹配信号模型：FFT的相位基准与自定义cosfunc的相位定义存在差异，需进行转换修正。
• 拟合环节硬编码频率：应复用已知频率或FFT提取的准确频率，避免固定值导致拟合曲线与原始信号错位。

修正后的完整代码 #

import numpy
import pylab as plt

def cosfunc(time, amplitude, omega, phase, offset):
    '''生成正弦波，相位单位为弧度'''
    return amplitude * numpy.cos(omega * time + phase) + offset

def get_cosine_approx(timeline, sine_data, known_freq):
    points_num = len(timeline)
    dt = timeline[1] - timeline[0]
    
    # 1. 计算直流偏移
    guess_offset = numpy.mean(sine_data)
    # 去除直流分量后的信号
    signal_no_dc = sine_data - guess_offset
    
    # 2. 计算FFT（针对去直流后的信号）
    fft_result = numpy.fft.fft(signal_no_dc)
    fft_freq = numpy.fft.fftfreq(points_num, dt)
    
    # 3. 定位已知频率对应的索引（考虑浮点精度，找最接近的频率点）
    target_freq_idx = numpy.argmin(numpy.abs(fft_freq - known_freq))
    
    # 4. 计算幅度：FFT幅值需除以采样点数再乘以2（正频率分量）
    guess_amplitude = numpy.abs(fft_result[target_freq_idx]) / points_num * 2
    
    # 5. 计算相位：FFT相位是exp(-jωt)的相位，需转换为cos(ωt+φ)的相位
    fft_phase = numpy.angle(fft_result[target_freq_idx])
    # 转换公式：cos(ωt+φ) = Re[exp(j(ωt+φ))] = Re[exp(jφ)exp(jωt)]，而FFT对应exp(-jωt)的系数，因此相位需取反
    guess_phase = -fft_phase
    # 可选：将相位归一化到[-π, π]区间
    guess_phase = numpy.mod(guess_phase + numpy.pi, 2*numpy.pi) - numpy.pi
    guess_phase_deg = guess_phase * 180. / numpy.pi
    
    # 输出结果
    print(f"已知频率: {known_freq:.2e} Hz")
    print(f"直流偏移: {guess_offset:.4f}")
    print(f"幅度: {guess_amplitude:.4f}")
    print(f"相位（弧度）: {guess_phase:.4f}")
    print(f"相位（角度）: {guess_phase_deg:.2f}°")
    
    # 绘图对比
    plt.figure(figsize=(10,6))
    plt.plot(timeline, sine_data, "ok", markersize=4, label="原始含噪信号")
    # 生成高密度时间轴用于拟合曲线
    new_timeline = numpy.linspace(timeline[0], timeline[-1], len(timeline)*100)
    fit_curve = cosfunc(new_timeline, guess_amplitude, 2*numpy.pi*known_freq, guess_phase, guess_offset)
    plt.plot(new_timeline, fit_curve, "r-", linewidth=1.5, label="拟合曲线")
    plt.xlabel("时间 (s)")
    plt.ylabel("幅值")
    plt.legend(loc="best")
    plt.grid(True, alpha=0.3)
    plt.show()
    
    return {
        "amp": guess_amplitude,
        "ph_rad": guess_phase,
        "ph_deg": guess_phase_deg,
        "offset": guess_offset
    }

# 测试参数
N = 256  # 采样点数
f = 56e9  # 已知频率（56GHz）
t = numpy.linspace(0.0, 100./f, N)  # 时间轴
omega = 2.*numpy.pi*f
offset = 0.5  # 加入直流偏移测试
phase = numpy.pi/4  # 加入相位测试（45度）
A = 1.8  # 幅度

# 生成含噪正弦波
clean_cos = cosfunc(t, A, omega, phase, offset)
noisy_cos = clean_cos + numpy.random.normal(0, 0.1, size=N)  # 添加高斯噪声

# 调用函数提取参数
result = get_cosine_approx(t, noisy_cos, f)

关键修改说明 #

1. 直流偏移计算：直接取信号的平均值，这是最直接的近似方法，也可通过FFT的0频分量（numpy.abs(fft_result[0])/points_num）得到一致结果。
2. 频率定位优化：因为已知目标频率，直接通过argmin找FFT频率轴上最接近的点，避免了FFT峰值偏移的问题（尤其在接近奈奎斯特频率时）。
3. 相位修正：FFT的基函数是exp(-jωt)，而我们的信号模型是cos(ωt+φ) = Re[exp(j(ωt+φ))]，两者的相位存在符号差异，因此需要将FFT相位取反才能匹配自定义模型的相位定义。
4. 去直流后再做FFT：避免直流分量对FFT低频段的干扰，提升幅度和相位计算的准确性。

内容的提问来源于stack exchange，提问作者Alex