ref: 8ac1f3e7c6bc84e110f7597c6530aeb003804fbb
parent: 4e4b8ccde1699fc66113109b0795bdd71a5d284e
author: robs <robs>
date: Wed Jan 28 06:58:10 EST 2009
sinc updates
--- a/src/effects_i_dsp.c
+++ b/src/effects_i_dsp.c
@@ -248,19 +248,21 @@
}
}
-double * lsx_make_lpf(int num_taps, double Fc, double beta, double scale)
+double * lsx_make_lpf(int num_taps, double Fc, double beta, double scale, sox_bool dc_norm)
{- double * h = malloc(num_taps * sizeof(*h)), sum = 0;
int i, m = num_taps - 1;
+ double * h = malloc(num_taps * sizeof(*h)), sum = 0;
+ double mult = scale / lsx_bessel_I_0(beta);
assert(Fc >= 0 && Fc <= 1);
+ lsx_debug("make_lpf(n=%i, Fc=%g beta=%g dc-norm=%i scale=%g)", num_taps, Fc, beta, dc_norm, scale); 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;
- sum += h[i] *= lsx_bessel_I_0(beta * sqrt(1 - y * y));
+ sum += h[i] *= lsx_bessel_I_0(beta * sqrt(1 - y * y)) * mult;
if (m - i != i)
sum += h[m - i] = h[i];
}
- for (i = 0; i < num_taps; ++i) h[i] *= scale / sum;
+ for (i = 0; dc_norm && i < num_taps; ++i) h[i] *= scale / sum;
return h;
}
@@ -298,12 +300,12 @@
if (k)
*num_taps = *num_taps * k - 1;
else k = 1;
- lsx_debug("%i %g %g %g %g", *num_taps, Fp, tr_bw, Fc, beta);- return lsx_make_lpf(*num_taps, (Fc - tr_bw) / k, beta, (double)k);
+ lsx_debug("%g %g %g", Fp, tr_bw, Fc);+ return lsx_make_lpf(*num_taps, (Fc - tr_bw) / k, beta, (double)k, sox_true);
}
#if 0
-static void printf_vector(FILE * f, double h[], int num_points)
+static void printf_vector(double h[], int num_points)
{int i;
printf(
@@ -310,13 +312,22 @@
"# name: b\n"
"# type: matrix\n"
"# rows: %i\n"
- "# columns: 1\n", num_points);
- for (i = 0; i < num_points; ++i)
- printf("%24.16e\n", h[i]);+ "# columns: 1\n", num_points / 2);
+ for (i = 0; i < num_points; i += 2)
+ printf("%24.16e\n", h[i + 1]);exit(0);
}
#endif
+static double safe_log(double x)
+{+ assert(x >= 0);
+ if (x)
+ return log(x);
+ lsx_debug("log(0)");+ return -26;
+}
+
void lsx_fir_to_phase(double * * h, int * len,
int * post_len, double phase0)
{@@ -329,7 +340,7 @@
for (i = 0; i < *len; ++i) work[i] = (*h)[i];
lsx_safe_rdft(work_len, 1, work); /* Cepstral: */
- work[0] = safe_log(work[0]), work[1] = safe_log(work[1]);
+ work[0] = safe_log(fabs(work[0])), work[1] = safe_log(fabs(work[1]));
for (i = 2; i < work_len; i += 2) {work[i] = safe_log(sqrt(sqr(work[i]) + sqr(work[i + 1])));
work[i + 1] = 0;
@@ -344,11 +355,12 @@
/* Some filters require phase unwrapping at this point. Ours give dis-
* continuities only in the stop band, so no need to unwrap in this case. */
+
#if 0
- UNPACK(work, work_len);
+ LSX_UNPACK(work, work_len);
unwrap_phase(work, work_len);
- PACK(work, work_len);
- printf_vector(work, work_len);
+ printf_vector(work, work_len + 2);
+ LSX_PACK(work, work_len);
#endif
for (i = 2; i < work_len; i += 2) /* Interpolate between linear & min phase */
--- a/src/fft4g.h
+++ b/src/fft4g.h
@@ -29,3 +29,7 @@
#define dft_br_len(l) (2 + (1 << (int)(log(l / 2 + .5) / log(2.)) / 2))
#define dft_sc_len(l) (l / 2)
+
+/* Over-allocate h by 2 to use these macros */
+#define LSX_PACK(h, n) h[1] = h[n]
+#define LSX_UNPACK(h, n) h[n] = h[1], h[n + 1] = h[1] = 0;
--- a/src/sinc.c
+++ b/src/sinc.c
@@ -21,9 +21,10 @@
#include <string.h>
typedef struct {- dft_filter_priv_t base;
- double att, phase, freq0, freq1, tbw0, tbw1;
- int num_taps[2];
+ dft_filter_priv_t base;
+ double att, beta, phase, Fc0, Fc1, tbw0, tbw1;
+ int num_taps[2];
+ sox_bool round;
} priv_t;
static int create(sox_effect_t * effp, int argc, char * * argv)
@@ -35,33 +36,36 @@
b->filter_ptr = &b->filter;
p->phase = 50;
- p->att = 120;
+ p->beta = -1;
while (i < 2) {int c = 1;
- while (c && (c = getopt(argc, argv, "+a:p:t:n:MIL")) != -1) switch (c) {- GETOPT_NUMERIC('a', att, 40 , 200)- GETOPT_NUMERIC('p', phase, 0, 100);- GETOPT_NUMERIC('n', num_taps[1], 10, 32767);- case 't': p->tbw1 = lsx_parse_frequency(optarg, &parse_ptr);
- if (p->tbw1 < 1 || *parse_ptr) return lsx_usage(effp);
- break;
+ while (c && (c = getopt(argc, argv, "+ra:b:p:MILt:n:")) != -1) switch (c) {+ char * parse_ptr;
+ case 'r': p->round = sox_true; break;
+ GETOPT_NUMERIC('a', att, 40 , 170)+ GETOPT_NUMERIC('b', beta, 0 , 18)+ GETOPT_NUMERIC('p', phase, 0, 100)case 'M': p->phase = 0; break;
case 'I': p->phase = 25; break;
case 'L': p->phase = 50; break;
+ GETOPT_NUMERIC('n', num_taps[1], 11, 32767)+ case 't': p->tbw1 = lsx_parse_frequency(optarg, &parse_ptr);
+ if (p->tbw1 < 1 || *parse_ptr) return lsx_usage(effp);
+ break;
default: c = 0;
}
- if (!i++ && !(p->tbw1 && p->num_taps[1])) {+ if ((p->att && p->beta >= 0) || (p->tbw1 && p->num_taps[1]))
+ return lsx_usage(effp);
+ if (!i || !p->Fc1)
p->tbw0 = p->tbw1, p->num_taps[0] = p->num_taps[1];
- if (optind < argc) {- parse_ptr = argv[optind++];
- if (*parse_ptr != '-')
- p->freq0 = lsx_parse_frequency(parse_ptr, &parse_ptr);
- if (*parse_ptr == '-')
- p->freq1 = lsx_parse_frequency(parse_ptr+1, &parse_ptr);
- }
+ if (!i++ && optind < argc) {+ if (*(parse_ptr = argv[optind++]) != '-')
+ p->Fc0 = lsx_parse_frequency(parse_ptr, &parse_ptr);
+ if (*parse_ptr == '-')
+ p->Fc1 = lsx_parse_frequency(parse_ptr + 1, &parse_ptr);
}
}
- return optind != argc || p->freq0 < 0 || p->freq1 < 0 || *parse_ptr ?
+ return optind != argc || p->Fc0 < 0 || p->Fc1 < 0 || *parse_ptr ?
lsx_usage(effp) : SOX_SUCCESS;
}
@@ -73,19 +77,23 @@
h[(n - 1) / 2] += 1;
}
-static double * lpf(sox_effect_t * effp, double Fc, double tbw, int * num_taps, double att, double phase)
+static double * lpf(double Fn, double Fc, double tbw, int * num_taps, double att, double * beta, double phase, sox_bool round)
{- double Fn = effp->in_signal.rate * .5;
- if (!Fc) {+ if ((Fc /= Fn) <= 0 || Fc >= 1) {*num_taps = 0;
return NULL;
}
- if (phase != 0 && phase != 50 && phase != 100)
- att *= 34./33; /* negate degradation with intermediate phase */
- *num_taps = *num_taps? *num_taps |= 1 :
- min(32767, lsx_lpf_num_taps(att, (tbw? tbw / Fn : .05) * .5, 0));
- lsx_report("num taps = %i", *num_taps);- return lsx_make_lpf(*num_taps, Fc / Fn, lsx_kaiser_beta(att), 1.);
+ att = att? att : 120;
+ att = phase && phase != 50 && phase != 100? 34./33 * att : att;
+ *beta = *beta < 0? lsx_kaiser_beta(att) : *beta;
+ if (!*num_taps) {+ int n = lsx_lpf_num_taps(att, (tbw? tbw / Fn : .05) * .5, 0);
+ *num_taps = range_limit(n, 11, 32767);
+ if (round)
+ *num_taps = 1 + 2 * (int)((int)((*num_taps / 2) * Fc + .5) / Fc + .5);
+ lsx_report("num taps = %i (from %i)", *num_taps, n);+ }
+ return lsx_make_lpf(*num_taps | 1, Fc, *beta, 1., sox_false);
}
static int start(sox_effect_t * effp)
@@ -94,11 +102,17 @@
dft_filter_filter_t * f = p->base.filter_ptr;
if (!f->num_taps) {+ double Fn = effp->in_signal.rate * .5;
double * h[2];
int i, longer;
- h[0] = lpf(effp, p->freq0, p->tbw0, &p->num_taps[0], p->att, p->phase);
- h[1] = lpf(effp, p->freq1, p->tbw1, &p->num_taps[1], p->att, p->phase);
+
+ if (p->Fc0 >= Fn || p->Fc1 >= Fn) {+ lsx_fail("filter frequency must be less than sample-rate / 2");+ return SOX_EOF;
+ }
+ h[0] = lpf(Fn, p->Fc0, p->tbw0, &p->num_taps[0], p->att, &p->beta, p->phase, p->round);
if (h[0]) invert(h[0], p->num_taps[0]);
+ h[1] = lpf(Fn, p->Fc1, p->tbw1, &p->num_taps[1], p->att, &p->beta, p->phase, p->round);
longer = p->num_taps[1] > p->num_taps[0];
f->num_taps = p->num_taps[longer];
if (h[0] && h[1]) {@@ -105,7 +119,7 @@
for (i = 0; i < p->num_taps[!longer]; ++i)
h[longer][i + (f->num_taps - p->num_taps[!longer])/2] += h[!longer][i];
free(h[!longer]);
- if (p->freq0 < p->freq1)
+ if (p->Fc0 < p->Fc1)
invert(h[longer], f->num_taps);
}
if (p->phase != 50)
@@ -115,9 +129,9 @@
100 - 100. * f->post_peak / (f->num_taps - 1));
if (effp->global_info->plot != sox_plot_off) {char title[100];
- sprintf(title, "SoX effect: sinc filter freq=%g-%g", p->freq0, p->freq1);
+ sprintf(title, "SoX effect: sinc filter freq=%g-%g", p->Fc0, p->Fc1? p->Fc1 : Fn);
lsx_plot_fir(h[longer], f->num_taps, effp->in_signal.rate,
- effp->global_info->plot, title, -p->att - 20, 10.);
+ effp->global_info->plot, title, -p->beta * 10 - 25, 5.);
return SOX_EOF;
}
f->dft_length = lsx_set_dft_length(f->num_taps);
@@ -135,7 +149,7 @@
static sox_effect_handler_t handler;
handler = *sox_dft_filter_effect_fn();
handler.name = "sinc";
- handler.usage = "[-a att] [-p phase|-M|-I|-L] [-t tbw|-n taps] freq1[-freq2 [-t ...|-n ...]]";
+ handler.usage = "[-a att] [-p phase|-M|-I|-L] [-t tbw|-n taps] [freqHP][-freqLP [-t ...|-n ...]]";
handler.getopts = create;
handler.start = start;
handler.priv_size = sizeof(priv_t);
--- a/src/sox_i.h
+++ b/src/sox_i.h
@@ -107,7 +107,7 @@
void lsx_apply_blackman_nutall(double h[], const int num_points);
double lsx_kaiser_beta(double att);
void lsx_apply_kaiser(double h[], const int num_points, double beta);
-double * lsx_make_lpf(int num_taps, double Fc, double beta, double scale);
+double * lsx_make_lpf(int num_taps, double Fc, double beta, double scale, sox_bool dc_norm);
int lsx_lpf_num_taps(double att, double tr_bw, int k);
double * lsx_design_lpf(
double Fp, /* End of pass-band; ~= 0.01dB point */
--- a/src/util.h
+++ b/src/util.h
@@ -114,7 +114,6 @@
#define sqr(a) ((a) * (a))
#define sign(x) ((x) < 0? -1 : 1)
-#define safe_log(x) ((x) <= 0? -99 : log(x)) /* exp --> 0 (for single prec.) */
/* Numerical Recipes in C, p. 284 */
#define ranqd1(x) ((x) = 1664525L * (x) + 1013904223L) /* int32_t x */