libs/fluidsynth/src/fluid_rev.c

   1 /*
   2
   3   Freeverb
   4
   5   Written by Jezar at Dreampoint, June 2000
   6   http://www.dreampoint.co.uk
   7   This code is public domain
   8
   9   Translated to C by Peter Hanappe, Mai 2001
  10 */
  11
  12 #include "fluid_rev.h"
  13
  14 /***************************************************************
  15  *
  16  *                           REVERB
  17  */
  18
  19 /* Denormalising:
  20  *
  21  * According to music-dsp thread 'Denormalise', Pentium processors
  22  * have a hardware 'feature', that is of interest here, related to
  23  * numeric underflow.  We have a recursive filter. The output decays
  24  * exponentially, if the input stops.  So the numbers get smaller and
  25  * smaller... At some point, they reach 'denormal' level.  This will
  26  * lead to drastic spikes in the CPU load.  The effect was reproduced
  27  * with the reverb - sometimes the average load over 10 s doubles!!.
  28  *
  29  * The 'undenormalise' macro fixes the problem: As soon as the number
  30  * is close enough to denormal level, the macro forces the number to
  31  * 0.0f.  The original macro is:
  32  *
  33  * #define undenormalise(sample) if(((*(unsigned int*)&sample)&0x7f800000)==0) sample=0.0f
  34  *
  35  * This will zero out a number when it reaches the denormal level.
  36  * Advantage: Maximum dynamic range Disadvantage: We'll have to check
  37  * every sample, expensive.  The alternative macro comes from a later
  38  * mail from Jon Watte. It will zap a number before it reaches
  39  * denormal level. Jon suggests to run it once per block instead of
  40  * every sample.
  41  */
  42
  43 # if defined(WITH_FLOATX)
  44 # define zap_almost_zero(sample) (((*(unsigned int*)&(sample))&0x7f800000) < 0x08000000)?0.0f:(sample)
  45 # else
  46 /* 1e-20 was chosen as an arbitrary (small) threshold. */
  47 #define zap_almost_zero(sample) fabs(sample)<1e-10 ? 0 : sample;
  48 #endif
  49
  50 /* Denormalising part II:
  51  *
  52  * Another method fixes the problem cheaper: Use a small DC-offset in
  53  * the filter calculations.  Now the signals converge not against 0,
  54  * but against the offset.  The constant offset is invisible from the
  55  * outside world (i.e. it does not appear at the output.  There is a
  56  * very small turn-on transient response, which should not cause
  57  * problems.
  58  */
  59
  60
  61 //#define DC_OFFSET 0
  62 #define DC_OFFSET 1e-8
  63 //#define DC_OFFSET 0.001f
  64 typedef struct _fluid_allpass fluid_allpass;
  65 typedef struct _fluid_comb fluid_comb;
  66
  67 struct _fluid_allpass {
  68   fluid_real_t feedback;
  69   fluid_real_t *buffer;
  70   int bufsize;
  71   int bufidx;
  72 };
  73
  74 void fluid_allpass_init(fluid_allpass* allpass);
  75 void fluid_allpass_setfeedback(fluid_allpass* allpass, fluid_real_t val);
  76 fluid_real_t fluid_allpass_getfeedback(fluid_allpass* allpass);
  77
  78 static void
  79 fluid_allpass_setbuffer(fluid_allpass* allpass, int size)
  80 {
  81   allpass->bufidx = 0;
  82   allpass->buffer = FLUID_ARRAY(fluid_real_t,size);
  83   allpass->bufsize = size;
  84 }
  85
  86 static void
  87 fluid_allpass_release(fluid_allpass* allpass)
  88 {
  89   FLUID_FREE(allpass->buffer);
  90 }
  91
  92 void
  93 fluid_allpass_init(fluid_allpass* allpass)
  94 {
  95   int i;
  96   int len = allpass->bufsize;
  97   fluid_real_t* buf = allpass->buffer;
  98   for (i = 0; i < len; i++) {
  99     buf[i] = DC_OFFSET; /* this is not 100 % correct. */
 100   }
 101 }
 102
 103 void
 104 fluid_allpass_setfeedback(fluid_allpass* allpass, fluid_real_t val)
 105 {
 106   allpass->feedback = val;
 107 }
 108
 109 fluid_real_t
 110 fluid_allpass_getfeedback(fluid_allpass* allpass)
 111 {
 112   return allpass->feedback;
 113 }
 114
 115 #define fluid_allpass_process(_allpass, _input) \
 116 { \
 117   fluid_real_t output; \
 118   fluid_real_t bufout; \
 119   bufout = _allpass.buffer[_allpass.bufidx]; \
 120   output = bufout-_input; \
 121   _allpass.buffer[_allpass.bufidx] = _input + (bufout * _allpass.feedback); \
 122   if (++_allpass.bufidx >= _allpass.bufsize) { \
 123     _allpass.bufidx = 0; \
 124   } \
 125   _input = output; \
 126 }
 127
 128 /*  fluid_real_t fluid_allpass_process(fluid_allpass* allpass, fluid_real_t input) */
 129 /*  { */
 130 /*    fluid_real_t output; */
 131 /*    fluid_real_t bufout; */
 132 /*    bufout = allpass->buffer[allpass->bufidx]; */
 133 /*    undenormalise(bufout); */
 134 /*    output = -input + bufout; */
 135 /*    allpass->buffer[allpass->bufidx] = input + (bufout * allpass->feedback); */
 136 /*    if (++allpass->bufidx >= allpass->bufsize) { */
 137 /*      allpass->bufidx = 0; */
 138 /*    } */
 139 /*    return output; */
 140 /*  } */
 141
 142 struct _fluid_comb {
 143   fluid_real_t feedback;
 144   fluid_real_t filterstore;
 145   fluid_real_t damp1;
 146   fluid_real_t damp2;
 147   fluid_real_t *buffer;
 148   int bufsize;
 149   int bufidx;
 150 };
 151
 152 void fluid_comb_setbuffer(fluid_comb* comb, int size);
 153 void fluid_comb_release(fluid_comb* comb);
 154 void fluid_comb_init(fluid_comb* comb);
 155 void fluid_comb_setdamp(fluid_comb* comb, fluid_real_t val);
 156 fluid_real_t fluid_comb_getdamp(fluid_comb* comb);
 157 void fluid_comb_setfeedback(fluid_comb* comb, fluid_real_t val);
 158 fluid_real_t fluid_comb_getfeedback(fluid_comb* comb);
 159
 160 void
 161 fluid_comb_setbuffer(fluid_comb* comb, int size)
 162 {
 163   comb->filterstore = 0;
 164   comb->bufidx = 0;
 165   comb->buffer = FLUID_ARRAY(fluid_real_t,size);
 166   comb->bufsize = size;
 167 }
 168
 169 void
 170 fluid_comb_release(fluid_comb* comb)
 171 {
 172   FLUID_FREE(comb->buffer);
 173 }
 174
 175 void
 176 fluid_comb_init(fluid_comb* comb)
 177 {
 178   int i;
 179   fluid_real_t* buf = comb->buffer;
 180   int len = comb->bufsize;
 181   for (i = 0; i < len; i++) {
 182     buf[i] = DC_OFFSET; /* This is not 100 % correct. */
 183   }
 184 }
 185
 186 void
 187 fluid_comb_setdamp(fluid_comb* comb, fluid_real_t val)
 188 {
 189   comb->damp1 = val;
 190   comb->damp2 = 1 - val;
 191 }
 192
 193 fluid_real_t
 194 fluid_comb_getdamp(fluid_comb* comb)
 195 {
 196   return comb->damp1;
 197 }
 198
 199 void
 200 fluid_comb_setfeedback(fluid_comb* comb, fluid_real_t val)
 201 {
 202   comb->feedback = val;
 203 }
 204
 205 fluid_real_t
 206 fluid_comb_getfeedback(fluid_comb* comb)
 207 {
 208   return comb->feedback;
 209 }
 210
 211 #define fluid_comb_process(_comb, _input, _output) \
 212 { \
 213   fluid_real_t _tmp = _comb.buffer[_comb.bufidx]; \
 214   _comb.filterstore = (_tmp * _comb.damp2) + (_comb.filterstore * _comb.damp1); \
 215   _comb.buffer[_comb.bufidx] = _input + (_comb.filterstore * _comb.feedback); \
 216   if (++_comb.bufidx >= _comb.bufsize) { \
 217     _comb.bufidx = 0; \
 218   } \
 219   _output += _tmp; \
 220 }
 221
 222 /* fluid_real_t fluid_comb_process(fluid_comb* comb, fluid_real_t input) */
 223 /* { */
 224 /*    fluid_real_t output; */
 225
 226 /*    output = comb->buffer[comb->bufidx]; */
 227 /*    undenormalise(output); */
 228 /*    comb->filterstore = (output * comb->damp2) + (comb->filterstore * comb->damp1); */
 229 /*    undenormalise(comb->filterstore); */
 230 /*    comb->buffer[comb->bufidx] = input + (comb->filterstore * comb->feedback); */
 231 /*    if (++comb->bufidx >= comb->bufsize) { */
 232 /*      comb->bufidx = 0; */
 233 /*    } */
 234
 235 /*    return output; */
 236 /* } */
 237
 238 #define numcombs 8
 239 #define numallpasses 4
 240 #define fixedgain 0.015f
 241 #define scalewet 3.0f
 242 #define scaledamp 1.0f
 243 #define scaleroom 0.28f
 244 #define offsetroom 0.7f
 245 #define initialroom 0.5f
 246 #define initialdamp 0.2f
 247 #define initialwet 1
 248 #define initialdry 0
 249 #define initialwidth 1
 250 #define stereospread 23
 251
 252 /*
 253  These values assume 44.1KHz sample rate
 254  they will probably be OK for 48KHz sample rate
 255  but would need scaling for 96KHz (or other) sample rates.
 256  The values were obtained by listening tests.
 257 */
 258 #define combtuningL1 1116
 259 #define combtuningR1 (1116 + stereospread)
 260 #define combtuningL2 1188
 261 #define combtuningR2 (1188 + stereospread)
 262 #define combtuningL3 1277
 263 #define combtuningR3 (1277 + stereospread)
 264 #define combtuningL4 1356
 265 #define combtuningR4 (1356 + stereospread)
 266 #define combtuningL5 1422
 267 #define combtuningR5 (1422 + stereospread)
 268 #define combtuningL6 1491
 269 #define combtuningR6 (1491 + stereospread)
 270 #define combtuningL7 1557
 271 #define combtuningR7 (1557 + stereospread)
 272 #define combtuningL8 1617
 273 #define combtuningR8 (1617 + stereospread)
 274 #define allpasstuningL1 556
 275 #define allpasstuningR1 (556 + stereospread)
 276 #define allpasstuningL2 441
 277 #define allpasstuningR2 (441 + stereospread)
 278 #define allpasstuningL3 341
 279 #define allpasstuningR3 (341 + stereospread)
 280 #define allpasstuningL4 225
 281 #define allpasstuningR4 (225 + stereospread)
 282
 283 struct _fluid_revmodel_t {
 284   fluid_real_t roomsize;
 285   fluid_real_t damp;
 286   fluid_real_t wet, wet1, wet2;
 287   fluid_real_t width;
 288   fluid_real_t gain;
 289   /*
 290    The following are all declared inline
 291    to remove the need for dynamic allocation
 292    with its subsequent error-checking messiness
 293   */
 294   /* Comb filters */
 295   fluid_comb combL[numcombs];
 296   fluid_comb combR[numcombs];
 297   /* Allpass filters */
 298   fluid_allpass allpassL[numallpasses];
 299   fluid_allpass allpassR[numallpasses];
 300 };
 301
 302 static void fluid_revmodel_update(fluid_revmodel_t* rev);
 303 static void fluid_revmodel_init(fluid_revmodel_t* rev);
 304 void fluid_set_revmodel_buffers(fluid_revmodel_t* rev, fluid_real_t sample_rate);
 305
 306 fluid_revmodel_t*
 307 new_fluid_revmodel(fluid_real_t sample_rate)
 308 {
 309   fluid_revmodel_t* rev;
 310   rev = FLUID_NEW(fluid_revmodel_t);
 311   if (rev == NULL) {
 312     return NULL;
 313   }
 314
 315   fluid_set_revmodel_buffers(rev, sample_rate);
 316
 317   /* Set default values */
 318   fluid_allpass_setfeedback(&rev->allpassL[0], 0.5f);
 319   fluid_allpass_setfeedback(&rev->allpassR[0], 0.5f);
 320   fluid_allpass_setfeedback(&rev->allpassL[1], 0.5f);
 321   fluid_allpass_setfeedback(&rev->allpassR[1], 0.5f);
 322   fluid_allpass_setfeedback(&rev->allpassL[2], 0.5f);
 323   fluid_allpass_setfeedback(&rev->allpassR[2], 0.5f);
 324   fluid_allpass_setfeedback(&rev->allpassL[3], 0.5f);
 325   fluid_allpass_setfeedback(&rev->allpassR[3], 0.5f);
 326
 327   rev->gain = fixedgain;
 328   fluid_revmodel_set(rev,FLUID_REVMODEL_SET_ALL,initialroom,initialdamp,initialwidth,initialwet);
 329
 330   return rev;
 331 }
 332
 333 void
 334 delete_fluid_revmodel(fluid_revmodel_t* rev)
 335 {
 336   int i;
 337   for (i = 0; i < numcombs;i++) {
 338     fluid_comb_release(&rev->combL[i]);
 339     fluid_comb_release(&rev->combR[i]);
 340   }
 341   for (i = 0; i < numallpasses; i++) {
 342     fluid_allpass_release(&rev->allpassL[i]);
 343     fluid_allpass_release(&rev->allpassR[i]);
 344   }
 345
 346   FLUID_FREE(rev);
 347 }
 348
 349 void
 350 fluid_set_revmodel_buffers(fluid_revmodel_t* rev, fluid_real_t sample_rate) {
 351
 352   float srfactor = sample_rate/44100.0f;
 353
 354   fluid_comb_setbuffer(&rev->combL[0], combtuningL1*srfactor);
 355   fluid_comb_setbuffer(&rev->combR[0], combtuningR1*srfactor);
 356   fluid_comb_setbuffer(&rev->combL[1], combtuningL2*srfactor);
 357   fluid_comb_setbuffer(&rev->combR[1], combtuningR2*srfactor);
 358   fluid_comb_setbuffer(&rev->combL[2], combtuningL3*srfactor);
 359   fluid_comb_setbuffer(&rev->combR[2], combtuningR3*srfactor);
 360   fluid_comb_setbuffer(&rev->combL[3], combtuningL4*srfactor);
 361   fluid_comb_setbuffer(&rev->combR[3], combtuningR4*srfactor);
 362   fluid_comb_setbuffer(&rev->combL[4], combtuningL5*srfactor);
 363   fluid_comb_setbuffer(&rev->combR[4], combtuningR5*srfactor);
 364   fluid_comb_setbuffer(&rev->combL[5], combtuningL6*srfactor);
 365   fluid_comb_setbuffer(&rev->combR[5], combtuningR6*srfactor);
 366   fluid_comb_setbuffer(&rev->combL[6], combtuningL7*srfactor);
 367   fluid_comb_setbuffer(&rev->combR[6], combtuningR7*srfactor);
 368   fluid_comb_setbuffer(&rev->combL[7], combtuningL8*srfactor);
 369   fluid_comb_setbuffer(&rev->combR[7], combtuningR8*srfactor);
 370   fluid_allpass_setbuffer(&rev->allpassL[0], allpasstuningL1*srfactor);
 371   fluid_allpass_setbuffer(&rev->allpassR[0], allpasstuningR1*srfactor);
 372   fluid_allpass_setbuffer(&rev->allpassL[1], allpasstuningL2*srfactor);
 373   fluid_allpass_setbuffer(&rev->allpassR[1], allpasstuningR2*srfactor);
 374   fluid_allpass_setbuffer(&rev->allpassL[2], allpasstuningL3*srfactor);
 375   fluid_allpass_setbuffer(&rev->allpassR[2], allpasstuningR3*srfactor);
 376   fluid_allpass_setbuffer(&rev->allpassL[3], allpasstuningL4*srfactor);
 377   fluid_allpass_setbuffer(&rev->allpassR[3], allpasstuningR4*srfactor);
 378
 379   /* Clear all buffers */
 380   fluid_revmodel_init(rev);
 381 }
 382
 383
 384 static void
 385 fluid_revmodel_init(fluid_revmodel_t* rev)
 386 {
 387   int i;
 388   for (i = 0; i < numcombs;i++) {
 389     fluid_comb_init(&rev->combL[i]);
 390     fluid_comb_init(&rev->combR[i]);
 391   }
 392   for (i = 0; i < numallpasses; i++) {
 393     fluid_allpass_init(&rev->allpassL[i]);
 394     fluid_allpass_init(&rev->allpassR[i]);
 395   }
 396 }
 397
 398 void
 399 fluid_revmodel_reset(fluid_revmodel_t* rev)
 400 {
 401   fluid_revmodel_init(rev);
 402 }
 403
 404 void
 405 fluid_revmodel_processreplace(fluid_revmodel_t* rev, fluid_real_t *in,
 406                              fluid_real_t *left_out, fluid_real_t *right_out)
 407 {
 408   int i, k = 0;
 409   fluid_real_t outL, outR, input;
 410
 411   for (k = 0; k < FLUID_BUFSIZE; k++) {
 412
 413     outL = outR = 0;
 414
 415     /* The original Freeverb code expects a stereo signal and 'input'
 416      * is set to the sum of the left and right input sample. Since
 417      * this code works on a mono signal, 'input' is set to twice the
 418      * input sample. */
 419     input = (2.0f * in[k] + DC_OFFSET) * rev->gain;
 420
 421     /* Accumulate comb filters in parallel */
 422     for (i = 0; i < numcombs; i++) {
 423       fluid_comb_process(rev->combL[i], input, outL);
 424       fluid_comb_process(rev->combR[i], input, outR);
 425     }
 426     /* Feed through allpasses in series */
 427     for (i = 0; i < numallpasses; i++) {
 428       fluid_allpass_process(rev->allpassL[i], outL);
 429       fluid_allpass_process(rev->allpassR[i], outR);
 430     }
 431
 432     /* Remove the DC offset */
 433     outL -= DC_OFFSET;
 434     outR -= DC_OFFSET;
 435
 436     /* Calculate output REPLACING anything already there */
 437     left_out[k] = outL * rev->wet1 + outR * rev->wet2;
 438     right_out[k] = outR * rev->wet1 + outL * rev->wet2;
 439   }
 440 }
 441
 442 void
 443 fluid_revmodel_processmix(fluid_revmodel_t* rev, fluid_real_t *in,
 444                          fluid_real_t *left_out, fluid_real_t *right_out)
 445 {
 446   int i, k = 0;
 447   fluid_real_t outL, outR, input;
 448
 449   for (k = 0; k < FLUID_BUFSIZE; k++) {
 450
 451     outL = outR = 0;
 452
 453     /* The original Freeverb code expects a stereo signal and 'input'
 454      * is set to the sum of the left and right input sample. Since
 455      * this code works on a mono signal, 'input' is set to twice the
 456      * input sample. */
 457     input = (2.0f * in[k] + DC_OFFSET) * rev->gain;
 458
 459     /* Accumulate comb filters in parallel */
 460     for (i = 0; i < numcombs; i++) {
 461             fluid_comb_process(rev->combL[i], input, outL);
 462             fluid_comb_process(rev->combR[i], input, outR);
 463     }
 464     /* Feed through allpasses in series */
 465     for (i = 0; i < numallpasses; i++) {
 466       fluid_allpass_process(rev->allpassL[i], outL);
 467       fluid_allpass_process(rev->allpassR[i], outR);
 468     }
 469
 470     /* Remove the DC offset */
 471     outL -= DC_OFFSET;
 472     outR -= DC_OFFSET;
 473
 474     /* Calculate output MIXING with anything already there */
 475     left_out[k] += outL * rev->wet1 + outR * rev->wet2;
 476     right_out[k] += outR * rev->wet1 + outL * rev->wet2;
 477   }
 478 }
 479
 480 static void
 481 fluid_revmodel_update(fluid_revmodel_t* rev)
 482 {
 483   /* Recalculate internal values after parameter change */
 484   int i;
 485
 486   rev->wet1 = rev->wet * (rev->width / 2.0f + 0.5f);
 487   rev->wet2 = rev->wet * ((1.0f - rev->width) / 2.0f);
 488
 489   for (i = 0; i < numcombs; i++) {
 490     fluid_comb_setfeedback(&rev->combL[i], rev->roomsize);
 491     fluid_comb_setfeedback(&rev->combR[i], rev->roomsize);
 492   }
 493
 494   for (i = 0; i < numcombs; i++) {
 495     fluid_comb_setdamp(&rev->combL[i], rev->damp);
 496     fluid_comb_setdamp(&rev->combR[i], rev->damp);
 497   }
 498 }
 499
 500 /**
 501  * Set one or more reverb parameters.
 502  * @param rev Reverb instance
 503  * @param set One or more flags from #fluid_revmodel_set_t indicating what
 504  *   parameters to set (#FLUID_REVMODEL_SET_ALL to set all parameters)
 505  * @param roomsize Reverb room size
 506  * @param damping Reverb damping
 507  * @param width Reverb width
 508  * @param level Reverb level
 509  */
 510 void
 511 fluid_revmodel_set(fluid_revmodel_t* rev, int set, float roomsize,
 512                    float damping, float width, float level)
 513 {
 514   if (set & FLUID_REVMODEL_SET_ROOMSIZE)
 515     rev->roomsize = (roomsize * scaleroom) + offsetroom;
 516
 517   if (set & FLUID_REVMODEL_SET_DAMPING)
 518     rev->damp = damping * scaledamp;
 519
 520   if (set & FLUID_REVMODEL_SET_WIDTH)
 521     rev->width = width;
 522
 523   if (set & FLUID_REVMODEL_SET_LEVEL)
 524   {
 525     fluid_clip(level, 0.0f, 1.0f);
 526     rev->wet = level * scalewet;
 527   }
 528
 529   fluid_revmodel_update (rev);
 530 }
 531
 532 void
 533 fluid_revmodel_samplerate_change(fluid_revmodel_t* rev, fluid_real_t sample_rate) {
 534   int i;
 535   for (i = 0; i < numcombs;i++) {
 536     fluid_comb_release(&rev->combL[i]);
 537     fluid_comb_release(&rev->combR[i]);
 538   }
 539   for (i = 0; i < numallpasses; i++) {
 540     fluid_allpass_release(&rev->allpassL[i]);
 541     fluid_allpass_release(&rev->allpassR[i]);
 542   }
 543   fluid_set_revmodel_buffers(rev, sample_rate);
 544 }