libaf/af_hrtf.c

   1 /*
   2  * Experimental audio filter that mixes 5.1 and 5.1 with matrix
   3  * encoded rear channels into headphone signal using FIR filtering
   4  * with HRTF.
   5  *
   6  * This file is part of MPlayer.
   7  *
   8  * MPlayer is free software; you can redistribute it and/or modify
   9  * it under the terms of the GNU General Public License as published by
  10  * the Free Software Foundation; either version 2 of the License, or
  11  * (at your option) any later version.
  12  *
  13  * MPlayer is distributed in the hope that it will be useful,
  14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  * GNU General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU General Public License along
  19  * with MPlayer; if not, write to the Free Software Foundation, Inc.,
  20  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
  21  */
  22
  23 //#include <stdio.h>
  24 #include <stdlib.h>
  25 #include <string.h>
  26 #include <inttypes.h>
  27
  28 #include <math.h>
  29
  30 #include "af.h"
  31 #include "dsp.h"
  32
  33 /* HRTF filter coefficients and adjustable parameters */
  34 #include "af_hrtf.h"
  35
  36 typedef struct af_hrtf_s {
  37     /* Lengths */
  38     int dlbuflen, hrflen, basslen;
  39     /* L, C, R, Ls, Rs channels */
  40     float *lf, *rf, *lr, *rr, *cf, *cr;
  41     const float *cf_ir, *af_ir, *of_ir, *ar_ir, *or_ir, *cr_ir;
  42     int cf_o, af_o, of_o, ar_o, or_o, cr_o;
  43     /* Bass */
  44     float *ba_l, *ba_r;
  45     float *ba_ir;
  46     /* Whether to matrix decode the rear center channel */
  47     int matrix_mode;
  48     /* How to decode the input:
  49        0 = 5/5+1 channels
  50        1 = 2 channels
  51        2 = matrix encoded 2 channels */
  52     int decode_mode;
  53     /* Full wave rectified (FWR) amplitudes and gain used to steer the
  54        active matrix decoding of front channels (variable names
  55        lpr/lmr means Lt + Rt, Lt - Rt) */
  56     float l_fwr, r_fwr, lpr_fwr, lmr_fwr;
  57     float adapt_l_gain, adapt_r_gain, adapt_lpr_gain, adapt_lmr_gain;
  58     /* Matrix input decoding require special FWR buffer, since the
  59        decoding is done in place. */
  60     float *fwrbuf_l, *fwrbuf_r, *fwrbuf_lr, *fwrbuf_rr;
  61     /* Rear channel delay buffer for matrix decoding */
  62     float *rear_dlbuf;
  63     /* Full wave rectified amplitude and gain used to steer the active
  64        matrix decoding of center rear channel */
  65     float lr_fwr, rr_fwr, lrprr_fwr, lrmrr_fwr;
  66     float adapt_lr_gain, adapt_rr_gain;
  67     float adapt_lrprr_gain, adapt_lrmrr_gain;
  68     /* Cyclic position on the ring buffer */
  69     int cyc_pos;
  70     int print_flag;
  71 } af_hrtf_t;
  72
  73 /* Convolution on a ring buffer
  74  *    nx:       length of the ring buffer
  75  *    nk:       length of the convolution kernel
  76  *    sx:       ring buffer
  77  *    sk:       convolution kernel
  78  *    offset:   offset on the ring buffer, can be
  79  */
  80 static float conv(const int nx, const int nk, const float *sx, const float *sk,
  81                   const int offset)
  82 {
  83     /* k = reminder of offset / nx */
  84     int k = offset >= 0 ? offset % nx : nx + (offset % nx);
  85
  86     if(nk + k <= nx)
  87         return af_filter_fir(nk, sx + k, sk);
  88     else
  89         return af_filter_fir(nk + k - nx, sx, sk + nx - k) +
  90             af_filter_fir(nx - k, sx + k, sk);
  91 }
  92
  93 /* Detect when the impulse response starts (significantly) */
  94 static int pulse_detect(const float *sx)
  95 {
  96     /* nmax must be the reference impulse response length (128) minus
  97        s->hrflen */
  98     const int nmax = 128 - HRTFFILTLEN;
  99     const float thresh = IRTHRESH;
 100     int i;
 101
 102     for(i = 0; i < nmax; i++)
 103         if(fabs(sx[i]) > thresh)
 104             return i;
 105     return 0;
 106 }
 107
 108 /* Fuzzy matrix coefficient transfer function to "lock" the matrix on
 109    a effectively passive mode if the gain is approximately 1 */
 110 static inline float passive_lock(float x)
 111 {
 112    const float x1 = x - 1;
 113    const float ax1s = fabs(x - 1) * (1.0 / MATAGCLOCK);
 114
 115    return x1 - x1 / (1 + ax1s * ax1s) + 1;
 116 }
 117
 118 /* Unified active matrix decoder for 2 channel matrix encoded surround
 119    sources */
 120 static inline void matrix_decode(short *in, const int k, const int il,
 121                           const int ir, const int decode_rear,
 122                           const int dlbuflen,
 123                           float l_fwr, float r_fwr,
 124                           float lpr_fwr, float lmr_fwr,
 125                           float *adapt_l_gain, float *adapt_r_gain,
 126                           float *adapt_lpr_gain, float *adapt_lmr_gain,
 127                           float *lf, float *rf, float *lr,
 128                           float *rr, float *cf)
 129 {
 130    const int kr = (k + MATREARDELAY) % dlbuflen;
 131    float l_gain = (l_fwr + r_fwr) /
 132       (1 + l_fwr + l_fwr);
 133    float r_gain = (l_fwr + r_fwr) /
 134       (1 + r_fwr + r_fwr);
 135    /* The 2nd axis has strong gain fluctuations, and therefore require
 136       limits.  The factor corresponds to the 1 / amplification of (Lt
 137       - Rt) when (Lt, Rt) is strongly correlated. (e.g. during
 138       dialogues).  It should be bigger than -12 dB to prevent
 139       distortion. */
 140    float lmr_lim_fwr = lmr_fwr > M9_03DB * lpr_fwr ?
 141       lmr_fwr : M9_03DB * lpr_fwr;
 142    float lpr_gain = (lpr_fwr + lmr_lim_fwr) /
 143       (1 + lpr_fwr + lpr_fwr);
 144    float lmr_gain = (lpr_fwr + lmr_lim_fwr) /
 145       (1 + lmr_lim_fwr + lmr_lim_fwr);
 146    float lmr_unlim_gain = (lpr_fwr + lmr_fwr) /
 147       (1 + lmr_fwr + lmr_fwr);
 148    float lpr, lmr;
 149    float l_agc, r_agc, lpr_agc, lmr_agc;
 150    float f, d_gain, c_gain, c_agc_cfk;
 151
 152 #if 0
 153    static int counter = 0;
 154    static FILE *fp_out;
 155
 156    if(counter == 0)
 157       fp_out = fopen("af_hrtf.log", "w");
 158    if(counter % 240 == 0)
 159       fprintf(fp_out, "%g %g %g %g %g ", counter * (1.0 / 48000),
 160               l_gain, r_gain, lpr_gain, lmr_gain);
 161 #endif
 162
 163    /*** AXIS NO. 1: (Lt, Rt) -> (C, Ls, Rs) ***/
 164    /* AGC adaption */
 165    d_gain = (fabs(l_gain - *adapt_l_gain) +
 166              fabs(r_gain - *adapt_r_gain)) * 0.5;
 167    f = d_gain * (1.0 / MATAGCTRIG);
 168    f = MATAGCDECAY - MATAGCDECAY / (1 + f * f);
 169    *adapt_l_gain = (1 - f) * *adapt_l_gain + f * l_gain;
 170    *adapt_r_gain = (1 - f) * *adapt_r_gain + f * r_gain;
 171    /* Matrix */
 172    l_agc = in[il] * passive_lock(*adapt_l_gain);
 173    r_agc = in[ir] * passive_lock(*adapt_r_gain);
 174    cf[k] = (l_agc + r_agc) * M_SQRT1_2;
 175    if(decode_rear) {
 176       lr[kr] = rr[kr] = (l_agc - r_agc) * M_SQRT1_2;
 177       /* Stereo rear channel is steered with the same AGC steering as
 178          the decoding matrix. Note this requires a fast updating AGC
 179          at the order of 20 ms (which is the case here). */
 180       lr[kr] *= (l_fwr + l_fwr) /
 181          (1 + l_fwr + r_fwr);
 182       rr[kr] *= (r_fwr + r_fwr) /
 183          (1 + l_fwr + r_fwr);
 184    }
 185
 186    /*** AXIS NO. 2: (Lt + Rt, Lt - Rt) -> (L, R) ***/
 187    lpr = (in[il] + in[ir]) * M_SQRT1_2;
 188    lmr = (in[il] - in[ir]) * M_SQRT1_2;
 189    /* AGC adaption */
 190    d_gain = fabs(lmr_unlim_gain - *adapt_lmr_gain);
 191    f = d_gain * (1.0 / MATAGCTRIG);
 192    f = MATAGCDECAY - MATAGCDECAY / (1 + f * f);
 193    *adapt_lpr_gain = (1 - f) * *adapt_lpr_gain + f * lpr_gain;
 194    *adapt_lmr_gain = (1 - f) * *adapt_lmr_gain + f * lmr_gain;
 195    /* Matrix */
 196    lpr_agc = lpr * passive_lock(*adapt_lpr_gain);
 197    lmr_agc = lmr * passive_lock(*adapt_lmr_gain);
 198    lf[k] = (lpr_agc + lmr_agc) * M_SQRT1_2;
 199    rf[k] = (lpr_agc - lmr_agc) * M_SQRT1_2;
 200
 201    /*** CENTER FRONT CANCELLATION ***/
 202    /* A heuristic approach exploits that Lt + Rt gain contains the
 203       information about Lt, Rt correlation.  This effectively reshapes
 204       the front and rear "cones" to concentrate Lt + Rt to C and
 205       introduce Lt - Rt in L, R. */
 206    /* 0.67677 is the emprical lower bound for lpr_gain. */
 207    c_gain = 8 * (*adapt_lpr_gain - 0.67677);
 208    c_gain = c_gain > 0 ? c_gain : 0;
 209    /* c_gain should not be too high, not even reaching full
 210       cancellation (~ 0.50 - 0.55 at current AGC implementation), or
 211       the center will s0und too narrow. */
 212    c_gain = MATCOMPGAIN / (1 + c_gain * c_gain);
 213    c_agc_cfk = c_gain * cf[k];
 214    lf[k] -= c_agc_cfk;
 215    rf[k] -= c_agc_cfk;
 216    cf[k] += c_agc_cfk + c_agc_cfk;
 217 #if 0
 218    if(counter % 240 == 0)
 219       fprintf(fp_out, "%g %g %g %g %g\n",
 220               *adapt_l_gain, *adapt_r_gain,
 221               *adapt_lpr_gain, *adapt_lmr_gain,
 222               c_gain);
 223    counter++;
 224 #endif
 225 }
 226
 227 static inline void update_ch(af_hrtf_t *s, short *in, const int k)
 228 {
 229     const int fwr_pos = (k + FWRDURATION) % s->dlbuflen;
 230     /* Update the full wave rectified total amplitude */
 231     /* Input matrix decoder */
 232     if(s->decode_mode == HRTF_MIX_MATRIX2CH) {
 233        s->l_fwr += abs(in[0]) - fabs(s->fwrbuf_l[fwr_pos]);
 234        s->r_fwr += abs(in[1]) - fabs(s->fwrbuf_r[fwr_pos]);
 235        s->lpr_fwr += abs(in[0] + in[1]) -
 236           fabs(s->fwrbuf_l[fwr_pos] + s->fwrbuf_r[fwr_pos]);
 237        s->lmr_fwr += abs(in[0] - in[1]) -
 238           fabs(s->fwrbuf_l[fwr_pos] - s->fwrbuf_r[fwr_pos]);
 239     }
 240     /* Rear matrix decoder */
 241     if(s->matrix_mode) {
 242        s->lr_fwr += abs(in[2]) - fabs(s->fwrbuf_lr[fwr_pos]);
 243        s->rr_fwr += abs(in[3]) - fabs(s->fwrbuf_rr[fwr_pos]);
 244        s->lrprr_fwr += abs(in[2] + in[3]) -
 245           fabs(s->fwrbuf_lr[fwr_pos] + s->fwrbuf_rr[fwr_pos]);
 246        s->lrmrr_fwr += abs(in[2] - in[3]) -
 247           fabs(s->fwrbuf_lr[fwr_pos] - s->fwrbuf_rr[fwr_pos]);
 248     }
 249
 250     switch (s->decode_mode) {
 251     case HRTF_MIX_51:
 252        /* 5/5+1 channel sources */
 253        s->lf[k] = in[0];
 254        s->cf[k] = in[4];
 255        s->rf[k] = in[1];
 256        s->fwrbuf_lr[k] = s->lr[k] = in[2];
 257        s->fwrbuf_rr[k] = s->rr[k] = in[3];
 258        break;
 259     case HRTF_MIX_MATRIX2CH:
 260        /* Matrix encoded 2 channel sources */
 261        s->fwrbuf_l[k] = in[0];
 262        s->fwrbuf_r[k] = in[1];
 263        matrix_decode(in, k, 0, 1, 1, s->dlbuflen,
 264                      s->l_fwr, s->r_fwr,
 265                      s->lpr_fwr, s->lmr_fwr,
 266                      &(s->adapt_l_gain), &(s->adapt_r_gain),
 267                      &(s->adapt_lpr_gain), &(s->adapt_lmr_gain),
 268                      s->lf, s->rf, s->lr, s->rr, s->cf);
 269        break;
 270     case HRTF_MIX_STEREO:
 271        /* Stereo sources */
 272        s->lf[k] = in[0];
 273        s->rf[k] = in[1];
 274        s->cf[k] = s->lr[k] = s->rr[k] = 0;
 275        break;
 276     }
 277
 278     /* We need to update the bass compensation delay line, too. */
 279     s->ba_l[k] = in[0] + in[4] + in[2];
 280     s->ba_r[k] = in[4] + in[1] + in[3];
 281 }
 282
 283 /* Initialization and runtime control */
 284 static int control(struct af_instance_s *af, int cmd, void* arg)
 285 {
 286     af_hrtf_t *s = af->setup;
 287     int test_output_res;
 288     char mode;
 289
 290     switch(cmd) {
 291     case AF_CONTROL_REINIT:
 292         af->data->rate   = ((af_data_t*)arg)->rate;
 293         if(af->data->rate != 48000) {
 294             // automatic samplerate adjustment in the filter chain
 295             // is not yet supported.
 296             mp_msg(MSGT_AFILTER, MSGL_ERR,
 297                    "[hrtf] ERROR: Sampling rate is not 48000 Hz (%d)!\n",
 298                    af->data->rate);
 299             return AF_ERROR;
 300         }
 301         af->data->nch    = ((af_data_t*)arg)->nch;
 302             if(af->data->nch == 2) {
 303                /* 2 channel input */
 304                if(s->decode_mode != HRTF_MIX_MATRIX2CH) {
 305                   /* Default behavior is stereo mixing. */
 306                   s->decode_mode = HRTF_MIX_STEREO;
 307                }
 308             }
 309             else if (af->data->nch < 5)
 310               af->data->nch = 5;
 311         af->data->format = AF_FORMAT_S16_NE;
 312         af->data->bps    = 2;
 313         test_output_res = af_test_output(af, (af_data_t*)arg);
 314         af->mul = 2.0 / af->data->nch;
 315         // after testing input set the real output format
 316         af->data->nch = 2;
 317         s->print_flag = 1;
 318         return test_output_res;
 319     case AF_CONTROL_COMMAND_LINE:
 320         sscanf((char*)arg, "%c", &mode);
 321         switch(mode) {
 322         case 'm':
 323             /* Use matrix rear decoding. */
 324             s->matrix_mode = 1;
 325             break;
 326         case 's':
 327             /* Input needs matrix decoding. */
 328             s->decode_mode = HRTF_MIX_MATRIX2CH;
 329             break;
 330         case '0':
 331             s->matrix_mode = 0;
 332             break;
 333         default:
 334             mp_msg(MSGT_AFILTER, MSGL_ERR,
 335                    "[hrtf] Mode is neither 'm', 's', nor '0' (%c).\n",
 336                    mode);
 337             return AF_ERROR;
 338         }
 339         s->print_flag = 1;
 340         return AF_OK;
 341     }
 342
 343     return AF_UNKNOWN;
 344 }
 345
 346 /* Deallocate memory */
 347 static void uninit(struct af_instance_s *af)
 348 {
 349     if(af->setup) {
 350         af_hrtf_t *s = af->setup;
 351
 352         free(s->lf);
 353         free(s->rf);
 354         free(s->lr);
 355         free(s->rr);
 356         free(s->cf);
 357         free(s->cr);
 358         free(s->ba_l);
 359         free(s->ba_r);
 360         free(s->ba_ir);
 361         free(s->fwrbuf_l);
 362         free(s->fwrbuf_r);
 363         free(s->fwrbuf_lr);
 364         free(s->fwrbuf_rr);
 365         free(af->setup);
 366     }
 367     if(af->data)
 368         free(af->data->audio);
 369     free(af->data);
 370 }
 371
 372 /* Filter data through filter
 373
 374 Two "tricks" are used to compensate the "color" of the KEMAR data:
 375
 376 1. The KEMAR data is refiltered to ensure that the front L, R channels
 377 on the same side of the ear are equalized (especially in the high
 378 frequencies).
 379
 380 2. A bass compensation is introduced to ensure that 0-200 Hz are not
 381 damped (without any real 3D acoustical image, however).
 382 */
 383 static af_data_t* play(struct af_instance_s *af, af_data_t *data)
 384 {
 385     af_hrtf_t *s = af->setup;
 386     short *in = data->audio; // Input audio data
 387     short *out = NULL; // Output audio data
 388     short *end = in + data->len / sizeof(short); // Loop end
 389     float common, left, right, diff, left_b, right_b;
 390     const int dblen = s->dlbuflen, hlen = s->hrflen, blen = s->basslen;
 391
 392     if(AF_OK != RESIZE_LOCAL_BUFFER(af, data))
 393         return NULL;
 394
 395     if(s->print_flag) {
 396         s->print_flag = 0;
 397         switch (s->decode_mode) {
 398         case HRTF_MIX_51:
 399           mp_msg(MSGT_AFILTER, MSGL_INFO,
 400                  "[hrtf] Using HRTF to mix %s discrete surround into "
 401                  "L, R channels\n", s->matrix_mode ? "5+1" : "5");
 402           break;
 403         case HRTF_MIX_STEREO:
 404           mp_msg(MSGT_AFILTER, MSGL_INFO,
 405                  "[hrtf] Using HRTF to mix stereo into "
 406                  "L, R channels\n");
 407           break;
 408         case HRTF_MIX_MATRIX2CH:
 409           mp_msg(MSGT_AFILTER, MSGL_INFO,
 410                  "[hrtf] Using active matrix to decode 2 channel "
 411                  "input, HRTF to mix %s matrix surround into "
 412                  "L, R channels\n", "3/2");
 413           break;
 414         default:
 415           mp_msg(MSGT_AFILTER, MSGL_WARN,
 416                  "[hrtf] bogus decode_mode: %d\n", s->decode_mode);
 417           break;
 418         }
 419
 420        if(s->matrix_mode)
 421           mp_msg(MSGT_AFILTER, MSGL_INFO,
 422                  "[hrtf] Using active matrix to decode rear center "
 423                  "channel\n");
 424     }
 425
 426     out = af->data->audio;
 427
 428     /* MPlayer's 5 channel layout (notation for the variable):
 429      *
 430      * 0: L (LF), 1: R (RF), 2: Ls (LR), 3: Rs (RR), 4: C (CF), matrix
 431      * encoded: Cs (CR)
 432      *
 433      * or: L = left, C = center, R = right, F = front, R = rear
 434      *
 435      * Filter notation:
 436      *
 437      *      CF
 438      * OF        AF
 439      *      Ear->
 440      * OR        AR
 441      *      CR
 442      *
 443      * or: C = center, A = same side, O = opposite, F = front, R = rear
 444      */
 445
 446     while(in < end) {
 447         const int k = s->cyc_pos;
 448
 449         update_ch(s, in, k);
 450
 451         /* Simulate a 7.5 ms -20 dB echo of the center channel in the
 452            front channels (like reflection from a room wall) - a kind of
 453            psycho-acoustically "cheating" to focus the center front
 454            channel, which is normally hard to be perceived as front */
 455         s->lf[k] += CFECHOAMPL * s->cf[(k + CFECHODELAY) % s->dlbuflen];
 456         s->rf[k] += CFECHOAMPL * s->cf[(k + CFECHODELAY) % s->dlbuflen];
 457
 458         switch (s->decode_mode) {
 459         case HRTF_MIX_51:
 460         case HRTF_MIX_MATRIX2CH:
 461            /* Mixer filter matrix */
 462            common = conv(dblen, hlen, s->cf, s->cf_ir, k + s->cf_o);
 463            if(s->matrix_mode) {
 464               /* In matrix decoding mode, the rear channel gain must be
 465                  renormalized, as there is an additional channel. */
 466               matrix_decode(in, k, 2, 3, 0, s->dlbuflen,
 467                             s->lr_fwr, s->rr_fwr,
 468                             s->lrprr_fwr, s->lrmrr_fwr,
 469                             &(s->adapt_lr_gain), &(s->adapt_rr_gain),
 470                             &(s->adapt_lrprr_gain), &(s->adapt_lrmrr_gain),
 471                             s->lr, s->rr, NULL, NULL, s->cr);
 472               common +=
 473                  conv(dblen, hlen, s->cr, s->cr_ir, k + s->cr_o) *
 474                  M1_76DB;
 475               left    =
 476                  ( conv(dblen, hlen, s->lf, s->af_ir, k + s->af_o) +
 477                    conv(dblen, hlen, s->rf, s->of_ir, k + s->of_o) +
 478                    (conv(dblen, hlen, s->lr, s->ar_ir, k + s->ar_o) +
 479                     conv(dblen, hlen, s->rr, s->or_ir, k + s->or_o)) *
 480                    M1_76DB + common);
 481               right   =
 482                  ( conv(dblen, hlen, s->rf, s->af_ir, k + s->af_o) +
 483                    conv(dblen, hlen, s->lf, s->of_ir, k + s->of_o) +
 484                    (conv(dblen, hlen, s->rr, s->ar_ir, k + s->ar_o) +
 485                     conv(dblen, hlen, s->lr, s->or_ir, k + s->or_o)) *
 486                    M1_76DB + common);
 487            } else {
 488               left    =
 489                  ( conv(dblen, hlen, s->lf, s->af_ir, k + s->af_o) +
 490                    conv(dblen, hlen, s->rf, s->of_ir, k + s->of_o) +
 491                    conv(dblen, hlen, s->lr, s->ar_ir, k + s->ar_o) +
 492                    conv(dblen, hlen, s->rr, s->or_ir, k + s->or_o) +
 493                    common);
 494               right   =
 495                  ( conv(dblen, hlen, s->rf, s->af_ir, k + s->af_o) +
 496                    conv(dblen, hlen, s->lf, s->of_ir, k + s->of_o) +
 497                    conv(dblen, hlen, s->rr, s->ar_ir, k + s->ar_o) +
 498                    conv(dblen, hlen, s->lr, s->or_ir, k + s->or_o) +
 499                    common);
 500            }
 501            break;
 502         case HRTF_MIX_STEREO:
 503            left    =
 504               ( conv(dblen, hlen, s->lf, s->af_ir, k + s->af_o) +
 505                 conv(dblen, hlen, s->rf, s->of_ir, k + s->of_o));
 506            right   =
 507               ( conv(dblen, hlen, s->rf, s->af_ir, k + s->af_o) +
 508                 conv(dblen, hlen, s->lf, s->of_ir, k + s->of_o));
 509            break;
 510         default:
 511             /* make gcc happy */
 512             left = 0.0;
 513             right = 0.0;
 514             break;
 515         }
 516
 517         /* Bass compensation for the lower frequency cut of the HRTF.  A
 518            cross talk of the left and right channel is introduced to
 519            match the directional characteristics of higher frequencies.
 520            The bass will not have any real 3D perception, but that is
 521            OK (note at 180 Hz, the wavelength is about 2 m, and any
 522            spatial perception is impossible). */
 523         left_b  = conv(dblen, blen, s->ba_l, s->ba_ir, k);
 524         right_b = conv(dblen, blen, s->ba_r, s->ba_ir, k);
 525         left  += (1 - BASSCROSS) * left_b  + BASSCROSS * right_b;
 526         right += (1 - BASSCROSS) * right_b + BASSCROSS * left_b;
 527         /* Also mix the LFE channel (if available) */
 528         if(data->nch >= 6) {
 529             left  += in[5] * M3_01DB;
 530             right += in[5] * M3_01DB;
 531         }
 532
 533         /* Amplitude renormalization. */
 534         left  *= AMPLNORM;
 535         right *= AMPLNORM;
 536
 537         switch (s->decode_mode) {
 538         case HRTF_MIX_51:
 539         case HRTF_MIX_STEREO:
 540            /* "Cheating": linear stereo expansion to amplify the 3D
 541               perception.  Note: Too much will destroy the acoustic space
 542               and may even result in headaches. */
 543            diff = STEXPAND2 * (left - right);
 544            out[0] = (int16_t)(left  + diff);
 545            out[1] = (int16_t)(right - diff);
 546            break;
 547         case HRTF_MIX_MATRIX2CH:
 548            /* Do attempt any stereo expansion with matrix encoded
 549               sources.  The L, R channels are already stereo expanded
 550               by the steering, any further stereo expansion will sound
 551               very unnatural. */
 552            out[0] = (int16_t)left;
 553            out[1] = (int16_t)right;
 554            break;
 555         }
 556
 557         /* Next sample... */
 558         in = &in[data->nch];
 559         out = &out[af->data->nch];
 560         (s->cyc_pos)--;
 561         if(s->cyc_pos < 0)
 562             s->cyc_pos += dblen;
 563     }
 564
 565     /* Set output data */
 566     data->audio = af->data->audio;
 567     data->len   = data->len / data->nch * 2;
 568     data->nch   = 2;
 569
 570     return data;
 571 }
 572
 573 static int allocate(af_hrtf_t *s)
 574 {
 575     if ((s->lf = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1;
 576     if ((s->rf = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1;
 577     if ((s->lr = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1;
 578     if ((s->rr = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1;
 579     if ((s->cf = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1;
 580     if ((s->cr = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1;
 581     if ((s->ba_l = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1;
 582     if ((s->ba_r = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1;
 583     if ((s->fwrbuf_l =
 584          malloc(s->dlbuflen * sizeof(float))) == NULL) return -1;
 585     if ((s->fwrbuf_r =
 586          malloc(s->dlbuflen * sizeof(float))) == NULL) return -1;
 587     if ((s->fwrbuf_lr =
 588          malloc(s->dlbuflen * sizeof(float))) == NULL) return -1;
 589     if ((s->fwrbuf_rr =
 590          malloc(s->dlbuflen * sizeof(float))) == NULL) return -1;
 591     return 0;
 592 }
 593
 594 /* Allocate memory and set function pointers */
 595 static int af_open(af_instance_t* af)
 596 {
 597     int i;
 598     af_hrtf_t *s;
 599     float fc;
 600
 601     af->control = control;
 602     af->uninit = uninit;
 603     af->play = play;
 604     af->mul = 1;
 605     af->data = calloc(1, sizeof(af_data_t));
 606     af->setup = calloc(1, sizeof(af_hrtf_t));
 607     if((af->data == NULL) || (af->setup == NULL))
 608         return AF_ERROR;
 609
 610     s = af->setup;
 611
 612     s->dlbuflen = DELAYBUFLEN;
 613     s->hrflen = HRTFFILTLEN;
 614     s->basslen = BASSFILTLEN;
 615
 616     s->cyc_pos = s->dlbuflen - 1;
 617     /* With a full (two axis) steering matrix decoder, s->matrix_mode
 618        should not be enabled lightly (it will also steer the Ls, Rs
 619        channels). */
 620     s->matrix_mode = 0;
 621     s->decode_mode = HRTF_MIX_51;
 622
 623     s->print_flag = 1;
 624
 625     if (allocate(s) != 0) {
 626         mp_msg(MSGT_AFILTER, MSGL_ERR, "[hrtf] Memory allocation error.\n");
 627         return AF_ERROR;
 628     }
 629
 630     for(i = 0; i < s->dlbuflen; i++)
 631         s->lf[i] = s->rf[i] = s->lr[i] = s->rr[i] = s->cf[i] =
 632             s->cr[i] = 0;
 633
 634     s->lr_fwr =
 635         s->rr_fwr = 0;
 636
 637     s->cf_ir = cf_filt + (s->cf_o = pulse_detect(cf_filt));
 638     s->af_ir = af_filt + (s->af_o = pulse_detect(af_filt));
 639     s->of_ir = of_filt + (s->of_o = pulse_detect(of_filt));
 640     s->ar_ir = ar_filt + (s->ar_o = pulse_detect(ar_filt));
 641     s->or_ir = or_filt + (s->or_o = pulse_detect(or_filt));
 642     s->cr_ir = cr_filt + (s->cr_o = pulse_detect(cr_filt));
 643
 644     if((s->ba_ir = malloc(s->basslen * sizeof(float))) == NULL) {
 645         mp_msg(MSGT_AFILTER, MSGL_ERR, "[hrtf] Memory allocation error.\n");
 646         return AF_ERROR;
 647     }
 648     fc = 2.0 * BASSFILTFREQ / (float)af->data->rate;
 649     if(af_filter_design_fir(s->basslen, s->ba_ir, &fc, LP | KAISER, 4 * M_PI) ==
 650        -1) {
 651         mp_msg(MSGT_AFILTER, MSGL_ERR, "[hrtf] Unable to design low-pass "
 652                "filter.\n");
 653         return AF_ERROR;
 654     }
 655     for(i = 0; i < s->basslen; i++)
 656         s->ba_ir[i] *= BASSGAIN;
 657
 658     return AF_OK;
 659 }
 660
 661 /* Description of this filter */
 662 af_info_t af_info_hrtf = {
 663     "HRTF Headphone",
 664     "hrtf",
 665     "ylai",
 666     "",
 667     AF_FLAGS_REENTRANT,
 668     af_open
 669 };