导数法检测PSD图峰值点的可靠性问题技术问询

问题背景 #

我正在开展一项信号处理任务，需在特定频率范围（Freqs[0]至Freqs[2]）内检测Power Spectral Density（PSD）图的峰值点。允许峰值处于目标频率±0.2Hz范围内。目前使用导数法（寻找斜率由正变负的位置）无法可靠捕获峰值，尤其是图中红、绿、蓝线标注区域；简单阈值法也无效，因为前段噪声与峰值幅度相同。

PSD图如下：

当前简化版代码：

def diff_f(f, a):
    return (f[a+1] - f[a]) / 0.2  # a是对应频率的索引

find_peak = []

# 检测Freqs[0]附近的峰值示例
if diff_f(f1, int(round((Freqs[0] / 0.2)) - 1)) > 0 and diff_f(f1, int(round((Freqs[0] / 0.2)))) < 0:
    find_peak.append(Freqs[0]) # 0.2是FFT的频率分辨率
# 对Freqs[1]、Freqs[2]及其他组合做类似检查

导数法存在的问题 #

• 单点斜率判断抗噪性差：仅依赖相邻两点的斜率变化，PSD图的局部噪声波动很容易导致误判或漏判，标注区域的峰值可能并非严格在目标频率索引处出现正负斜率反转。
• 检查范围过窄：代码只固定检查目标频率对应索引及其前一个位置的斜率，而峰值可能落在±0.2Hz范围内的其他索引上，直接导致漏检。
• 未验证局部极大性：斜率由正变负只能说明该点是局部下降点，但无法保证是周围区域的最大值，容易把噪声波动误识别为峰值。

改进方案 #

1. 扩大目标频率的检查窗口 #

既然允许峰值在±0.2Hz范围内（对应0.2Hz分辨率下前后各1个索引），可在窗口内寻找局部极大值，再验证斜率变化：

def find_local_peak_in_window(psd_data, target_freq, freq_res=0.2):
    center_idx = int(round(target_freq / freq_res))
    # 定义±0.2Hz的检查窗口，可按需调整范围
    window_start = max(0, center_idx - 1)
    window_end = min(len(psd_data)-1, center_idx + 1)
    
    # 找到窗口内最大值的索引
    window_values = psd_data[window_start:window_end+1]
    max_val_idx_in_window = window_values.index(max(window_values))
    peak_idx = window_start + max_val_idx_in_window
    
    # 验证斜率变化：左升右降
    left_slope = (psd_data[peak_idx] - psd_data[peak_idx-1])/freq_res if peak_idx > 0 else float('inf')
    right_slope = (psd_data[peak_idx+1] - psd_data[peak_idx])/freq_res if peak_idx < len(psd_data)-1 else -float('inf')
    
    if left_slope >= 0 and right_slope <= 0:
        return peak_idx * freq_res
    return None

# 批量检测目标频率的峰值
find_peak = []
for freq in Freqs[:3]:
    peak = find_local_peak_in_window(f1, freq)
    if peak is not None:
        find_peak.append(peak)

2. 先平滑PSD数据再检测 #

针对噪声干扰，先对PSD数据做平滑处理（比如移动平均），减少局部波动的影响：

import numpy as np

def smooth_psd(psd_data, window_size=3):
    # 移动平均平滑
    return np.convolve(psd_data, np.ones(window_size)/window_size, mode='same')

# 平滑后再执行峰值检测
smoothed_f1 = smooth_psd(f1)
find_peak = []
for freq in Freqs[:3]:
    peak = find_local_peak_in_window(smoothed_f1, freq)
    if peak is not None:
        find_peak.append(peak)

3. 结合邻域比较验证峰值 #

除了斜率变化，要求峰值点幅度大于邻域内所有点（比如前后2个点），确保是真正的局部极大值：
在find_local_peak_in_window函数中添加邻域验证逻辑：

# 验证该点是否为邻域内的最大值
is_local_max = True
# 检查前后各2个点的范围
for i in range(max(0, peak_idx-2), min(len(psd_data), peak_idx+3)):
    if i == peak_idx:
        continue
    if psd_data[i] > psd_data[peak_idx]:
        is_local_max = False
        break

if left_slope >=0 and right_slope <=0 and is_local_max:
    return peak_idx * freq_res

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