C语言实现DWT信号去噪与PyWavelets效果不符求排查

我基于C语言实现了离散小波变换（DWT），目标是复现Python PyWavelets库的信号去噪功能，但在多阶小波变换时出现了严重问题：自研代码不仅无法有效去除低频噪声，还过度滤除了信号的峰值。（注：以下结果中蓝色为原始信号，橙色为去噪后信号）

• 1阶变换结果：原始峰值被轻微削弱，噪声未有效滤除
• 2阶变换结果：峰值进一步被压平，低频噪声依然存在
• 3阶变换结果：信号几乎失去原始峰值特征，整体过度平滑
• 4阶变换结果：信号严重失真，原始峰值基本消失
• PyWavelets 4阶变换结果：去噪信号完整保留原始峰值，同时噪声被有效滤除，干净度远优于自研实现

恳请熟悉小波变换的开发者帮忙排查以下代码问题：

void forward(float* arrTime, float* out, int size) {

    int h = size >> 1; // .. -> 8 -> 4 -> 2 .. shrinks in each step by half wavelength
    for( int i = 0; i < h; i++ ) {

        out[ i ] = out[ i + h ] = 0; // set to zero before sum up

        for( int j = 0; j < _motherWavelength; j++ ) {

            int k = ( i << 1 ) + j; // k = ( i * 2 ) + j;
            while( k >= size)
                k -= size; // circulate over arrays if scaling and wavelet are are larger

            out[ i ] += arrTime[ k ] * _scalingDeCom[ j ]; // low pass filter for the energy (approximation)
            out[ i + h ] += arrTime[ k ] * _waveletDeCom[ j ]; // high pass filter for the details

        } // Sorting each step in patterns of: { scaling coefficients | wavelet coefficients }

    } // h = 2^(p-1) | p = { 1, 2, .., N } .. shrinks in each step by half wavelength

} // forward

void threshold(float* coeffs, int size)
{
    // soft thresholding
    float sigma = 10e-4;
    float thres = sigma * sqrt(2 * log(size));
    for (int i = (size >> 1); i < size; i++)
    {
        if (abs(coeffs[i]) < thres)
        {
            coeffs[i] = 0;
        } else {
            coeffs[i] = (coeffs[i]/abs(coeffs[i])) * (abs(coeffs[i]) - thres);
        }


    }
}

void reverse(float* arrHilb, float* out, int size) {

    for( int i = 0; i < (sizeof(out)/sizeof(float)); i++ )
        out[ i ] = 0; // set to zero before sum up

    int h = size >> 1; // .. -> 8 -> 4 -> 2 .. shrinks in each step by half wavelength
    for( int i = 0; i < h; i++ ) {

        for( int j = 0; j < _motherWavelength; j++ ) {

            int k = ( i << 1 ) + j; // k = ( i * 2 ) + j;
            while( k >= size )
                k -= size; // circulate over arrays if scaling and wavelet are larger

            // adding up energy from low pass (approximation) and details from high pass filter
            out[ k ] +=
                    ( arrHilb[ i ] * _scalingReCon[ j ] )
                    + ( arrHilb[ i + h ] * _waveletReCon[ j ] ); // looks better with brackets

        } // Reconstruction from patterns of: { scaling coefficients | wavelet coefficients }

    } // h = 2^(p-1) | p = { 1, 2, .., N } .. shrink in each step by half wavelength

} // reverse

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