static void MatrixDecode(const float *in, const int k, const int il,
    const int ir, bool decode_rear,
    const int _dlbuflen,
    float _l_fwr, float _r_fwr,
    float _lpr_fwr, float _lmr_fwr,
    float *_adapt_l_gain, float *_adapt_r_gain,
    float *_adapt_lpr_gain, float *_adapt_lmr_gain,
    float *_lf, float *_rf, float *_lr,
    float *_rr, float *_cf)
{
    static const float M9_03DB = 0.3535533906f;
    static const float MATAGCTRIG = 8.0f;   /* (Fuzzy) AGC trigger */
    static const float MATAGCDECAY = 1.0f;  /* AGC baseline decay rate (1/samp.) */
    static const float MATCOMPGAIN = 0.37f; /* Cross talk compensation gain,  0.50 - 0.55 is full cancellation. */

    const int kr = (k + olddelay) % _dlbuflen;
    float l_gain = (_l_fwr + _r_fwr) / (1 + _l_fwr + _l_fwr);
    float r_gain = (_l_fwr + _r_fwr) / (1 + _r_fwr + _r_fwr);
    // The 2nd axis has strong gain fluctuations, and therefore require
    // limits.  The factor corresponds to the 1 / amplification of (Lt
    // - Rt) when (Lt, Rt) is strongly correlated. (e.g. during
    // dialogues).  It should be bigger than -12 dB to prevent
    // distortion.
    float lmr_lim_fwr = _lmr_fwr > M9_03DB * _lpr_fwr ? _lmr_fwr : M9_03DB * _lpr_fwr;
    float lpr_gain = (_lpr_fwr + lmr_lim_fwr) / (1 + _lpr_fwr + _lpr_fwr);
    float lmr_gain = (_lpr_fwr + lmr_lim_fwr) / (1 + lmr_lim_fwr + lmr_lim_fwr);
    float lmr_unlim_gain = (_lpr_fwr + _lmr_fwr) / (1 + _lmr_fwr + _lmr_fwr);
    float lpr, lmr;
    float l_agc, r_agc, lpr_agc, lmr_agc;
    float f, d_gain, c_gain, c_agc_cfk;

    /*** AXIS NO. 1: (Lt, Rt) -> (C, Ls, Rs) ***/
    /* AGC adaption */
    d_gain = (fabs(l_gain - *_adapt_l_gain) + fabs(r_gain - *_adapt_r_gain)) * 0.5f;
    f = d_gain * (1.0f / MATAGCTRIG);
    f = MATAGCDECAY - MATAGCDECAY / (1 + f * f);
    *_adapt_l_gain = (1 - f) * *_adapt_l_gain + f * l_gain;
    *_adapt_r_gain = (1 - f) * *_adapt_r_gain + f * r_gain;
    /* Matrix */
    l_agc = in[il] * PassiveLock(*_adapt_l_gain);
    r_agc = in[ir] * PassiveLock(*_adapt_r_gain);
    _cf[k] = (l_agc + r_agc) * (float)M_SQRT1_2;
    if (decode_rear)
    {
        _lr[kr] = _rr[kr] = (l_agc - r_agc) * (float)M_SQRT1_2;
        // Stereo rear channel is steered with the same AGC steering as
        // the decoding matrix. Note this requires a fast updating AGC
        // at the order of 20 ms (which is the case here).
        _lr[kr] *= (_l_fwr + _l_fwr) / (1 + _l_fwr + _r_fwr);
        _rr[kr] *= (_r_fwr + _r_fwr) / (1 + _l_fwr + _r_fwr);
    }

    /*** AXIS NO. 2: (Lt + Rt, Lt - Rt) -> (L, R) ***/
    lpr = (in[il] + in[ir]) * (float)M_SQRT1_2;
    lmr = (in[il] - in[ir]) * (float)M_SQRT1_2;
    /* AGC adaption */
    d_gain = fabs(lmr_unlim_gain - *_adapt_lmr_gain);
    f = d_gain * (1.0f / MATAGCTRIG);
    f = MATAGCDECAY - MATAGCDECAY / (1 + f * f);
    *_adapt_lpr_gain = (1 - f) * *_adapt_lpr_gain + f * lpr_gain;
    *_adapt_lmr_gain = (1 - f) * *_adapt_lmr_gain + f * lmr_gain;
    /* Matrix */
    lpr_agc = lpr * PassiveLock(*_adapt_lpr_gain);
    lmr_agc = lmr * PassiveLock(*_adapt_lmr_gain);
    _lf[k] = (lpr_agc + lmr_agc) * (float)M_SQRT1_2;
    _rf[k] = (lpr_agc - lmr_agc) * (float)M_SQRT1_2;

    /*** CENTER FRONT CANCELLATION ***/
    // A heuristic approach exploits that Lt + Rt gain contains the
    // information about Lt, Rt correlation.  This effectively reshapes
    // the front and rear "cones" to concentrate Lt + Rt to C and
    // introduce Lt - Rt in L, R.
    /* 0.67677 is the empirical lower bound for lpr_gain. */
    c_gain = 8 * (*_adapt_lpr_gain - 0.67677f);
    c_gain = c_gain > 0 ? c_gain : 0;
    // c_gain should not be too high, not even reaching full
    // cancellation (~ 0.50 - 0.55 at current AGC implementation), or
    // the center will sound too narrow. */
    c_gain = MATCOMPGAIN / (1 + c_gain * c_gain);
    c_agc_cfk = c_gain * _cf[k];
    _lf[k] -= c_agc_cfk;
    _rf[k] -= c_agc_cfk;
    _cf[k] += c_agc_cfk + c_agc_cfk;
}

void DPL2Decode(float *samples, int numsamples, float *out)
{
    static const unsigned int FWRDURATION = 240; // FWR average duration (samples)
    static const int cfg_delay = 0;
    static const unsigned int fmt_freq = 48000;
    static const unsigned int fmt_nchannels = 2; // input channels

    int cur = 0;

    if (olddelay != cfg_delay || oldfreq != fmt_freq)
    {
        Done();
        olddelay = cfg_delay;
        oldfreq = fmt_freq;
        dlbuflen = std::max(FWRDURATION, (fmt_freq * cfg_delay / 1000)); //+(len7000-1);
        cyc_pos = dlbuflen - 1;
        fwrbuf_l.resize(dlbuflen);
        fwrbuf_r.resize(dlbuflen);
        lf.resize(dlbuflen);
        rf.resize(dlbuflen);
        lr.resize(dlbuflen);
        rr.resize(dlbuflen);
        cf.resize(dlbuflen);
        cr.resize(dlbuflen);
        filter_coefs_lfe = CalculateCoefficients125HzLowpass(fmt_freq);
        lfe_pos = 0;
        memset(LFE_buf, 0, sizeof(LFE_buf));
    }

    float *in = samples; // Input audio data
    float *end = in + numsamples * fmt_nchannels; // Loop end

    while (in < end)
    {
        const int k = cyc_pos;

        const int fwr_pos = (k + FWRDURATION) % dlbuflen;
        /* Update the full wave rectified total amplitude */
        /* Input matrix decoder */
        l_fwr += fabs(in[0]) - fabs(fwrbuf_l[fwr_pos]);
        r_fwr += fabs(in[1]) - fabs(fwrbuf_r[fwr_pos]);
        lpr_fwr += fabs(in[0] + in[1]) - fabs(fwrbuf_l[fwr_pos] + fwrbuf_r[fwr_pos]);
        lmr_fwr += fabs(in[0] - in[1]) - fabs(fwrbuf_l[fwr_pos] - fwrbuf_r[fwr_pos]);

        /* Matrix encoded 2 channel sources */
        fwrbuf_l[k] = in[0];
        fwrbuf_r[k] = in[1];
        MatrixDecode(in, k, 0, 1, true, dlbuflen,
            l_fwr, r_fwr,
            lpr_fwr, lmr_fwr,
            &adapt_l_gain, &adapt_r_gain,
            &adapt_lpr_gain, &adapt_lmr_gain,
            &lf[0], &rf[0], &lr[0], &rr[0], &cf[0]);

        out[cur + 0] = lf[k];
        out[cur + 1] = rf[k];
        out[cur + 2] = cf[k];
        LFE_buf[lfe_pos] = (lf[k] + rf[k] + 2.0f * cf[k] + lr[k] + rr[k]) / 2.0f;
        out[cur + 3] = FIRFilter(LFE_buf, lfe_pos, len125, len125, filter_coefs_lfe);
        lfe_pos++;
        if (lfe_pos == len125)
        {
            lfe_pos = 0;
        }
        out[cur + 4] = lr[k];
        out[cur + 5] = rr[k];
        // Next sample...
        in += 2;
        cur += 6;
        cyc_pos--;
        if (cyc_pos < 0)
        {
            cyc_pos += dlbuflen;
        }
    }
}