blob: bc92285f2c0511446d411f746d397b58b20eaf02 [file] [log] [blame]
 /********************************************************************** resamplesubs.c Real-time library interface by Dominic Mazzoni Based on resample-1.7: http://www-ccrma.stanford.edu/~jos/resample/ License: LGPL - see the file LICENSE.txt for more information This file provides Kaiser-windowed low-pass filter support, including a function to create the filter coefficients, and two functions to apply the filter at a particular point. **********************************************************************/ /* Definitions */ #include "resample_defs.h" #include "filterkit.h" #include #include #include #include /* LpFilter() * * reference: "Digital Filters, 2nd edition" * R.W. Hamming, pp. 178-179 * * Izero() computes the 0th order modified bessel function of the first kind. * (Needed to compute Kaiser window). * * LpFilter() computes the coeffs of a Kaiser-windowed low pass filter with * the following characteristics: * * c[] = array in which to store computed coeffs * frq = roll-off frequency of filter * N = Half the window length in number of coeffs * Beta = parameter of Kaiser window * Num = number of coeffs before 1/frq * * Beta trades the rejection of the lowpass filter against the transition * width from passband to stopband. Larger Beta means a slower * transition and greater stopband rejection. See Rabiner and Gold * (Theory and Application of DSP) under Kaiser windows for more about * Beta. The following table from Rabiner and Gold gives some feel * for the effect of Beta: * * All ripples in dB, width of transition band = D*N where N = window length * * BETA D PB RIP SB RIP * 2.120 1.50 +-0.27 -30 * 3.384 2.23 0.0864 -40 * 4.538 2.93 0.0274 -50 * 5.658 3.62 0.00868 -60 * 6.764 4.32 0.00275 -70 * 7.865 5.0 0.000868 -80 * 8.960 5.7 0.000275 -90 * 10.056 6.4 0.000087 -100 */ #define IzeroEPSILON 1E-21 /* Max error acceptable in Izero */ static double Izero(double x) { double sum, u, halfx, temp; int n; sum = u = n = 1; halfx = x/2.0; do { temp = halfx/(double)n; n += 1; temp *= temp; u *= temp; sum += u; } while (u >= IzeroEPSILON*sum); return(sum); } void lrsLpFilter(double c[], int N, double frq, double Beta, int Num) { double IBeta, temp, temp1, inm1; int i; /* Calculate ideal lowpass filter impulse response coefficients: */ c[0] = 2.0*frq; for (i=1; i