libs/fluidsynth/src/fluid_iir_filter.c

   1 /* FluidSynth - A Software Synthesizer
   2  *
   3  * Copyright (C) 2003  Peter Hanappe and others.
   4  *
   5  * This library is free software; you can redistribute it and/or
   6  * modify it under the terms of the GNU Lesser General Public License
   7  * as published by the Free Software Foundation; either version 2.1 of
   8  * the License, or (at your option) any later version.
   9  *
  10  * This library is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * Lesser General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU Lesser General Public
  16  * License along with this library; if not, write to the Free
  17  * Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  18  * 02110-1301, USA
  19  */
  20
  21 #include "fluid_iir_filter.h"
  22 #include "fluid_sys.h"
  23 #include "fluid_conv.h"
  24
  25 /**
  26  * Applies a low- or high-pass filter with variable cutoff frequency and quality factor
  27  * for a given biquad transfer function:
  28  *          b0 + b1*z^-1 + b2*z^-2
  29  *  H(z) = ------------------------
  30  *          a0 + a1*z^-1 + a2*z^-2
  31  *
  32  * Also modifies filter state accordingly.
  33  * @param iir_filter Filter parameter
  34  * @param dsp_buf Pointer to the synthesized audio data
  35  * @param count Count of samples in dsp_buf
  36  */
  37 /*
  38  * Variable description:
  39  * - dsp_a1, dsp_a2: Filter coefficients for the the previously filtered output signal
  40  * - dsp_b0, dsp_b1, dsp_b2: Filter coefficients for input signal
  41  * - coefficients normalized to a0
  42  *
  43  * A couple of variables are used internally, their results are discarded:
  44  * - dsp_i: Index through the output buffer
  45  * - dsp_centernode: delay line for the IIR filter
  46  * - dsp_hist1: same
  47  * - dsp_hist2: same
  48  */
  49 void
  50 fluid_iir_filter_apply(fluid_iir_filter_t *iir_filter,
  51                        fluid_real_t *dsp_buf, int count)
  52 {
  53     if(iir_filter->type == FLUID_IIR_DISABLED || iir_filter->q_lin == 0)
  54     {
  55         return;
  56     }
  57     else
  58     {
  59         /* IIR filter sample history */
  60         fluid_real_t dsp_hist1 = iir_filter->hist1;
  61         fluid_real_t dsp_hist2 = iir_filter->hist2;
  62
  63         /* IIR filter coefficients */
  64         fluid_real_t dsp_a1 = iir_filter->a1;
  65         fluid_real_t dsp_a2 = iir_filter->a2;
  66         fluid_real_t dsp_b02 = iir_filter->b02;
  67         fluid_real_t dsp_b1 = iir_filter->b1;
  68         int dsp_filter_coeff_incr_count = iir_filter->filter_coeff_incr_count;
  69
  70         fluid_real_t dsp_centernode;
  71         int dsp_i;
  72
  73         /* filter (implement the voice filter according to SoundFont standard) */
  74
  75         /* Check for denormal number (too close to zero). */
  76         if(fabs(dsp_hist1) < 1e-20)
  77         {
  78             dsp_hist1 = 0.0f;    /* FIXME JMG - Is this even needed? */
  79         }
  80
  81         /* Two versions of the filter loop. One, while the filter is
  82         * changing towards its new setting. The other, if the filter
  83         * doesn't change.
  84         */
  85
  86         if(dsp_filter_coeff_incr_count > 0)
  87         {
  88             fluid_real_t dsp_a1_incr = iir_filter->a1_incr;
  89             fluid_real_t dsp_a2_incr = iir_filter->a2_incr;
  90             fluid_real_t dsp_b02_incr = iir_filter->b02_incr;
  91             fluid_real_t dsp_b1_incr = iir_filter->b1_incr;
  92
  93
  94             /* Increment is added to each filter coefficient filter_coeff_incr_count times. */
  95             for(dsp_i = 0; dsp_i < count; dsp_i++)
  96             {
  97                 /* The filter is implemented in Direct-II form. */
  98                 dsp_centernode = dsp_buf[dsp_i] - dsp_a1 * dsp_hist1 - dsp_a2 * dsp_hist2;
  99                 dsp_buf[dsp_i] = dsp_b02 * (dsp_centernode + dsp_hist2) + dsp_b1 * dsp_hist1;
 100                 dsp_hist2 = dsp_hist1;
 101                 dsp_hist1 = dsp_centernode;
 102
 103                 if(dsp_filter_coeff_incr_count-- > 0)
 104                 {
 105                     fluid_real_t old_b02 = dsp_b02;
 106                     dsp_a1 += dsp_a1_incr;
 107                     dsp_a2 += dsp_a2_incr;
 108                     dsp_b02 += dsp_b02_incr;
 109                     dsp_b1 += dsp_b1_incr;
 110
 111                     /* Compensate history to avoid the filter going havoc with large frequency changes */
 112                     if(iir_filter->compensate_incr && fabs(dsp_b02) > 0.001)
 113                     {
 114                         fluid_real_t compensate = old_b02 / dsp_b02;
 115                         dsp_hist1 *= compensate;
 116                         dsp_hist2 *= compensate;
 117                     }
 118                 }
 119             } /* for dsp_i */
 120         }
 121         else /* The filter parameters are constant.  This is duplicated to save time. */
 122         {
 123             for(dsp_i = 0; dsp_i < count; dsp_i++)
 124             {
 125                 /* The filter is implemented in Direct-II form. */
 126                 dsp_centernode = dsp_buf[dsp_i] - dsp_a1 * dsp_hist1 - dsp_a2 * dsp_hist2;
 127                 dsp_buf[dsp_i] = dsp_b02 * (dsp_centernode + dsp_hist2) + dsp_b1 * dsp_hist1;
 128                 dsp_hist2 = dsp_hist1;
 129                 dsp_hist1 = dsp_centernode;
 130             }
 131         }
 132
 133         iir_filter->hist1 = dsp_hist1;
 134         iir_filter->hist2 = dsp_hist2;
 135         iir_filter->a1 = dsp_a1;
 136         iir_filter->a2 = dsp_a2;
 137         iir_filter->b02 = dsp_b02;
 138         iir_filter->b1 = dsp_b1;
 139         iir_filter->filter_coeff_incr_count = dsp_filter_coeff_incr_count;
 140
 141         fluid_check_fpe("voice_filter");
 142     }
 143 }
 144
 145
 146 DECLARE_FLUID_RVOICE_FUNCTION(fluid_iir_filter_init)
 147 {
 148     fluid_iir_filter_t *iir_filter = obj;
 149     enum fluid_iir_filter_type type = param[0].i;
 150     enum fluid_iir_filter_flags flags = param[1].i;
 151
 152     iir_filter->type = type;
 153     iir_filter->flags = flags;
 154
 155     if(type != FLUID_IIR_DISABLED)
 156     {
 157         fluid_iir_filter_reset(iir_filter);
 158     }
 159 }
 160
 161 void
 162 fluid_iir_filter_reset(fluid_iir_filter_t *iir_filter)
 163 {
 164     iir_filter->hist1 = 0;
 165     iir_filter->hist2 = 0;
 166     iir_filter->last_fres = -1.;
 167     iir_filter->q_lin = 0;
 168     iir_filter->filter_startup = 1;
 169 }
 170
 171 DECLARE_FLUID_RVOICE_FUNCTION(fluid_iir_filter_set_fres)
 172 {
 173     fluid_iir_filter_t *iir_filter = obj;
 174     fluid_real_t fres = param[0].real;
 175
 176     iir_filter->fres = fres;
 177     iir_filter->last_fres = -1.;
 178 }
 179
 180 static fluid_real_t fluid_iir_filter_q_from_dB(fluid_real_t q_dB)
 181 {
 182     /* The generator contains 'centibels' (1/10 dB) => divide by 10 to
 183      * obtain dB */
 184     q_dB /= 10.0f;
 185
 186     /* Range: SF2.01 section 8.1.3 # 8 (convert from cB to dB => /10) */
 187     fluid_clip(q_dB, 0.0f, 96.0f);
 188
 189     /* Short version: Modify the Q definition in a way, that a Q of 0
 190      * dB leads to no resonance hump in the freq. response.
 191      *
 192      * Long version: From SF2.01, page 39, item 9 (initialFilterQ):
 193      * "The gain at the cutoff frequency may be less than zero when
 194      * zero is specified".  Assume q_dB=0 / q_lin=1: If we would leave
 195      * q as it is, then this results in a 3 dB hump slightly below
 196      * fc. At fc, the gain is exactly the DC gain (0 dB).  What is
 197      * (probably) meant here is that the filter does not show a
 198      * resonance hump for q_dB=0. In this case, the corresponding
 199      * q_lin is 1/sqrt(2)=0.707.  The filter should have 3 dB of
 200      * attenuation at fc now.  In this case Q_dB is the height of the
 201      * resonance peak not over the DC gain, but over the frequency
 202      * response of a non-resonant filter.  This idea is implemented as
 203      * follows: */
 204     q_dB -= 3.01f;
 205
 206     /* The 'sound font' Q is defined in dB. The filter needs a linear
 207        q. Convert. */
 208     return pow(10.0f, q_dB / 20.0f);
 209 }
 210
 211 DECLARE_FLUID_RVOICE_FUNCTION(fluid_iir_filter_set_q)
 212 {
 213     fluid_iir_filter_t *iir_filter = obj;
 214     fluid_real_t q = param[0].real;
 215     int flags = iir_filter->flags;
 216
 217     if(flags & FLUID_IIR_Q_ZERO_OFF && q <= 0.0)
 218     {
 219         q = 0;
 220     }
 221     else if(flags & FLUID_IIR_Q_LINEAR)
 222     {
 223         /* q is linear (only for user-defined filter)
 224          * increase to avoid Q being somewhere between zero and one,
 225          * which results in some strange amplified lowpass signal
 226          */
 227         q++;
 228     }
 229     else
 230     {
 231         q = fluid_iir_filter_q_from_dB(q);
 232     }
 233
 234     iir_filter->q_lin = q;
 235     iir_filter->filter_gain = 1.0;
 236
 237     if(!(flags & FLUID_IIR_NO_GAIN_AMP))
 238     {
 239         /* SF 2.01 page 59:
 240          *
 241          *  The SoundFont specs ask for a gain reduction equal to half the
 242          *  height of the resonance peak (Q).  For example, for a 10 dB
 243          *  resonance peak, the gain is reduced by 5 dB.  This is done by
 244          *  multiplying the total gain with sqrt(1/Q).  `Sqrt' divides dB
 245          *  by 2 (100 lin = 40 dB, 10 lin = 20 dB, 3.16 lin = 10 dB etc)
 246          *  The gain is later factored into the 'b' coefficients
 247          *  (numerator of the filter equation).  This gain factor depends
 248          *  only on Q, so this is the right place to calculate it.
 249          */
 250         iir_filter->filter_gain /= sqrt(q);
 251     }
 252
 253     /* The synthesis loop will have to recalculate the filter coefficients. */
 254     iir_filter->last_fres = -1.;
 255 }
 256
 257 static FLUID_INLINE void
 258 fluid_iir_filter_calculate_coefficients(fluid_iir_filter_t *iir_filter,
 259                                         int transition_samples,
 260                                         fluid_real_t output_rate)
 261 {
 262     /* FLUID_IIR_Q_LINEAR may switch the filter off by setting Q==0 */
 263     if(iir_filter->q_lin == 0)
 264     {
 265         return;
 266     }
 267     else
 268     {
 269         /*
 270          * Those equations from Robert Bristow-Johnson's `Cookbook
 271          * formulae for audio EQ biquad filter coefficients', obtained
 272          * from Harmony-central.com / Computer / Programming. They are
 273          * the result of the bilinear transform on an analogue filter
 274          * prototype. To quote, `BLT frequency warping has been taken
 275          * into account for both significant frequency relocation and for
 276          * bandwidth readjustment'. */
 277
 278         fluid_real_t omega = (fluid_real_t)(2.0 * M_PI) *
 279                                             (iir_filter->last_fres / output_rate);
 280         fluid_real_t sin_coeff = (fluid_real_t) sin(omega);
 281         fluid_real_t cos_coeff = (fluid_real_t) cos(omega);
 282         fluid_real_t alpha_coeff = sin_coeff / (2.0f * iir_filter->q_lin);
 283         fluid_real_t a0_inv = 1.0f / (1.0f + alpha_coeff);
 284
 285         /* Calculate the filter coefficients. All coefficients are
 286          * normalized by a0. Think of `a1' as `a1/a0'.
 287          *
 288          * Here a couple of multiplications are saved by reusing common expressions.
 289          * The original equations should be:
 290          *  iir_filter->b0=(1.-cos_coeff)*a0_inv*0.5*iir_filter->filter_gain;
 291          *  iir_filter->b1=(1.-cos_coeff)*a0_inv*iir_filter->filter_gain;
 292          *  iir_filter->b2=(1.-cos_coeff)*a0_inv*0.5*iir_filter->filter_gain; */
 293
 294         /* "a" coeffs are same for all 3 available filter types */
 295         fluid_real_t a1_temp = -2.0f * cos_coeff * a0_inv;
 296         fluid_real_t a2_temp = (1.0f - alpha_coeff) * a0_inv;
 297
 298         fluid_real_t b02_temp, b1_temp;
 299
 300         switch(iir_filter->type)
 301         {
 302         case FLUID_IIR_HIGHPASS:
 303             b1_temp = (1.0f + cos_coeff) * a0_inv * iir_filter->filter_gain;
 304
 305             /* both b0 -and- b2 */
 306             b02_temp = b1_temp * 0.5f;
 307
 308             b1_temp *= -1.0f;
 309             break;
 310
 311         case FLUID_IIR_LOWPASS:
 312             b1_temp = (1.0f - cos_coeff) * a0_inv * iir_filter->filter_gain;
 313
 314             /* both b0 -and- b2 */
 315             b02_temp = b1_temp * 0.5f;
 316             break;
 317
 318         default:
 319             /* filter disabled, should never get here */
 320             return;
 321         }
 322
 323         iir_filter->compensate_incr = 0;
 324
 325         if(iir_filter->filter_startup || (transition_samples == 0))
 326         {
 327             /* The filter is calculated, because the voice was started up.
 328              * In this case set the filter coefficients without delay.
 329              */
 330             iir_filter->a1 = a1_temp;
 331             iir_filter->a2 = a2_temp;
 332             iir_filter->b02 = b02_temp;
 333             iir_filter->b1 = b1_temp;
 334             iir_filter->filter_coeff_incr_count = 0;
 335             iir_filter->filter_startup = 0;
 336 //       printf("Setting initial filter coefficients.\n");
 337         }
 338         else
 339         {
 340
 341             /* The filter frequency is changed.  Calculate an increment
 342              * factor, so that the new setting is reached after one buffer
 343              * length. x_incr is added to the current value FLUID_BUFSIZE
 344              * times. The length is arbitrarily chosen. Longer than one
 345              * buffer will sacrifice some performance, though.  Note: If
 346              * the filter is still too 'grainy', then increase this number
 347              * at will.
 348              */
 349
 350             iir_filter->a1_incr = (a1_temp - iir_filter->a1) / transition_samples;
 351             iir_filter->a2_incr = (a2_temp - iir_filter->a2) / transition_samples;
 352             iir_filter->b02_incr = (b02_temp - iir_filter->b02) / transition_samples;
 353             iir_filter->b1_incr = (b1_temp - iir_filter->b1) / transition_samples;
 354
 355             if(fabs(iir_filter->b02) > 0.0001)
 356             {
 357                 fluid_real_t quota = b02_temp / iir_filter->b02;
 358                 iir_filter->compensate_incr = quota < 0.5 || quota > 2;
 359             }
 360
 361             /* Have to add the increments filter_coeff_incr_count times. */
 362             iir_filter->filter_coeff_incr_count = transition_samples;
 363         }
 364
 365         fluid_check_fpe("voice_write filter calculation");
 366     }
 367 }
 368
 369
 370 void fluid_iir_filter_calc(fluid_iir_filter_t *iir_filter,
 371                            fluid_real_t output_rate,
 372                            fluid_real_t fres_mod)
 373 {
 374     fluid_real_t fres;
 375
 376     /* calculate the frequency of the resonant filter in Hz */
 377     fres = fluid_ct2hz(iir_filter->fres + fres_mod);
 378
 379     /* FIXME - Still potential for a click during turn on, can we interpolate
 380        between 20khz cutoff and 0 Q? */
 381
 382     /* I removed the optimization of turning the filter off when the
 383      * resonance frequence is above the maximum frequency. Instead, the
 384      * filter frequency is set to a maximum of 0.45 times the sampling
 385      * rate. For a 44100 kHz sampling rate, this amounts to 19845
 386      * Hz. The reason is that there were problems with anti-aliasing when the
 387      * synthesizer was run at lower sampling rates. Thanks to Stephan
 388      * Tassart for pointing me to this bug. By turning the filter on and
 389      * clipping the maximum filter frequency at 0.45*srate, the filter
 390      * is used as an anti-aliasing filter. */
 391
 392     if(fres > 0.45f * output_rate)
 393     {
 394         fres = 0.45f * output_rate;
 395     }
 396     else if(fres < 5)
 397     {
 398         fres = 5;
 399     }
 400
 401     /* if filter enabled and there is a significant frequency change.. */
 402     if(iir_filter->type != FLUID_IIR_DISABLED && fabs(fres - iir_filter->last_fres) > 0.01)
 403     {
 404         /* The filter coefficients have to be recalculated (filter
 405          * parameters have changed). Recalculation for various reasons is
 406          * forced by setting last_fres to -1.  The flag filter_startup
 407          * indicates, that the DSP loop runs for the first time, in this
 408          * case, the filter is set directly, instead of smoothly fading
 409          * between old and new settings. */
 410         iir_filter->last_fres = fres;
 411         fluid_iir_filter_calculate_coefficients(iir_filter, FLUID_BUFSIZE,
 412                                                 output_rate);
 413     }
 414
 415
 416     fluid_check_fpe("voice_write DSP coefficients");
 417
 418 }
 419