Speex维纳滤波与超几何分布的改写

目前在研究降噪算法

里面有一段关于维纳滤波和超几何分布计算耗时比较大,看原来的代码PDIV32_16,SHL32,EXTEND32等等这些宏定义看着很烦躁,下面直接把原来的格式用float的C语言方式直接改写一下

原文如下

        //N = 128; M = 24
        for (i = 0; i < N; i++)
        {
            spx_word32_t MM;
            spx_word32_t theta;
            spx_word16_t prior_ratio;
            spx_word16_t tmp;
            spx_word16_t p;
            spx_word16_t g;

            /* Wiener filter gain  维纳滤波 prior_snr=   prior_snr /(prior_snr +1)   */
            prior_ratio = PDIV32_16(SHL32(EXTEND32(st->prior[i]), 15), ADD16(st->prior[i], SHL32(1, SNR_SHIFT)));
            //超几何分布增益参数 theta= prior_snr *(1+post_snr)；
            theta = MULT16_32_P15(prior_ratio, QCONST32(1.f, EXPIN_SHIFT) + SHL32(EXTEND32(st->post[i]), EXPIN_SHIFT - SNR_SHIFT));

            /* Optimal estimator for loudness domain
            超几何分布增益 MM=exp(-theta/2)*[(1+theta)*I0(theta/2)+theta*I1(theta/2)];其中I0和I1是贝塞尔函数
            */
            MM = hypergeom_gain(theta);
            /* EM gain with bound   增益 g=min(1,prior_ratio*mm)  */
            g = EXTRACT16(MIN32(Q15_ONE, MULT16_32_Q15(prior_ratio, MM)));
            /* Interpolated speech probability of presence */
            p = st->gain2[i]; //gain2是由临界频率计算后的增益扩展到线性频域后的增益

            /* Constrain the gain to be close to the Bark scale gain 限制Gain值在Bark域增益
             约束增益:  如果 g/3>st->gain  则 g=3*st->gain
             */
            if (MULT16_16_Q15(QCONST16(.333f, 15), g) > st->gain[i])
            {
                g = MULT16_16(3, st->gain[i]);
            }

            st->gain[i] = g;  //  如果 gain<gain_floor  则  gain =gain_floor

            /* Save old power spectrum */
            st->old_ps[i] = MULT16_32_P15(QCONST16(.2f, 15), st->old_ps[i]) + MULT16_32_P15(MULT16_16_P15(QCONST16(.8f, 15), SQR16_Q15(st->gain[i])), ps[i]);

            /* Apply gain floor */
            if (st->gain[i] < st->gain_floor[i])
            {
                st->gain[i] = st->gain_floor[i];
            }

            st->gain2[i] = powf(st->gain[i], p) * powf(st->gain_floor[i], 1.f - p); //=0.112ms

        }

改写后如下

        //N = 128; M = 24
        for (i = 0; i < N; i++)
        {
            spx_word32_t MM;
            spx_word32_t theta;
            spx_word16_t prior_ratio;
            spx_word16_t tmp;
            spx_word16_t p;
            spx_word16_t g;

            /* 维纳滤波 */
            prior_ratio = st->prior[i] / (st->prior[i] + 1.0f);
            //超几何分布增益参数
            theta = prior_ratio * (1.0f + st->post[i]);
            /* Optimal estimator for loudness domain
            超几何分布增益 MM=exp(-theta/2)*[(1+theta)*I0(theta/2)+theta*I1(theta/2)];其中I0和I1是贝塞尔函数
            */
            MM = hypergeom_gain(theta);
            /* EM gain with bound   增益 g=min(1,prior_ratio*mm)  */
            g = MIN32(Q15_ONE, prior_ratio * MM);
            /* Interpolated speech probability of presence 计算权重或者叫概率 */
            p = st->gain2[i]; //gain2是由临界频率计算后的增益扩展到线性频域后的增益

            /* Constrain the gain to be close to the Bark scale gain 限制Gain值在Bark域增益
             约束增益:  如果 g/3>st->gain  则 g=3*st->gain
             */
            if ( (0.333f * g) > st->gain[i])
            {
                g = 3.0f * st->gain[i];
            }

            st->gain[i] = g;  //  如果 gain<gain_floor  则  gain =gain_floor

            /* Save old power spectrum */
            st->old_ps[i] = 0.2f * st->old_ps[i] + 0.8f * st->gain[i] * st->gain[i] * ps[i];

            /* Apply gain floor */
            if (st->gain[i] < st->gain_floor[i])
            {
                st->gain[i] = st->gain_floor[i];
            }

            /* Use this if you want a log-domain MMSE estimator instead */
            /*最终幅度谱增益 gain2={p*sqrt(g)+(1-p)*sqrt(st->gain_floor  )}^2 */
            gpio_output(GPIO37, 1);
            st->gain2[i] = powf(st->gain[i], p) * powf(st->gain_floor[i], 1.f - p); //=0.112ms占用85%
            gpio_output(GPIO37, 0);
        }