PPG信号运动伪影去除：小波阈值代码无效果问题求助

我研读文献后得知，可通过小波变换对系数施加阈值的方式去除PPG信号伪影。参照相关文献编写Python代码后，发现代码对输入信号完全没有作用。需求是对第1至3层小波系数施加阈值，且各层因幅值差异使用不同阈值，原代码如下：

coeff = pywt.wavedec(filterd2, wavelet='db4', level=3)
final_y=[]

for k in range(1, 4):
    temp = np.zeros(coeff[k].shape)
    for i in range(0,len(coeff[k])): 
        # print(f'gama={gamaathresh},alpha={alpha}')
        sigma = (1/0.6745) * madev(coeff[k])
        alphathresh = sigma * np.sqrt(2 * np.log(len(coeff[k])))
        gamaathresh=0.8*alphathresh
        alpha=-(1/gamaathresh)*(-0.95)
        if coeff[k][i]<=gamaathresh or coeff[k][i]>=(-gamaathresh) :
            # print(f'x<thresh= {x[i]}')
            #final_y.append(coeff[1])
            temp[i] = coeff[k][i]
        else:
            coeff[k][i]=(1-np.exp(-alpha*coeff[k][i]))/(1+np.exp(-alpha*coeff[k][i]))*alphathresh
            # print(f'x>thresh= {x[i]}')
            #final_y.append(coeff[1])
            temp[i] = coeff[k][i]
            # print(coeff[1].shape)
    final_y.append(temp)
print(size(final_y))
plt.plot(filterd2)
plt.show()
signal=pywt.waverec(final_y,wavelet='db4')
plt.plot(signal)

代码无效的核心问题： #

• 重构系数结构不完整：pywt.waverec要求传入包含近似系数（coeff[0]）和所有细节系数的完整列表，原代码只放入了第1-3层细节系数，漏掉了近似系数，导致重构逻辑完全错误。
• 阈值判断逻辑错误：条件coeff[k][i]<=gamaathresh or coeff[k][i]>=(-gamaathresh)覆盖了所有可能的数值（任何数要么小于等于正阈值，要么大于等于负阈值），else分支永远不会执行，等于没有做任何阈值处理。正确逻辑应为判断系数绝对值是否超过阈值。
• 系数处理无实际修改：无论if还是else分支，最终都将原系数赋值给temp[i]，没有对超过阈值的系数进行有效修改。
• 阈值参数重复计算：在每个系数的循环内重复计算sigma、阈值等参数，这些参数是针对整层系数的，应放在内层循环外，避免冗余计算。

修正后的代码示例： #

import pywt
import numpy as np
import matplotlib.pyplot as plt

def madev(d, axis=None):
    # 实现中位数绝对偏差计算，需确保与文献一致
    return np.median(np.abs(d - np.median(d, axis)), axis)

# 假设filterd2为输入PPG信号
coeff = pywt.wavedec(filterd2, wavelet='db4', level=3)
final_y = [coeff[0]]  # 保留近似系数，重构必须包含

for k in range(1, 4):
    # 针对当前层一次性计算阈值参数
    sigma = (1/0.6745) * madev(coeff[k])
    alphathresh = sigma * np.sqrt(2 * np.log(len(coeff[k])))
    gamaathresh = 0.8 * alphathresh
    alpha = 0.95 / gamaathresh  # 简化负负得正的计算
    
    temp = np.copy(coeff[k])
    for i in range(len(temp)):
        # 绝对值超过阈值时执行非线性阈值处理
        if abs(temp[i]) > gamaathresh:
            temp[i] = (1 - np.exp(-alpha * temp[i])) / (1 + np.exp(-alpha * temp[i])) * alphathresh
    final_y.append(temp)

# 重构信号
signal = pywt.waverec(final_y, wavelet='db4')

# 对比可视化
plt.figure(figsize=(12, 6))
plt.subplot(211)
plt.plot(filterd2, label='原始信号')
plt.legend()
plt.title('原始PPG信号')
plt.subplot(212)
plt.plot(signal, label='去伪影后信号')
plt.legend()
plt.title('小波阈值去伪影后信号')
plt.tight_layout()
plt.show()

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