Python中rPPG信号SNR计算逻辑与掩码作用咨询

以下是涉及的实现代码：

def pow2db(x):
    """returns the corresponding decibel (dB) value for a power value x.

    The relationship between power and decibels is:

    .. math::    X_{dB} = 10 * \log_{10}(x)

    .. doctest::

        >>> from spectrum import pow2db
        >>> x = pow2db(0.1)
        >>> x
        -10.0
    """
    return 10 * log10(x)

# calculate SNR from Pulse - Signal
def SNR(bvp, fps, hr):

    #BVPSNR Estimates the signal-to-noise ratio of the blood volume pulse
    # signal. Adapted from the method by G. de Haan, TBME, 2013.
    # SNR calculated as the ratio (in dB) of power contained within +/- 0.1 Hz
    # of the reference heart rate frequency and +/- 0.2 of its first
    # harmonic and sum of all other power between 0.7 and 2 Hz.
    # Adapted from the method by G. de Haan, TBME, 2013
    #
    #   Inputs:
    #       BVP                     = A BVP timeseries.
    #       FS                      = The sample rate of the BVP time series (Hz/fps).
    #       HR                      = The reference heart rate (Hz/fps).
    #       PlotTF                  = Boolean to turn plotting results on or off.
    #
    #  Outputs:
    #       SNR                     = Blood Volume Pulse Signal-to-Noise Ratio.
    #

    HR_F = hr/60
    Nyquist = fps/2
    FResBPM = 0.5
    N = (60*2*Nyquist)/FResBPM
    print(N)

    f_set, Pxx_den = welch(bvp, fps, window='hann', nperseg=N)
    GTMask1 = (f_set >= HR_F-0.15) & (f_set <= HR_F+0.15)
    GTMask2 = (f_set >= HR_F-0.3) & (f_set <= HR_F+0.3)

    SPower = sum(Pxx_den[GTMask1|GTMask2])
    FMask2 = (f_set >= 0.7)&(f_set <= 2)
    AllPower = sum(Pxx_den[FMask2])

    figure100 = plt.figure(100)
    plt.semilogy(f_set,Pxx_den)
    plt.title('welch Power Spectrum for SNR calculation / ROIs')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Power (dB)')

    return pow2db(SPower/(AllPower-SPower))

代码逐段逻辑拆解 #

1. 工具函数与参数预处理 #

pow2db是标准的功率值转分贝工具函数，遵循功率转分贝的通用公式X_dB = 10 * log10(x)，无特殊定制逻辑。
进入SNR函数后的参数处理规则：

• HR_F = hr/60：把输入单位为次/分钟的参考心率，转换为单位为Hz的基频值，比如66次/分的心率对应1.1Hz基频，也就是你观察到的1.1Hz目标点
• Welch功率谱参数计算：设置功率谱的频率分辨率为0.5BPM，据此计算Welch方法的分段长度N，保证输出的功率谱频率精度满足计算要求
• 调用welch方法加汉宁窗计算输入BVP信号的功率谱密度，得到频率轴序列f_set和对应频率点的功率值序列Pxx_den

2. 频段掩码设计逻辑 #

需要特别注意：这份代码的实际实现和顶部注释、原始论文方法存在不一致，不要完全按照注释理解逻辑：

• 原始de Haan 2013论文的规则是：信号功率包含两部分，一是心率基频f0为中心±0.1Hz范围内的功率，二是一次谐波2*f0为中心±0.2Hz范围内的功率
• 你贴的代码里两个掩码的中心都设置为了基频HR_F：
  • GTMask1：筛选基频左右0.15Hz范围内的频率点
  • GTMask2：筛选基频左右0.3Hz范围内的频率点
    两个范围是包含关系（±0.15Hz完全落在±0.3Hz区间内），做或运算后实际取的是基频为中心±0.3Hz的整个区间，没有覆盖谐波频段，属于代码实现的笔误。
• FMask2：筛选0.7Hz_{2Hz的频率范围，对应42}120次/分的生理合理心率区间，所有功率计算都限定在这个范围内，排除低频漂移、高频噪声等非生理相关信号的干扰。

3. 功率计算与SNR换算 #

• SPower：所有被判定为有效脉搏信号的频段功率总和，在这份代码的实际实现中，就是基频周围±0.3Hz区间的功率和
• AllPower：0.7~2Hz整个生理频段内的总功率
• 噪声功率为AllPower - SPower，即生理频段内扣除信号带之后的剩余功率
• 最终SNR计算：先求信号功率和噪声功率的比值，再调用pow2db转换为dB单位输出。

对你原有理解的修正 #

• 代码不是分别计算±0.15Hz、±0.3Hz两个区间的SNR，两个掩码是共同圈定信号频段的范围
• 不存在“判定目标频段无同等强度其他信号”的逻辑，计算规则是硬划分：指定范围内的功率全部算信号，范围外生理带内的功率全部算噪声，直接做比值
• 你观察到的1.1Hz是参考心率对应的基频位置，所有信号带都是围绕这个中心划定的。

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