与50kHz tone burst卷积信号功率谱峰值在零频的问题求解

问题背景 #

目前在计算与中心频率50kHz的tone burst信号卷积后的信号功率谱，相关参数如下：

• 采样率：30MHz
• 采样点数：30000
• 关注频段：0~500kHz

卷积后的时域信号如下图：

已对信号做去均值、去趋势处理，但计算得到的功率谱峰值出现在零频附近，没有出现预期的50kHz峰值，原始功率谱结果如下图：

原始计算代码如下：

fsig = 50e3
N = 3e4
for i in range(1):
    sampling_rate = 30e6
    data = df0.iloc[:, i] - np.mean(df0.iloc[:, i])
    fourier_transform = np.fft.rfft(data)
    abs_fourier_transform = np.abs(fourier_transform)
    power_spectrum = np.square(abs_fourier_transform)
    frequency = np.linspace(0, sampling_rate / 2, len(power_spectrum))
    frequency /= 1e3
    power_spectrum = (power_spectrum - np.min(power_spectrum)) / (np.max(power_spectrum) - np.min(power_spectrum))
    fig, ax = plt.subplots()
    ax.plot(frequency, power_spectrum)
    ax.set_xlim([0, 500])
    plt.xlabel('Frequency [kHz]')
    plt.ylabel('Normalised amplitude')
    plt.show()

问题原因分析 #

功率谱峰值集中在零频、看不到50kHz峰值的问题，主要由以下几个常见原因导致：

1. 归一化方式与坐标选择不合理：你使用了线性坐标+min-max归一化的方式展示功率谱，如果信号本身存在幅度稍高的低频分量（比如tone burst的包络频谱、系统低频噪声、少量残留直流），归一化后50kHz的峰值会被直接压缩到接近0的基线位置，完全无法分辨。
2. 频谱泄漏未抑制：tone burst是有限时长的脉冲信号，直接做FFT相当于加了矩形窗，会产生很高的旁瓣，50kHz的载波能量会泄漏到周围频点，导致主峰幅度被拉低，更容易被低频分量掩盖。
3. 频率轴生成方式存在隐患：你用np.linspace生成频率轴，虽然在当前偶数采样点数下数值偏差不大，但如果后续采样点数变化很容易出现频率轴偏移，更稳妥的方式是用FFT库自带的频率轴生成函数。

解决方法 #

按照以下步骤修改代码，即可观测到50kHz处的峰值：

1. 给输入信号加汉宁窗抑制频谱泄漏
2. 放弃min-max归一化，改用对数坐标（dB）展示功率谱，提升动态范围，让小幅度的峰也能被观测到
3. 用np.fft.rfftfreq生成准确的频率轴

修正后代码 #

import numpy as np
import matplotlib.pyplot as plt

fsig = 50e3
N = int(3e4)
sampling_rate = 30e6

for i in range(1):
    # 读取并去均值
    data = df0.iloc[:, i] - np.mean(df0.iloc[:, i])
    # 加汉宁窗抑制频谱泄漏
    data = data * np.hanning(len(data))
    # 做FFT
    fourier_transform = np.fft.rfft(data)
    abs_fourier_transform = np.abs(fourier_transform)
    power_spectrum = np.square(abs_fourier_transform)
    # 生成准确的频率轴
    frequency = np.fft.rfftfreq(n=len(data), d=1/sampling_rate)
    frequency /= 1e3
    # 转dB归一化，动态范围更高
    power_spectrum_dB = 10 * np.log10(power_spectrum / np.max(power_spectrum))
    
    # 绘图
    fig, ax = plt.subplots()
    ax.plot(frequency, power_spectrum_dB)
    ax.set_xlim([0, 500])
    plt.xlabel('Frequency [kHz]')
    plt.ylabel('Normalised amplitude [dB]')
    plt.show()

如果修改后还是看不到50kHz的峰，建议先检查卷积后的时域波形，确认确实存在50kHz的周期振荡，排除卷积参数配置错误、tone burst中心频率偏差、采样率设置错误等前置问题。

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