Alc/filters/nfc.c

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