Speex Wiener filter and rewriting of hypergeometric distribution

At present, noise reduction algorithms are being studied

There is a paragraph about Wiener filtering and hypergeometric distribution, which takes a lot of time to calculate , Look at the original code PDIV32_16,SHL32,EXTEND32 Wait, these macro definitions look annoying , Let's use the original format directly float Of C Rewrite the language directly

The text is as follows

        //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   Wiener filtering  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)));
            // Hypergeometric distribution gain parameter  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
             Hypergeometric distribution gain  MM=exp(-theta/2)*[(1+theta)*I0(theta/2)+theta*I1(theta/2)]; among I0 and I1 It's a Bessel function 
            */
            MM = hypergeom_gain(theta);
            /* EM gain with bound    gain  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 It is the gain after the gain calculated from the critical frequency is extended to the linear frequency domain 

            /* Constrain the gain to be close to the Bark scale gain  Limit Gain Values in Bark Domain gain 
              Constraint gain :   If  g/3>st->gain   be  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;  //   If  gain<gain_floor   be   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

        }

It was rewritten as follows

        //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 filtering  */
            prior_ratio = st->prior[i] / (st->prior[i] + 1.0f);
            // Hypergeometric distribution gain parameter 
            theta = prior_ratio * (1.0f + st->post[i]);
            /* Optimal estimator for loudness domain
             Hypergeometric distribution gain  MM=exp(-theta/2)*[(1+theta)*I0(theta/2)+theta*I1(theta/2)]; among I0 and I1 It's a Bessel function 
            */
            MM = hypergeom_gain(theta);
            /* EM gain with bound    gain  g=min(1,prior_ratio*mm)  */
            g = MIN32(Q15_ONE, prior_ratio * MM);
            /* Interpolated speech probability of presence  Calculate the weight or probability  */
            p = st->gain2[i]; //gain2 It is the gain after the gain calculated from the critical frequency is extended to the linear frequency domain 

            /* Constrain the gain to be close to the Bark scale gain  Limit Gain Values in Bark Domain gain 
              Constraint gain :   If  g/3>st->gain   be  g=3*st->gain
             */
            if ( (0.333f * g) > st->gain[i])
            {
                g = 3.0f * st->gain[i];
            }

            st->gain[i] = g;  //   If  gain<gain_floor   be   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 */
            /* Final amplitude spectral gain  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 Occupy 85%
            gpio_output(GPIO37, 0);
        }