ps_dec.c

tcpmp.src.0.72RC1 优秀的多媒体播放器TCPMP的源代码

                }

                /* accumulate energy */
#ifdef FIXED_POINT
                /* NOTE: all input is scaled by 2^(-5) because of fixed point QMF
                 * meaning that P will be scaled by 2^(-10) compared to floating point version
                 */
                in_re = ((abs(RE(inputLeft))+(1<<(REAL_BITS-1)))>>REAL_BITS);
                in_im = ((abs(IM(inputLeft))+(1<<(REAL_BITS-1)))>>REAL_BITS);
                P[n][bk] += in_re*in_re + in_im*in_im;
#else
                P[n][bk] += MUL_R(RE(inputLeft),RE(inputLeft)) + MUL_R(IM(inputLeft),IM(inputLeft));
#endif
            }
        }
    }

#if 0
    for (n = 0; n < 32; n++)
    {
        for (bk = 0; bk < 34; bk++)
        {
#ifdef FIXED_POINT
            printf("%d %d: %d\n", n, bk, P[n][bk] /*/(float)REAL_PRECISION*/);
#else
            printf("%d %d: %f\n", n, bk, P[n][bk]/1024.0);
#endif
        }
    }
#endif

    /* calculate transient reduction ratio for each parameter band b(k) */
    for (bk = 0; bk < ps->nr_par_bands; bk++)
    {
        for (n = ps->border_position[0]; n < ps->border_position[ps->num_env]; n++)
        {
            const real_t gamma = COEF_CONST(1.5);

            ps->P_PeakDecayNrg[bk] = MUL_F(ps->P_PeakDecayNrg[bk], ps->alpha_decay);
            if (ps->P_PeakDecayNrg[bk] < P[n][bk])
                ps->P_PeakDecayNrg[bk] = P[n][bk];

            /* apply smoothing filter to peak decay energy */
            P_SmoothPeakDecayDiffNrg = ps->P_SmoothPeakDecayDiffNrg_prev[bk];
            P_SmoothPeakDecayDiffNrg += MUL_F((ps->P_PeakDecayNrg[bk] - P[n][bk] - ps->P_SmoothPeakDecayDiffNrg_prev[bk]), ps->alpha_smooth);
            ps->P_SmoothPeakDecayDiffNrg_prev[bk] = P_SmoothPeakDecayDiffNrg;

            /* apply smoothing filter to energy */
            nrg = ps->P_prev[bk];
            nrg += MUL_F((P[n][bk] - ps->P_prev[bk]), ps->alpha_smooth);
            ps->P_prev[bk] = nrg;

            /* calculate transient ratio */
            if (MUL_C(P_SmoothPeakDecayDiffNrg, gamma) <= nrg)
            {
                G_TransientRatio[n][bk] = REAL_CONST(1.0);
            } else {
                G_TransientRatio[n][bk] = DIV_R(nrg, (MUL_C(P_SmoothPeakDecayDiffNrg, gamma)));
            }
        }
    }

#if 0
    for (n = 0; n < 32; n++)
    {
        for (bk = 0; bk < 34; bk++)
        {
#ifdef FIXED_POINT
            printf("%d %d: %f\n", n, bk, G_TransientRatio[n][bk]/(float)REAL_PRECISION);
#else
            printf("%d %d: %f\n", n, bk, G_TransientRatio[n][bk]);
#endif
        }
    }
#endif

    /* apply stereo decorrelation filter to the signal */
    for (gr = 0; gr < ps->num_groups; gr++)
    {
        if (gr < ps->num_hybrid_groups)
            maxsb = ps->group_border[gr] + 1;
        else
            maxsb = ps->group_border[gr + 1];

        /* QMF channel */
        for (sb = ps->group_border[gr]; sb < maxsb; sb++)
        {
            real_t g_DecaySlope;
            real_t g_DecaySlope_filt[NO_ALLPASS_LINKS];

            /* g_DecaySlope: [0..1] */
            if (gr < ps->num_hybrid_groups || sb <= ps->decay_cutoff)
            {
                g_DecaySlope = FRAC_CONST(1.0);
            } else {
                int8_t decay = ps->decay_cutoff - sb;
                if (decay <= -20 /* -1/DECAY_SLOPE */)
                {
                    g_DecaySlope = 0;
                } else {
                    /* decay(int)*decay_slope(frac) = g_DecaySlope(frac) */
                    g_DecaySlope = FRAC_CONST(1.0) + DECAY_SLOPE * decay;
                }
            }

            /* calculate g_DecaySlope_filt for every m multiplied by filter_a[m] */
            for (m = 0; m < NO_ALLPASS_LINKS; m++)
            {
                g_DecaySlope_filt[m] = MUL_F(g_DecaySlope, filter_a[m]);
            }


            /* set delay indices */
            temp_delay = ps->saved_delay;
            for (n = 0; n < NO_ALLPASS_LINKS; n++)
                temp_delay_ser[n] = ps->delay_buf_index_ser[n];

            for (n = ps->border_position[0]; n < ps->border_position[ps->num_env]; n++)
            {
                complex_t tmp, tmp0, R0;

                if (gr < ps->num_hybrid_groups)
                {
                    /* hybrid filterbank input */
                    RE(inputLeft) = QMF_RE(X_hybrid_left[n][sb]);
                    IM(inputLeft) = QMF_IM(X_hybrid_left[n][sb]);
                } else {
                    /* QMF filterbank input */
                    RE(inputLeft) = QMF_RE(X_left[n][sb]);
                    IM(inputLeft) = QMF_IM(X_left[n][sb]);
                }

                if (sb > ps->nr_allpass_bands && gr >= ps->num_hybrid_groups)
                {
                    /* delay */

                    /* never hybrid subbands here, always QMF subbands */
                    RE(tmp) = RE(ps->delay_Qmf[ps->delay_buf_index_delay[sb]][sb]);
                    IM(tmp) = IM(ps->delay_Qmf[ps->delay_buf_index_delay[sb]][sb]);
                    RE(R0) = RE(tmp);
                    IM(R0) = IM(tmp);
                    RE(ps->delay_Qmf[ps->delay_buf_index_delay[sb]][sb]) = RE(inputLeft);
                    IM(ps->delay_Qmf[ps->delay_buf_index_delay[sb]][sb]) = IM(inputLeft);
                } else {
                    /* allpass filter */
                    uint8_t m;
                    complex_t Phi_Fract;

                    /* fetch parameters */
                    if (gr < ps->num_hybrid_groups)
                    {
                        /* select data from the hybrid subbands */
                        RE(tmp0) = RE(ps->delay_SubQmf[temp_delay][sb]);
                        IM(tmp0) = IM(ps->delay_SubQmf[temp_delay][sb]);

                        RE(ps->delay_SubQmf[temp_delay][sb]) = RE(inputLeft);
                        IM(ps->delay_SubQmf[temp_delay][sb]) = IM(inputLeft);

                        RE(Phi_Fract) = RE(Phi_Fract_SubQmf[sb]);
                        IM(Phi_Fract) = IM(Phi_Fract_SubQmf[sb]);
                    } else {
                        /* select data from the QMF subbands */
                        RE(tmp0) = RE(ps->delay_Qmf[temp_delay][sb]);
                        IM(tmp0) = IM(ps->delay_Qmf[temp_delay][sb]);

                        RE(ps->delay_Qmf[temp_delay][sb]) = RE(inputLeft);
                        IM(ps->delay_Qmf[temp_delay][sb]) = IM(inputLeft);

                        RE(Phi_Fract) = RE(Phi_Fract_Qmf[sb]);
                        IM(Phi_Fract) = IM(Phi_Fract_Qmf[sb]);
                    }

                    /* z^(-2) * Phi_Fract[k] */
                    ComplexMult(&RE(tmp), &IM(tmp), RE(tmp0), IM(tmp0), RE(Phi_Fract), IM(Phi_Fract));

                    RE(R0) = RE(tmp);
                    IM(R0) = IM(tmp);
                    for (m = 0; m < NO_ALLPASS_LINKS; m++)
                    {
                        complex_t Q_Fract_allpass, tmp2;

                        /* fetch parameters */
                        if (gr < ps->num_hybrid_groups)
                        {
                            /* select data from the hybrid subbands */
                            RE(tmp0) = RE(ps->delay_SubQmf_ser[m][temp_delay_ser[m]][sb]);
                            IM(tmp0) = IM(ps->delay_SubQmf_ser[m][temp_delay_ser[m]][sb]);

                            if (ps->use34hybrid_bands)
                            {
                                RE(Q_Fract_allpass) = RE(Q_Fract_allpass_SubQmf34[sb][m]);
                                IM(Q_Fract_allpass) = IM(Q_Fract_allpass_SubQmf34[sb][m]);
                            } else {
                                RE(Q_Fract_allpass) = RE(Q_Fract_allpass_SubQmf20[sb][m]);
                                IM(Q_Fract_allpass) = IM(Q_Fract_allpass_SubQmf20[sb][m]);
                            }
                        } else {
                            /* select data from the QMF subbands */
                            RE(tmp0) = RE(ps->delay_Qmf_ser[m][temp_delay_ser[m]][sb]);
                            IM(tmp0) = IM(ps->delay_Qmf_ser[m][temp_delay_ser[m]][sb]);

                            RE(Q_Fract_allpass) = RE(Q_Fract_allpass_Qmf[sb][m]);
                            IM(Q_Fract_allpass) = IM(Q_Fract_allpass_Qmf[sb][m]);
                        }

                        /* delay by a fraction */
                        /* z^(-d(m)) * Q_Fract_allpass[k,m] */
                        ComplexMult(&RE(tmp), &IM(tmp), RE(tmp0), IM(tmp0), RE(Q_Fract_allpass), IM(Q_Fract_allpass));

                        /* -a(m) * g_DecaySlope[k] */
                        RE(tmp) += -MUL_F(g_DecaySlope_filt[m], RE(R0));
                        IM(tmp) += -MUL_F(g_DecaySlope_filt[m], IM(R0));

                        /* -a(m) * g_DecaySlope[k] * Q_Fract_allpass[k,m] * z^(-d(m)) */
                        RE(tmp2) = RE(R0) + MUL_F(g_DecaySlope_filt[m], RE(tmp));
                        IM(tmp2) = IM(R0) + MUL_F(g_DecaySlope_filt[m], IM(tmp));

                        /* store sample */
                        if (gr < ps->num_hybrid_groups)
                        {
                            RE(ps->delay_SubQmf_ser[m][temp_delay_ser[m]][sb]) = RE(tmp2);
                            IM(ps->delay_SubQmf_ser[m][temp_delay_ser[m]][sb]) = IM(tmp2);
                        } else {
                            RE(ps->delay_Qmf_ser[m][temp_delay_ser[m]][sb]) = RE(tmp2);
                            IM(ps->delay_Qmf_ser[m][temp_delay_ser[m]][sb]) = IM(tmp2);
                        }

                        /* store for next iteration (or as output value if last iteration) */
                        RE(R0) = RE(tmp);
                        IM(R0) = IM(tmp);
                    }
                }

                /* select b(k) for reading the transient ratio */
                bk = (~NEGATE_IPD_MASK) & ps->map_group2bk[gr];

                /* duck if a past transient is found */
                RE(R0) = MUL_R(G_TransientRatio[n][bk], RE(R0));
                IM(R0) = MUL_R(G_TransientRatio[n][bk], IM(R0));

                if (gr < ps->num_hybrid_groups)
                {
                    /* hybrid */
                    QMF_RE(X_hybrid_right[n][sb]) = RE(R0);
                    QMF_IM(X_hybrid_right[n][sb]) = IM(R0);
                } else {
                    /* QMF */
                    QMF_RE(X_right[n][sb]) = RE(R0);
                    QMF_IM(X_right[n][sb]) = IM(R0);
                }

                /* Update delay buffer index */
                if (++temp_delay >= 2)
                {
                    temp_delay = 0;
                }

                /* update delay indices */
                if (sb > ps->nr_allpass_bands && gr >= ps->num_hybrid_groups)
                {
                    /* delay_D depends on the samplerate, it can hold the values 14 and 1 */
                    if (++ps->delay_buf_index_delay[sb] >= ps->delay_D[sb])
                    {
                        ps->delay_buf_index_delay[sb] = 0;
                    }
                }

                for (m = 0; m < NO_ALLPASS_LINKS; m++)
                {
                    if (++temp_delay_ser[m] >= ps->num_sample_delay_ser[m])
                    {
                        temp_delay_ser[m] = 0;
                    }
                }
            }
        }
    }

    /* update delay indices */
    ps->saved_delay = temp_delay;
    for (m = 0; m < NO_ALLPASS_LINKS; m++)
        ps->delay_buf_index_ser[m] = temp_delay_ser[m];
}

#ifdef FIXED_POINT
#define step(shift) \
    if ((0x40000000l >> shift) + root <= value)       \
    {                                                 \
        value -= (0x40000000l >> shift) + root;       \
        root = (root >> 1) | (0x40000000l >> shift);  \
    } else {                                          \
        root = root >> 1;                             \
    }

/* fixed point square root approximation */
static real_t ps_sqrt(real_t value)
{
    real_t root = 0;

    step( 0); step( 2); step( 4); step( 6);
    step( 8); step(10); step(12); step(14);
    step(16); step(18); step(20); step(22);
    step(24); step(26); step(28); step(30);

    if (root < value)
        ++root;

    root <<= (REAL_BITS/2);

    return root;
}
#else
#define ps_sqrt(A) sqrt(A)
#endif

static const real_t ipdopd_cos_tab[] = {
    FRAC_CONST(1.000000000000000),
    FRAC_CONST(0.707106781186548),
    FRAC_CONST(0.000000000000000),
    FRAC_CONST(-0.707106781186547),
    FRAC_CONST(-1.000000000000000),
    FRAC_CONST(-0.707106781186548),
    FRAC_CONST(-0.000000000000000),
    FRAC_CONST(0.707106781186547),
    FRAC_CONST(1.000000000000000)
};

static const real_t ipdopd_sin_tab[] = {
    FRAC_CONST(0.000000000000000),
    FRAC_CONST(0.707106781186547),
    FRAC_CONST(1.000000000000000),
    FRAC_CONST(0.707106781186548),
    FRAC_CONST(0.000000000000000),
    FRAC_CONST(-0.707106781186547),
    FRAC_CONST(-1.000000000000000),
    FRAC_CONST(-0.707106781186548),
    FRAC_CONST(-0.000000000000000)
};

static void ps_mix_phase(ps_info *ps, qmf_t X_left[38][64], qmf_t X_right[38][64],
                         qmf_t X_hybrid_left[32][32], qmf_t X_hybrid_right[32][32])
{
    uint8_t n;
    uint8_t gr;
    uint8_t bk = 0;
    uint8_t sb, maxsb;
    uint8_t env;
    uint8_t nr_ipdopd_par;
    complex_t h11, h12, h21, h22;
    complex_t H11, H12, H21, H22;
    complex_t deltaH11, deltaH12, deltaH21, deltaH22;
    complex_t tempLeft;
    complex_t tempRight;
    complex_t phaseLeft;
    complex_t phaseRight;
    real_t L;
    const real_t *sf_iid;
    uint8_t no_iid_steps;

    if (ps->iid_mode >= 3)
    {