Alc/nfcfilter.c

   1
   2 #include "config.h"
   3
   4 #include "nfcfilter.h"
   5
   6 #include "alu.h"
   7
   8
   9 /* Near-field control filters are the basis for handling the near-field effect.
  10  * The near-field effect is a bass-boost present in the directional components
  11  * of a recorded signal, created as a result of the wavefront curvature (itself
  12  * a function of sound distance). Proper reproduction dictates this be
  13  * compensated for using a bass-cut given the playback speaker distance, to
  14  * avoid excessive bass in the playback.
  15  *
  16  * For real-time rendered audio, emulating the near-field effect based on the
  17  * sound source's distance, and subsequently compensating for it at output
  18  * based on the speaker distances, can create a more realistic perception of
  19  * sound distance beyond a simple 1/r attenuation.
  20  *
  21  * These filters do just that. Each one applies a low-shelf filter, created as
  22  * the combination of a bass-boost for a given sound source distance (near-
  23  * field emulation) along with a bass-cut for a given control/speaker distance
  24  * (near-field compensation).
  25  *
  26  * Note that it is necessary to apply a cut along with the boost, since the
  27  * boost alone is unstable in higher-order ambisonics as it causes an infinite
  28  * DC gain (even first-order ambisonics requires there to be no DC offset for
  29  * the boost to work). Consequently, ambisonics requires a control parameter to
  30  * be used to avoid an unstable boost-only filter. NFC-HOA defines this control
  31  * as a reference delay, calculated with:
  32  *
  33  * reference_delay = control_distance / speed_of_sound
  34  *
  35  * This means w0 (for input) or w1 (for output) should be set to:
  36  *
  37  * wN = 1 / (reference_delay * sample_rate)
  38  *
  39  * when dealing with NFC-HOA content. For FOA input content, which does not
  40  * specify a reference_delay variable, w0 should be set to 0 to apply only
  41  * near-field compensation for output. It's important that w1 be a finite,
  42  * positive, non-0 value or else the bass-boost will become unstable again.
  43  * Also, w0 should not be too large compared to w1, to avoid excessively loud
  44  * low frequencies.
  45  */
  46
  47 static const float B[4][3] = {
  48     {    0.0f                             },
  49     {    1.0f                             },
  50     {    3.0f,     3.0f                   },
  51     { 3.6778f,  6.4595f, 2.3222f          },
  52   /*{ 4.2076f, 11.4877f, 5.7924f, 9.1401f }*/
  53 };
  54
  55 static void NfcFilterCreate1(struct NfcFilter1 *nfc, const float w0, const float w1)
  56 {
  57     float b_00, g_0;
  58     float r;
  59
  60     nfc->g = 1.0f;
  61     nfc->coeffs[0] = 1.0f;
  62
  63     /* Calculate bass-boost coefficients. */
  64     r = 0.5f * w0;
  65     b_00 = B[1][0] * r;
  66     g_0 = 1.0f + b_00;
  67
  68     nfc->coeffs[0] *= g_0;
  69     nfc->coeffs[1] = (2.0f * b_00) / g_0;
  70
  71     /* Calculate bass-cut coefficients. */
  72     r = 0.5f * w1;
  73     b_00 = B[1][0] * r;
  74     g_0 = 1.0f + b_00;
  75
  76     nfc->g /= g_0;
  77     nfc->coeffs[0] /= g_0;
  78     nfc->coeffs[1+1] = (2.0f * b_00) / g_0;
  79 }
  80
  81 static void NfcFilterAdjust1(struct NfcFilter1 *nfc, const float w0)
  82 {
  83     float b_00, g_0;
  84     float r;
  85
  86     r = 0.5f * w0;
  87     b_00 = B[1][0] * r;
  88     g_0 = 1.0f + b_00;
  89
  90     nfc->coeffs[0] = nfc->g * g_0;
  91     nfc->coeffs[1] = (2.0f * b_00) / g_0;
  92 }
  93
  94
  95 static void NfcFilterCreate2(struct NfcFilter2 *nfc, const float w0, const float w1)
  96 {
  97     float b_10, b_11, g_1;
  98     float r;
  99
 100     nfc->g = 1.0f;
 101     nfc->coeffs[0] = 1.0f;
 102
 103     /* Calculate bass-boost coefficients. */
 104     r = 0.5f * w0;
 105     b_10 = B[2][0] * r;
 106     b_11 = B[2][1] * r * r;
 107     g_1 = 1.0f + b_10 + b_11;
 108
 109     nfc->coeffs[0] *= g_1;
 110     nfc->coeffs[1] = ((2.0f * b_10) + (4.0f * b_11)) / g_1;
 111     nfc->coeffs[2] = (4.0f * b_11) / g_1;
 112
 113     /* Calculate bass-cut coefficients. */
 114     r = 0.5f * w1;
 115     b_10 = B[2][0] * r;
 116     b_11 = B[2][1] * r * r;
 117     g_1 = 1.0f + b_10 + b_11;
 118
 119     nfc->g /= g_1;
 120     nfc->coeffs[0] /= g_1;
 121     nfc->coeffs[2+1] = ((2.0f * b_10) + (4.0f * b_11)) / g_1;
 122     nfc->coeffs[2+2] = (4.0f * b_11) / g_1;
 123 }
 124
 125 static void NfcFilterAdjust2(struct NfcFilter2 *nfc, const float w0)
 126 {
 127     float b_10, b_11, g_1;
 128     float r;
 129
 130     r = 0.5f * w0;
 131     b_10 = B[2][0] * r;
 132     b_11 = B[2][1] * r * r;
 133     g_1 = 1.0f + b_10 + b_11;
 134
 135     nfc->coeffs[0] = nfc->g * g_1;
 136     nfc->coeffs[1] = ((2.0f * b_10) + (4.0f * b_11)) / g_1;
 137     nfc->coeffs[2] = (4.0f * b_11) / g_1;
 138 }
 139
 140
 141 static void NfcFilterCreate3(struct NfcFilter3 *nfc, const float w0, const float w1)
 142 {
 143     float b_10, b_11, g_1;
 144     float b_00, g_0;
 145     float r;
 146
 147     nfc->g = 1.0f;
 148     nfc->coeffs[0] = 1.0f;
 149
 150     /* Calculate bass-boost coefficients. */
 151     r = 0.5f * w0;
 152     b_10 = B[3][0] * r;
 153     b_11 = B[3][1] * r * r;
 154     g_1 = 1.0f + b_10 + b_11;
 155
 156     nfc->coeffs[0] *= g_1;
 157     nfc->coeffs[1] = ((2.0f * b_10) + (4.0f * b_11)) / g_1;
 158     nfc->coeffs[2] = (4.0f * b_11) / g_1;
 159
 160     b_00 = B[3][2] * r;
 161     g_0 = 1.0f + b_00;
 162
 163     nfc->coeffs[0] *= g_0;
 164     nfc->coeffs[2+1] = (2.0f * b_00) / g_0;
 165
 166     /* Calculate bass-cut coefficients. */
 167     r = 0.5f * w1;
 168     b_10 = B[3][0] * r;
 169     b_11 = B[3][1] * r * r;
 170     g_1 = 1.0f + b_10 + b_11;
 171
 172     nfc->g /= g_1;
 173     nfc->coeffs[0] /= g_1;
 174     nfc->coeffs[3+1] = ((2.0f * b_10) + (4.0f * b_11)) / g_1;
 175     nfc->coeffs[3+2] = (4.0f * b_11) / g_1;
 176
 177     b_00 = B[3][2] * r;
 178     g_0 = 1.0f + b_00;
 179
 180     nfc->g /= g_0;
 181     nfc->coeffs[0] /= g_0;
 182     nfc->coeffs[3+2+1] = (2.0f * b_00) / g_0;
 183 }
 184
 185 static void NfcFilterAdjust3(struct NfcFilter3 *nfc, const float w0)
 186 {
 187     float b_10, b_11, g_1;
 188     float b_00, g_0;
 189     float r;
 190
 191     r = 0.5f * w0;
 192     b_10 = B[3][0] * r;
 193     b_11 = B[3][1] * r * r;
 194     g_1 = 1.0f + b_10 + b_11;
 195
 196     nfc->coeffs[0] = nfc->g * g_1;
 197     nfc->coeffs[1] = ((2.0f * b_10) + (4.0f * b_11)) / g_1;
 198     nfc->coeffs[2] = (4.0f * b_11) / g_1;
 199
 200     b_00 = B[3][2] * r;
 201     g_0 = 1.0f + b_00;
 202
 203     nfc->coeffs[0] *= g_0;
 204     nfc->coeffs[2+1] = (2.0f * b_00) / g_0;
 205 }
 206
 207
 208 void NfcFilterCreate(NfcFilter *nfc, const float w0, const float w1)
 209 {
 210     memset(nfc, 0, sizeof(*nfc));
 211     NfcFilterCreate1(&nfc->first, w0, w1);
 212     NfcFilterCreate2(&nfc->second, w0, w1);
 213     NfcFilterCreate3(&nfc->third, w0, w1);
 214 }
 215
 216 void NfcFilterAdjust(NfcFilter *nfc, const float w0)
 217 {
 218     NfcFilterAdjust1(&nfc->first, w0);
 219     NfcFilterAdjust2(&nfc->second, w0);
 220     NfcFilterAdjust3(&nfc->third, w0);
 221 }
 222
 223
 224 void NfcFilterUpdate1(NfcFilter *nfc, ALfloat *restrict dst, const float *restrict src, const int count)
 225 {
 226     const float b0 = nfc->first.coeffs[0];
 227     const float a0 = nfc->first.coeffs[1];
 228     const float a1 = nfc->first.coeffs[2];
 229     float z1 = nfc->first.history[0];
 230     int i;
 231
 232     for(i = 0;i < count;i++)
 233     {
 234         float out = src[i] * b0;
 235         float y;
 236
 237         y = out - (a1*z1);
 238         out = y + (a0*z1);
 239         z1 += y;
 240
 241         dst[i] = out;
 242     }
 243     nfc->first.history[0] = z1;
 244 }
 245
 246 void NfcFilterUpdate2(NfcFilter *nfc, ALfloat *restrict dst, const float *restrict src, const int count)
 247 {
 248     const float b0 = nfc->second.coeffs[0];
 249     const float a00 = nfc->second.coeffs[1];
 250     const float a01 = nfc->second.coeffs[2];
 251     const float a10 = nfc->second.coeffs[3];
 252     const float a11 = nfc->second.coeffs[4];
 253     float z1 = nfc->second.history[0];
 254     float z2 = nfc->second.history[1];
 255     int i;
 256
 257     for(i = 0;i < count;i++)
 258     {
 259         float out = src[i] * b0;
 260         float y;
 261
 262         y = out - (a10*z1) - (a11*z2);
 263         out = y + (a00*z1) + (a01*z2);
 264         z2 += z1;
 265         z1 += y;
 266
 267         dst[i] = out;
 268     }
 269     nfc->second.history[0] = z1;
 270     nfc->second.history[1] = z2;
 271 }
 272
 273 void NfcFilterUpdate3(NfcFilter *nfc, ALfloat *restrict dst, const float *restrict src, const int count)
 274 {
 275     const float b0 = nfc->third.coeffs[0];
 276     const float a00 = nfc->third.coeffs[1];
 277     const float a01 = nfc->third.coeffs[2];
 278     const float a02 = nfc->third.coeffs[3];
 279     const float a10 = nfc->third.coeffs[4];
 280     const float a11 = nfc->third.coeffs[5];
 281     const float a12 = nfc->third.coeffs[6];
 282     float z1 = nfc->third.history[0];
 283     float z2 = nfc->third.history[1];
 284     float z3 = nfc->third.history[2];
 285     int i;
 286
 287     for(i = 0;i < count;i++)
 288     {
 289         float out = src[i] * b0;
 290         float y;
 291
 292         y = out - (a10*z1) - (a11*z2);
 293         out = y + (a00*z1) + (a01*z2);
 294         z2 += z1;
 295         z1 += y;
 296
 297         y = out - (a12*z3);
 298         out = y + (a02*z3);
 299         z3 += y;
 300
 301         dst[i] = out;
 302     }
 303     nfc->third.history[0] = z1;
 304     nfc->third.history[1] = z2;
 305     nfc->third.history[2] = z3;
 306 }
 307
 308 #if 0 /* Original methods the above are derived from. */
 309 static void NfcFilterCreate(NfcFilter *nfc, const ALsizei order, const float src_dist, const float ctl_dist, const float rate)
 310 {
 311     static const float B[4][5] = {
 312         {                                     },
 313         {    1.0f                             },
 314         {    3.0f,     3.0f                   },
 315         { 3.6778f,  6.4595f, 2.3222f          },
 316         { 4.2076f, 11.4877f, 5.7924f, 9.1401f }
 317     };
 318     float w0 = SPEEDOFSOUNDMETRESPERSEC / (src_dist * rate);
 319     float w1 = SPEEDOFSOUNDMETRESPERSEC / (ctl_dist * rate);
 320     ALsizei i;
 321     float r;
 322
 323     nfc->g = 1.0f;
 324     nfc->coeffs[0] = 1.0f;
 325
 326     /* NOTE: Slight adjustment from the literature to raise the center
 327      * frequency a bit (0.5 -> 1.0).
 328      */
 329     r = 1.0f * w0;
 330     for(i = 0; i < (order-1);i += 2)
 331     {
 332         float b_10 = B[order][i  ] * r;
 333         float b_11 = B[order][i+1] * r * r;
 334         float g_1 = 1.0f + b_10 + b_11;
 335
 336         nfc->b[i] = b_10;
 337         nfc->b[i + 1] = b_11;
 338         nfc->coeffs[0] *= g_1;
 339         nfc->coeffs[i+1] = ((2.0f * b_10) + (4.0f * b_11)) / g_1;
 340         nfc->coeffs[i+2] = (4.0f * b_11) / g_1;
 341     }
 342     if(i < order)
 343     {
 344         float b_00 = B[order][i] * r;
 345         float g_0 = 1.0f + b_00;
 346
 347         nfc->b[i] = b_00;
 348         nfc->coeffs[0] *= g_0;
 349         nfc->coeffs[i+1] = (2.0f * b_00) / g_0;
 350     }
 351
 352     r = 1.0f * w1;
 353     for(i = 0;i < (order-1);i += 2)
 354     {
 355         float b_10 = B[order][i  ] * r;
 356         float b_11 = B[order][i+1] * r * r;
 357         float g_1 = 1.0f + b_10 + b_11;
 358
 359         nfc->g /= g_1;
 360         nfc->coeffs[0] /= g_1;
 361         nfc->coeffs[order+i+1] = ((2.0f * b_10) + (4.0f * b_11)) / g_1;
 362         nfc->coeffs[order+i+2] = (4.0f * b_11) / g_1;
 363     }
 364     if(i < order)
 365     {
 366         float b_00 = B[order][i] * r;
 367         float g_0 = 1.0f + b_00;
 368
 369         nfc->g /= g_0;
 370         nfc->coeffs[0] /= g_0;
 371         nfc->coeffs[order+i+1] = (2.0f * b_00) / g_0;
 372     }
 373
 374     for(i = 0; i < MAX_AMBI_ORDER; i++)
 375         nfc->history[i] = 0.0f;
 376 }
 377
 378 static void NfcFilterAdjust(NfcFilter *nfc, const float distance)
 379 {
 380     int i;
 381
 382     nfc->coeffs[0] = nfc->g;
 383
 384     for(i = 0;i < (nfc->order-1);i += 2)
 385     {
 386         float b_10 = nfc->b[i] / distance;
 387         float b_11 = nfc->b[i+1] / (distance * distance);
 388         float g_1 = 1.0f + b_10 + b_11;
 389
 390         nfc->coeffs[0] *= g_1;
 391         nfc->coeffs[i+1] = ((2.0f * b_10) + (4.0f * b_11)) / g_1;
 392         nfc->coeffs[i+2] = (4.0f * b_11) / g_1;
 393     }
 394     if(i < nfc->order)
 395     {
 396         float b_00 = nfc->b[i] / distance;
 397         float g_0 = 1.0f + b_00;
 398
 399         nfc->coeffs[0] *= g_0;
 400         nfc->coeffs[i+1] = (2.0f * b_00) / g_0;
 401     }
 402 }
 403
 404 static float NfcFilterUpdate(const float in, NfcFilter *nfc)
 405 {
 406     int i;
 407     float out = in * nfc->coeffs[0];
 408
 409     for(i = 0;i < (nfc->order-1);i += 2)
 410     {
 411         float y = out - (nfc->coeffs[nfc->order+i+1] * nfc->history[i]) -
 412                         (nfc->coeffs[nfc->order+i+2] * nfc->history[i+1]) + 1.0e-30f;
 413         out = y + (nfc->coeffs[i+1]*nfc->history[i]) + (nfc->coeffs[i+2]*nfc->history[i+1]);
 414
 415         nfc->history[i+1] += nfc->history[i];
 416         nfc->history[i] += y;
 417     }
 418     if(i < nfc->order)
 419     {
 420         float y = out - (nfc->coeffs[nfc->order+i+1] * nfc->history[i]) + 1.0e-30f;
 421
 422         out = y + (nfc->coeffs[i+1] * nfc->history[i]);
 423         nfc->history[i] += y;
 424     }
 425
 426     return out;
 427 }
 428 #endif