ref: 1f3679dfcf9f07dda9a568e62e8004bebd1df4cd
parent: 18bb7971d61b882545c0a0623caf799bd96de464
author: robs <robs>
date: Tue Jun 10 13:54:22 EDT 2008
tune filters
--- a/src/effects_i.c
+++ b/src/effects_i.c
@@ -206,3 +206,14 @@
}
return NULL;
}
+
+double bessel_I_0(double x)
+{+ double term = 1, sum = 1, last_sum, x2 = x / 2;
+ int i = 1;
+ do {+ double y = x2 / i++;
+ last_sum = sum, sum += term *= y * y;
+ } while (sum != last_sum);
+ return sum;
+}
--- a/src/filter.c
+++ b/src/filter.c
@@ -19,8 +19,7 @@
*
*
* REMARKS: (Stan Brooks speaking)
- * This code is heavily based on the resample.c code which was
- * apparently itself a rewrite (by Lance Norskog?) of code originally
+ * This code is a rewrite of filterkit.c from the `resample' package
* by Julius O. Smith, now distributed under the LGPL.
*/
@@ -31,6 +30,7 @@
#define ISCALE 0x10000
#define BUFFSIZE 8192
+#define MAXNWING (80<<7)
/* Private data for Lerp via LCM file */
typedef struct {@@ -45,9 +45,138 @@
double *X, *Y;/* I/O buffers */
} priv_t;
-/* lsx_makeFilter() declared in resample.c */
-extern int
-lsx_makeFilter(double Fp[], long Nwing, double Froll, double Beta, long Num, int Normalize);
+/* lsx_makeFilter is used by resample.c */
+int lsx_makeFilter(double Fp[], long Nwing, double Froll, double Beta, long Num, int Normalize);
+
+/* 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;
+ long 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);
+}
+
+static void LpFilter(double *c, long N, double frq, double Beta, long Num)
+{+ long i;
+
+ /* Calculate filter coeffs: */
+ c[0] = frq;
+ for (i=1; i<N; i++) {+ double x = M_PI*(double)i/(double)(Num);
+ c[i] = sin(x*frq)/x;
+ }
+
+ if (Beta>2) { /* Apply Kaiser window to filter coeffs: */+ double IBeta = 1.0/Izero(Beta);
+ for (i=1; i<N; i++) {+ double x = (double)i / (double)(N);
+ c[i] *= Izero(Beta*sqrt(1.0-x*x)) * IBeta;
+ }
+ } else { /* Apply Nuttall window: */+ for(i = 0; i < N; i++) {+ double x = M_PI*i / N;
+ c[i] *= 0.36335819 + 0.4891775*cos(x) + 0.1365995*cos(2*x) + 0.0106411*cos(3*x);
+ }
+ }
+}
+
+int lsx_makeFilter(double Imp[], long Nwing, double Froll, double Beta,
+ long Num, int Normalize)
+{+ double *ImpR;
+ long Mwing, i;
+
+ if (Nwing > MAXNWING) /* Check for valid parameters */
+ return(-1);
+ if ((Froll<=0) || (Froll>1))
+ return(-2);
+
+ /* it does help accuracy a bit to have the window stop at
+ * a zero-crossing of the sinc function */
+ Mwing = floor((double)Nwing/(Num/Froll))*(Num/Froll) +0.5;
+ if (Mwing==0)
+ return(-4);
+
+ ImpR = lsx_malloc(sizeof(double) * Mwing);
+
+ /* Design a Nuttall or Kaiser windowed Sinc low-pass filter */
+ LpFilter(ImpR, Mwing, Froll, Beta, Num);
+
+ if (Normalize) { /* 'correct' the DC gain of the lowpass filter */+ long Dh;
+ double DCgain;
+ DCgain = 0;
+ Dh = Num; /* Filter sampling period for factors>=1 */
+ for (i=Dh; i<Mwing; i+=Dh)
+ DCgain += ImpR[i];
+ DCgain = 2*DCgain + ImpR[0]; /* DC gain of real coefficients */
+ sox_debug("DCgain err=%.12f",DCgain-1.0);+
+ DCgain = 1.0/DCgain;
+ for (i=0; i<Mwing; i++)
+ Imp[i] = ImpR[i]*DCgain;
+
+ } else {+ for (i=0; i<Mwing; i++)
+ Imp[i] = ImpR[i];
+ }
+ free(ImpR);
+ for (i=Mwing; i<=Nwing; i++) Imp[i] = 0;
+ /* Imp[Mwing] and Imp[-1] needed for quadratic interpolation */
+ Imp[-1] = Imp[1];
+
+ return(Mwing);
+}
static void FiltWin(priv_t * f, long Nx);
--- a/src/polyphas.c
+++ b/src/polyphas.c
@@ -621,7 +621,7 @@
" -width n window width in samples [default 1024]\n"
"\n"
" -cutoff float frequency cutoff for base bandwidth [default 0.95]",
- SOX_EFF_RATE,
+ SOX_EFF_RATE | SOX_EFF_DEPRECATED,
sox_poly_getopts,
sox_poly_start,
sox_poly_flow,
--- a/src/rabbit.c
+++ b/src/rabbit.c
@@ -185,7 +185,7 @@
{ static sox_effect_handler_t handler = {"rabbit", "[-c0|-c1|-c2|-c3|-c4] [rate]",
- SOX_EFF_RATE | SOX_EFF_MCHAN,
+ SOX_EFF_RATE | SOX_EFF_MCHAN | SOX_EFF_DEPRECATED,
getopts, start, flow, drain, stop, NULL, sizeof(priv_t)
};
--- a/src/rate.c
+++ b/src/rate.c
@@ -187,29 +187,19 @@
}
}
-static double bessel_I_0(double x)
-{- double term = 1, sum = 1, last_sum, x2 = x / 2;
- int i = 1;
- do {- double y = x2 / i++;
- last_sum = sum, sum += term *= y * y;
- } while (sum != last_sum);
- return sum;
-}
-
static double * make_lpf(int num_taps, double Fc, double beta, double scale)
{- double * h = malloc(num_taps * sizeof(*h));
+ double * h = malloc(num_taps * sizeof(*h)), sum = 0;
int i, m = num_taps - 1;
assert(Fc >= 0 && Fc <= 1);
- scale /= bessel_I_0(beta);
for (i = 0; i <= m / 2; ++i) {double x = M_PI * (i - .5 * m), y = 2. * i / m - 1;
h[i] = x? sin(Fc * x) / x : Fc;
- h[i] *= bessel_I_0(beta * sqrt(1 - y * y)) * scale;
- h[m - i] = h[i];
+ sum += h[i] *= bessel_I_0(beta * sqrt(1 - y * y));
+ if (m - i != i)
+ sum += h[m - i] = h[i];
}
+ for (i = 0; i < num_taps; ++i) h[i] *= scale / sum;
return h;
}
@@ -380,7 +370,7 @@
static filter_t const filters[] = { {2 * array_length(half_fir_coefs_low) - 1, half_fir_coefs_low, 0,0}, {0, NULL, .986, 110}, {0, NULL, .986, 125},- {0, NULL, .996, 156}, {0, NULL, .999, 156},+ {0, NULL, .996, 165}, {0, NULL, .999, 165},};
filter_t const * f = &filters[quality - Low];
assert((size_t)(quality - Low) < array_length(filters));
@@ -499,7 +489,7 @@
static int create(sox_effect_t * effp, int argc, char **argv)
{priv_t * p = (priv_t *) effp->priv;
- char dummy, * found_at, * bws = "-qlmhvu", * nums = "01";
+ char dummy, * found_at, * bws = "-qlmhvu", * nums = "123";
p->quality = p->coef_interp = -1;
p->shared_ptr = &p->shared;
--- a/src/rate_filters.h
+++ b/src/rate_filters.h
@@ -18,13 +18,14 @@
/* Generated by m4 */
static const sample_t half_fir_coefs_25[] = {- 4.9622034521467345e-001, 3.1227031280072459e-001, 3.5198986166898451e-003,
- -8.9221427418232149e-002, -2.8372758530229671e-003, 3.9139593785072563e-002,
- 1.9671926007707326e-003, -1.7250462144826118e-002, -1.1596954982514773e-003,
- 6.8588796907120527e-003, 5.7031347966698574e-004, -2.3044720814877507e-003,
- -2.2679857326247756e-004, 6.0962502224911946e-004, 6.9133289864740816e-005,
- -1.1323492804087647e-004, -1.4554678571602767e-005, 1.1197369453748703e-005,
- 1.6171529855434563e-006,
+ 4.9866643051942178e-001, 3.1333582318860204e-001, 1.2567743716165585e-003,
+ -9.2035726038137103e-002, -1.0507348255277846e-003, 4.2764945027796687e-002,
+ 7.7661461450703555e-004, -2.0673365323361139e-002, -5.0429677622613805e-004,
+ 9.4223774565849357e-003, 2.8491539998284476e-004, -3.8562347294894628e-003,
+ -1.3803431143314762e-004, 1.3634218103234187e-003, 5.6110366313398705e-005,
+ -3.9872042837864422e-004, -1.8501044952475473e-005, 9.0580351350892191e-005,
+ 4.6764104835321042e-006, -1.4284332593063177e-005, -8.1340436298087893e-007,
+ 1.1833367010222812e-006, 7.3979325233687461e-008,
};
static const sample_t half_fir_coefs_low[] = {4.2759802493108773e-001, 3.0939308096100709e-001, 6.9285325719540158e-002,
@@ -44,7 +45,7 @@
-6.0893901283459912e-006, 1.4040363940567877e-005, 4.9834316581482789e-006,
};
#define FUNCTION half_sample_25
-#define CONVOLVE _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
+#define CONVOLVE _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
#define COEFS half_fir_coefs_25
assert_static(!((array_length(COEFS)- 1) & 1), HALF_FIR_LENGTH_25 );
#include "rate_half_fir.h"
@@ -61,18 +62,25 @@
#include "rate_poly_fir0.h"
#define FUNCTION d100_1
#define COEF_INTERP 1
-#define PHASE_BITS 8
+#define PHASE_BITS 9
#define FIR_LENGTH d100_l
#define CONVOLVE poly_fir_convolve_d100
#include "rate_poly_fir.h"
-#define d100_1_b 8
+#define d100_1_b 9
#define FUNCTION d100_2
#define COEF_INTERP 2
+#define PHASE_BITS 7
+#define FIR_LENGTH d100_l
+#define CONVOLVE poly_fir_convolve_d100
+#include "rate_poly_fir.h"
+#define d100_2_b 7
+#define FUNCTION d100_3
+#define COEF_INTERP 3
#define PHASE_BITS 6
#define FIR_LENGTH d100_l
#define CONVOLVE poly_fir_convolve_d100
#include "rate_poly_fir.h"
-#define d100_2_b 6
+#define d100_3_b 6
#define d120_l 30
#define poly_fir_convolve_d120 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
#define FUNCTION d120_0
@@ -88,11 +96,18 @@
#define d120_1_b 10
#define FUNCTION d120_2
#define COEF_INTERP 2
+#define PHASE_BITS 9
+#define FIR_LENGTH d120_l
+#define CONVOLVE poly_fir_convolve_d120
+#include "rate_poly_fir.h"
+#define d120_2_b 9
+#define FUNCTION d120_3
+#define COEF_INTERP 3
#define PHASE_BITS 7
#define FIR_LENGTH d120_l
#define CONVOLVE poly_fir_convolve_d120
#include "rate_poly_fir.h"
-#define d120_2_b 7
+#define d120_3_b 7
#define d150_l 38
#define poly_fir_convolve_d150 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
#define FUNCTION d150_0
@@ -101,18 +116,25 @@
#include "rate_poly_fir0.h"
#define FUNCTION d150_1
#define COEF_INTERP 1
-#define PHASE_BITS 11
+#define PHASE_BITS 12
#define FIR_LENGTH d150_l
#define CONVOLVE poly_fir_convolve_d150
#include "rate_poly_fir.h"
-#define d150_1_b 11
+#define d150_1_b 12
#define FUNCTION d150_2
#define COEF_INTERP 2
+#define PHASE_BITS 10
+#define FIR_LENGTH d150_l
+#define CONVOLVE poly_fir_convolve_d150
+#include "rate_poly_fir.h"
+#define d150_2_b 10
+#define FUNCTION d150_3
+#define COEF_INTERP 3
#define PHASE_BITS 8
#define FIR_LENGTH d150_l
#define CONVOLVE poly_fir_convolve_d150
#include "rate_poly_fir.h"
-#define d150_2_b 8
+#define d150_3_b 8
#define U100_l 42
#define poly_fir_convolve_U100 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
#define FUNCTION U100_0
@@ -121,18 +143,25 @@
#include "rate_poly_fir0.h"
#define FUNCTION U100_1
#define COEF_INTERP 1
-#define PHASE_BITS 9
+#define PHASE_BITS 10
#define FIR_LENGTH U100_l
#define CONVOLVE poly_fir_convolve_U100
#include "rate_poly_fir.h"
-#define U100_1_b 9
+#define U100_1_b 10
#define FUNCTION U100_2
#define COEF_INTERP 2
+#define PHASE_BITS 8
+#define FIR_LENGTH U100_l
+#define CONVOLVE poly_fir_convolve_U100
+#include "rate_poly_fir.h"
+#define U100_2_b 8
+#define FUNCTION U100_3
+#define COEF_INTERP 3
#define PHASE_BITS 6
#define FIR_LENGTH U100_l
#define CONVOLVE poly_fir_convolve_U100
#include "rate_poly_fir.h"
-#define U100_2_b 6
+#define U100_3_b 6
#define u100_l 10
#define poly_fir_convolve_u100 _ _ _ _ _ _ _ _ _ _
#define FUNCTION u100_0
@@ -141,18 +170,25 @@
#include "rate_poly_fir0.h"
#define FUNCTION u100_1
#define COEF_INTERP 1
-#define PHASE_BITS 7
+#define PHASE_BITS 9
#define FIR_LENGTH u100_l
#define CONVOLVE poly_fir_convolve_u100
#include "rate_poly_fir.h"
-#define u100_1_b 7
+#define u100_1_b 9
#define FUNCTION u100_2
#define COEF_INTERP 2
-#define PHASE_BITS 5
+#define PHASE_BITS 7
#define FIR_LENGTH u100_l
#define CONVOLVE poly_fir_convolve_u100
#include "rate_poly_fir.h"
-#define u100_2_b 5
+#define u100_2_b 7
+#define FUNCTION u100_3
+#define COEF_INTERP 3
+#define PHASE_BITS 6
+#define FIR_LENGTH u100_l
+#define CONVOLVE poly_fir_convolve_u100
+#include "rate_poly_fir.h"
+#define u100_3_b 6
#define u120_l 14
#define poly_fir_convolve_u120 _ _ _ _ _ _ _ _ _ _ _ _ _ _
#define FUNCTION u120_0
@@ -161,18 +197,25 @@
#include "rate_poly_fir0.h"
#define FUNCTION u120_1
#define COEF_INTERP 1
-#define PHASE_BITS 8
+#define PHASE_BITS 10
#define FIR_LENGTH u120_l
#define CONVOLVE poly_fir_convolve_u120
#include "rate_poly_fir.h"
-#define u120_1_b 8
+#define u120_1_b 10
#define FUNCTION u120_2
#define COEF_INTERP 2
+#define PHASE_BITS 8
+#define FIR_LENGTH u120_l
+#define CONVOLVE poly_fir_convolve_u120
+#include "rate_poly_fir.h"
+#define u120_2_b 8
+#define FUNCTION u120_3
+#define COEF_INTERP 3
#define PHASE_BITS 6
#define FIR_LENGTH u120_l
#define CONVOLVE poly_fir_convolve_u120
#include "rate_poly_fir.h"
-#define u120_2_b 6
+#define u120_3_b 6
#define u150_l 20
#define poly_fir_convolve_u150 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
#define FUNCTION u150_0
@@ -181,27 +224,34 @@
#include "rate_poly_fir0.h"
#define FUNCTION u150_1
#define COEF_INTERP 1
-#define PHASE_BITS 10
+#define PHASE_BITS 11
#define FIR_LENGTH u150_l
#define CONVOLVE poly_fir_convolve_u150
#include "rate_poly_fir.h"
-#define u150_1_b 10
+#define u150_1_b 11
#define FUNCTION u150_2
#define COEF_INTERP 2
+#define PHASE_BITS 9
+#define FIR_LENGTH u150_l
+#define CONVOLVE poly_fir_convolve_u150
+#include "rate_poly_fir.h"
+#define u150_2_b 9
+#define FUNCTION u150_3
+#define COEF_INTERP 3
#define PHASE_BITS 7
#define FIR_LENGTH u150_l
#define CONVOLVE poly_fir_convolve_u150
#include "rate_poly_fir.h"
-#define u150_2_b 7
+#define u150_3_b 7
typedef struct {int phase_bits; stage_fn_t fn;} poly_fir1_t;-typedef struct {int num_coefs; double pass, stop, att; poly_fir1_t interp[3];} poly_fir_t;+typedef struct {int num_coefs; double pass, stop, att; poly_fir1_t interp[4];} poly_fir_t; static poly_fir_t const poly_firs[] = {- {d100_l, .75,1.5, 108, {{0, d100_0}, {d100_1_b, d100_1}, {d100_2_b, d100_2}}},- {d120_l, 1, 1.5, 133, {{0, d120_0}, {d120_1_b, d120_1}, {d120_2_b, d120_2}}},- {d150_l, 1, 1.5, 165, {{0, d150_0}, {d150_1_b, d150_1}, {d150_2_b, d150_2}}},- {U100_l, .724, 1, 105, {{0, U100_0}, {U100_1_b, U100_1}, {U100_2_b, U100_2}}},- {u100_l, .3, 1.5, 107, {{0, u100_0}, {u100_1_b, u100_1}, {u100_2_b, u100_2}}},- {u120_l, .5, 1.5, 125, {{0, u120_0}, {u120_1_b, u120_1}, {u120_2_b, u120_2}}},- {u150_l, .5, 1.5, 174, {{0, u150_0}, {u150_1_b, u150_1}, {u150_2_b, u150_2}}},+ {d100_l, .75,1.5, 108, {{0, d100_0}, {d100_1_b, d100_1}, {d100_2_b, d100_2}, {d100_3_b, d100_3}}},+ {d120_l, 1, 1.5, 133, {{0, d120_0}, {d120_1_b, d120_1}, {d120_2_b, d120_2}, {d120_3_b, d120_3}}},+ {d150_l, 1, 1.5, 165, {{0, d150_0}, {d150_1_b, d150_1}, {d150_2_b, d150_2}, {d150_3_b, d150_3}}},+ {U100_l, .724, 1, 105, {{0, U100_0}, {U100_1_b, U100_1}, {U100_2_b, U100_2}, {U100_3_b, U100_3}}},+ {u100_l, .3, 1.5, 107, {{0, u100_0}, {u100_1_b, u100_1}, {u100_2_b, u100_2}, {u100_3_b, u100_3}}},+ {u120_l, .5, 1.5, 125, {{0, u120_0}, {u120_1_b, u120_1}, {u120_2_b, u120_2}, {u120_3_b, u120_3}}},+ {u150_l, .5, 1.5, 174, {{0, u150_0}, {u150_1_b, u150_1}, {u150_2_b, u150_2}, {u150_3_b, u150_3}}},};
--- a/src/sox_i.h
+++ b/src/sox_i.h
@@ -48,6 +48,8 @@
sox_sample_t lsx_gcd(sox_sample_t a, sox_sample_t b);
sox_sample_t lsx_lcm(sox_sample_t a, sox_sample_t b);
+double bessel_I_0(double x);
+
#ifndef HAVE_STRCASECMP
int strcasecmp(const char *s1, const char *s2);
int strncasecmp(char const * s1, char const * s2, size_t n);