ref: a15706d1e3fa2e7e4c8588434485cbd9526707ca
parent: 8025a99ac69ba3354f478316e5c209772eea3ee0
author: robs <robs>
date: Sat Dec 13 12:43:54 EST 2008
moving some things around
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -36,7 +36,7 @@
compandt filter overdrive reverse tremolo
contrast flanger pad silence trim
dcshift input pan skeleff vol
- delay loudness phaser speed
+ delay loudness phaser speed effects_i_dsp
dither mcompand polyphas splice
)
set(formats_srcs
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -273,7 +273,7 @@
polyphas.c rabbit.c rate.c rate_filters.h rate_half_fir.h \
rate_poly_fir0.h rate_poly_fir.h remix.c repeat.c resample.c reverb.c \
reverse.c silence.c skeleff.c speed.c splice.c stat.c swap.c stretch.c \
- synth.c tempo.c tremolo.c trim.c vol.c overdrive.c
+ synth.c tempo.c tremolo.c trim.c vol.c overdrive.c effects_i_dsp.c
if HAVE_PNG
libsox_la_SOURCES += spectrogram.c
endif
--- a/src/effects_i.c
+++ b/src/effects_i.c
@@ -18,12 +18,7 @@
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
-#ifdef NDEBUG /* Enable assert always. */
-#undef NDEBUG /* Must undef above assert.h or other that might include it. */
-#endif
-
#include "sox_i.h"
-#include <assert.h>
#include <string.h>
#undef lsx_fail
@@ -68,56 +63,6 @@
return a * (b / lsx_gcd(a, b));
}
-/* Numerical Recipes cubic spline */
-
-void lsx_prepare_spline3(double const * x, double const * y, int n,
- double start_1d, double end_1d, double * y_2d)
-{- double p, qn, sig, un, * u = lsx_malloc((n - 1) * sizeof(*u));
- int i;
-
- if (start_1d == HUGE_VAL)
- y_2d[0] = u[0] = 0; /* Start with natural spline or */
- else { /* set the start first derivative */- y_2d[0] = -.5;
- u[0] = (3 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - start_1d);
- }
-
- for (i = 1; i < n - 1; ++i) {- sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
- p = sig * y_2d[i - 1] + 2;
- y_2d[i] = (sig - 1) / p;
- u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) -
- (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
- u[i] = (6 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
- }
- if (end_1d == HUGE_VAL)
- qn = un = 0; /* End with natural spline or */
- else { /* set the end first derivative */- qn = .5;
- un = 3 / (x[n - 1] - x[n - 2]) * (end_1d - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
- }
- y_2d[n - 1] = (un - qn * u[n - 2]) / (qn * y_2d[n - 2] + 1);
- for (i = n - 2; i >= 0; --i)
- y_2d[i] = y_2d[i] * y_2d[i + 1] + u[i];
- free(u);
-}
-
-double lsx_spline3(double const * x, double const * y, double const * y_2d,
- int n, double x1)
-{- int t, i[2] = {0, 0};- double d, a, b;
-
- for (i[1] = n - 1; i[1] - i[0] > 1; t = (i[1] + i[0]) >> 1, i[x[t] > x1] = t);
- d = x[i[1]] - x[i[0]];
- assert(d != 0);
- a = (x[i[1]] - x1) / d;
- b = (x1 - x[i[0]]) / d;
- return a * y[i[0]] + b * y[i[1]] +
- ((a * a * a - a) * y_2d[i[0]] + (b * b * b - b) * y_2d[i[1]]) * d * d / 6;
-}
-
lsx_enum_item const lsx_wave_enum[] = {LSX_ENUM_ITEM(SOX_WAVE_,SINE)
LSX_ENUM_ITEM(SOX_WAVE_,TRIANGLE)
@@ -270,6 +215,8 @@
#if 0
+#include <assert.h>
+
#define TEST(st, samp, len) \
str = st; \
next = lsx_parsesamples(10000, str, &samples, 't'); \
@@ -374,178 +321,6 @@
}
}
return result < 0 ? -1 : result;
-}
-
-double lsx_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;
-}
-
-int lsx_set_dft_length(int num_taps) /* Set to 4 x nearest power of 2 */
-{- int result, n = num_taps;
- for (result = 8; n > 2; result <<= 1, n >>= 1);
- result = range_limit(result, 4096, 131072);
- assert(num_taps * 2 < result);
- return result;
-}
-
-#include "fft4g.h"
-int * lsx_fft_br;
-double * lsx_fft_sc;
-static int fft_len;
-#ifdef HAVE_OPENMP
-static omp_lock_t fft_cache_lock;
-#endif
-
-static void init_fft_cache(void)
-{- omp_init_lock(&fft_cache_lock);
-}
-
-static void clear_fft_cache(void)
-{- omp_destroy_lock(&fft_cache_lock);
- free(lsx_fft_br);
- free(lsx_fft_sc);
- lsx_fft_sc = NULL;
- lsx_fft_br = NULL;
- fft_len = 0;
-}
-
-static void update_fft_cache(int len)
-{- omp_set_lock(&fft_cache_lock);
- if (len > fft_len) {- int old_n = fft_len;
- fft_len = len;
- lsx_fft_br = lsx_realloc(lsx_fft_br, dft_br_len(fft_len) * sizeof(*lsx_fft_br));
- lsx_fft_sc = lsx_realloc(lsx_fft_sc, dft_sc_len(fft_len) * sizeof(*lsx_fft_sc));
- if (!old_n)
- lsx_fft_br[0] = 0;
- }
-}
-
-static sox_bool is_power_of_2(int x)
-{- return !(x < 2 || (x & (x - 1)));
-}
-
-void lsx_safe_rdft(int len, int type, double * d)
-{- assert(is_power_of_2(len));
- update_fft_cache(len);
- lsx_rdft(len, type, d, lsx_fft_br, lsx_fft_sc);
- omp_unset_lock(&fft_cache_lock);
-}
-
-void lsx_safe_cdft(int len, int type, double * d)
-{- assert(is_power_of_2(len));
- update_fft_cache(len);
- lsx_cdft(len, type, d, lsx_fft_br, lsx_fft_sc);
- omp_unset_lock(&fft_cache_lock);
-}
-
-void lsx_power_spectrum(int n, double const * in, double * out)
-{- int i;
- double * work = lsx_memdup(in, n * sizeof(*work));
- lsx_safe_rdft(n, 1, work);
- out[0] = sqr(work[0]);
- for (i = 2; i < n; i += 2)
- out[i >> 1] = sqr(work[i]) + sqr(work[i + 1]);
- out[i >> 1] = sqr(work[1]);
- free(work);
-}
-
-void lsx_power_spectrum_f(int n, float const * in, float * out)
-{- int i;
- double * work = lsx_malloc(n * sizeof(*work));
- for (i = 0; i< n; ++i) work[i] = in[i];
- lsx_safe_rdft(n, 1, work);
- out[0] = sqr(work[0]);
- for (i = 2; i < n; i += 2)
- out[i >> 1] = sqr(work[i]) + sqr(work[i + 1]);
- out[i >> 1] = sqr(work[1]);
- free(work);
-}
-
-void lsx_apply_hann_f(float h[], const int num_points)
-{- int i, m = num_points - 1;
- for (i = 0; i < num_points; ++i) {- double x = 2 * M_PI * i / m;
- h[i] *= .5 - .5 * cos(x);
- }
-}
-
-void lsx_apply_hann(double h[], const int num_points)
-{- int i, m = num_points - 1;
- for (i = 0; i < num_points; ++i) {- double x = 2 * M_PI * i / m;
- h[i] *= .5 - .5 * cos(x);
- }
-}
-
-void lsx_apply_hamming(double h[], const int num_points)
-{- int i, m = num_points - 1;
- for (i = 0; i < num_points; ++i) {- double x = 2 * M_PI * i / m;
- h[i] *= .53836 - .46164 * cos(x);
- }
-}
-
-void lsx_apply_bartlett(double h[], const int num_points)
-{- int i, m = num_points - 1;
- for (i = 0; i < num_points; ++i) {- h[i] *= 2. / m * (m / 2. - fabs(i - m / 2.));
- }
-}
-
-void lsx_apply_blackman(double h[], const int num_points, double alpha /*.16*/)
-{- int i, m = num_points - 1;
- for (i = 0; i < num_points; ++i) {- double x = 2 * M_PI * i / m;
- h[i] *= (1 - alpha) *.5 - .5 * cos(x) + alpha * .5 * cos(2 * x);
- }
-}
-
-void lsx_apply_blackman_nutall(double h[], const int num_points)
-{- int i, m = num_points - 1;
- for (i = 0; i < num_points; ++i) {- double x = 2 * M_PI * i / m;
- h[i] *= .3635819 - .4891775 * cos(x) + .1365995 * cos(2 * x) - .0106411 * cos(3 * x);
- }
-}
-
-double lsx_kaiser_beta(double att)
-{- if (att > 100 ) return .1117 * att - 1.11;
- if (att > 50 ) return .1102 * (att - 8.7);
- if (att > 20.96) return .58417 * pow(att -20.96, .4) + .07886 * (att - 20.96);
- return 0;
-}
-
-void lsx_apply_kaiser(double h[], const int num_points, double beta)
-{- int i, m = num_points - 1;
- for (i = 0; i <= m; ++i) {- double x = 2. * i / m - 1;
- h[i] *= lsx_bessel_I_0(beta * sqrt(1 - x * x)) / lsx_bessel_I_0(beta);
- }
}
int lsx_effects_init(void)
--- /dev/null
+++ b/src/effects_i_dsp.c
@@ -1,0 +1,359 @@
+/* libSoX internal DSP functions.
+ * All public functions & data are prefixed with lsx_ .
+ *
+ * Copyright (c) 2008 robs@users.sourceforge.net
+ *
+ * This library is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation; either version 2.1 of the License, or (at
+ * your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this library; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#ifdef NDEBUG /* Enable assert always. */
+#undef NDEBUG /* Must undef above assert.h or other that might include it. */
+#endif
+
+#include "sox_i.h"
+#include <assert.h>
+#include <string.h>
+
+/* Numerical Recipes cubic spline */
+
+void lsx_prepare_spline3(double const * x, double const * y, int n,
+ double start_1d, double end_1d, double * y_2d)
+{+ double p, qn, sig, un, * u = lsx_malloc((n - 1) * sizeof(*u));
+ int i;
+
+ if (start_1d == HUGE_VAL)
+ y_2d[0] = u[0] = 0; /* Start with natural spline or */
+ else { /* set the start first derivative */+ y_2d[0] = -.5;
+ u[0] = (3 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - start_1d);
+ }
+
+ for (i = 1; i < n - 1; ++i) {+ sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
+ p = sig * y_2d[i - 1] + 2;
+ y_2d[i] = (sig - 1) / p;
+ u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) -
+ (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
+ u[i] = (6 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
+ }
+ if (end_1d == HUGE_VAL)
+ qn = un = 0; /* End with natural spline or */
+ else { /* set the end first derivative */+ qn = .5;
+ un = 3 / (x[n - 1] - x[n - 2]) * (end_1d - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
+ }
+ y_2d[n - 1] = (un - qn * u[n - 2]) / (qn * y_2d[n - 2] + 1);
+ for (i = n - 2; i >= 0; --i)
+ y_2d[i] = y_2d[i] * y_2d[i + 1] + u[i];
+ free(u);
+}
+
+double lsx_spline3(double const * x, double const * y, double const * y_2d,
+ int n, double x1)
+{+ int t, i[2] = {0, 0};+ double d, a, b;
+
+ for (i[1] = n - 1; i[1] - i[0] > 1; t = (i[1] + i[0]) >> 1, i[x[t] > x1] = t);
+ d = x[i[1]] - x[i[0]];
+ assert(d != 0);
+ a = (x[i[1]] - x1) / d;
+ b = (x1 - x[i[0]]) / d;
+ return a * y[i[0]] + b * y[i[1]] +
+ ((a * a * a - a) * y_2d[i[0]] + (b * b * b - b) * y_2d[i[1]]) * d * d / 6;
+}
+
+double lsx_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;
+}
+
+int lsx_set_dft_length(int num_taps) /* Set to 4 x nearest power of 2 */
+{+ int result, n = num_taps;
+ for (result = 8; n > 2; result <<= 1, n >>= 1);
+ result = range_limit(result, 4096, 131072);
+ assert(num_taps * 2 < result);
+ return result;
+}
+
+#include "fft4g.h"
+int * lsx_fft_br;
+double * lsx_fft_sc;
+static int fft_len;
+#ifdef HAVE_OPENMP
+static omp_lock_t fft_cache_lock;
+#endif
+
+void init_fft_cache(void)
+{+ omp_init_lock(&fft_cache_lock);
+}
+
+void clear_fft_cache(void)
+{+ omp_destroy_lock(&fft_cache_lock);
+ free(lsx_fft_br);
+ free(lsx_fft_sc);
+ lsx_fft_sc = NULL;
+ lsx_fft_br = NULL;
+ fft_len = 0;
+}
+
+static void update_fft_cache(int len)
+{+ omp_set_lock(&fft_cache_lock);
+ if (len > fft_len) {+ int old_n = fft_len;
+ fft_len = len;
+ lsx_fft_br = lsx_realloc(lsx_fft_br, dft_br_len(fft_len) * sizeof(*lsx_fft_br));
+ lsx_fft_sc = lsx_realloc(lsx_fft_sc, dft_sc_len(fft_len) * sizeof(*lsx_fft_sc));
+ if (!old_n)
+ lsx_fft_br[0] = 0;
+ }
+}
+
+static sox_bool is_power_of_2(int x)
+{+ return !(x < 2 || (x & (x - 1)));
+}
+
+void lsx_safe_rdft(int len, int type, double * d)
+{+ assert(is_power_of_2(len));
+ update_fft_cache(len);
+ lsx_rdft(len, type, d, lsx_fft_br, lsx_fft_sc);
+ omp_unset_lock(&fft_cache_lock);
+}
+
+void lsx_safe_cdft(int len, int type, double * d)
+{+ assert(is_power_of_2(len));
+ update_fft_cache(len);
+ lsx_cdft(len, type, d, lsx_fft_br, lsx_fft_sc);
+ omp_unset_lock(&fft_cache_lock);
+}
+
+void lsx_power_spectrum(int n, double const * in, double * out)
+{+ int i;
+ double * work = lsx_memdup(in, n * sizeof(*work));
+ lsx_safe_rdft(n, 1, work);
+ out[0] = sqr(work[0]);
+ for (i = 2; i < n; i += 2)
+ out[i >> 1] = sqr(work[i]) + sqr(work[i + 1]);
+ out[i >> 1] = sqr(work[1]);
+ free(work);
+}
+
+void lsx_power_spectrum_f(int n, float const * in, float * out)
+{+ int i;
+ double * work = lsx_malloc(n * sizeof(*work));
+ for (i = 0; i< n; ++i) work[i] = in[i];
+ lsx_safe_rdft(n, 1, work);
+ out[0] = sqr(work[0]);
+ for (i = 2; i < n; i += 2)
+ out[i >> 1] = sqr(work[i]) + sqr(work[i + 1]);
+ out[i >> 1] = sqr(work[1]);
+ free(work);
+}
+
+void lsx_apply_hann_f(float h[], const int num_points)
+{+ int i, m = num_points - 1;
+ for (i = 0; i < num_points; ++i) {+ double x = 2 * M_PI * i / m;
+ h[i] *= .5 - .5 * cos(x);
+ }
+}
+
+void lsx_apply_hann(double h[], const int num_points)
+{+ int i, m = num_points - 1;
+ for (i = 0; i < num_points; ++i) {+ double x = 2 * M_PI * i / m;
+ h[i] *= .5 - .5 * cos(x);
+ }
+}
+
+void lsx_apply_hamming(double h[], const int num_points)
+{+ int i, m = num_points - 1;
+ for (i = 0; i < num_points; ++i) {+ double x = 2 * M_PI * i / m;
+ h[i] *= .53836 - .46164 * cos(x);
+ }
+}
+
+void lsx_apply_bartlett(double h[], const int num_points)
+{+ int i, m = num_points - 1;
+ for (i = 0; i < num_points; ++i) {+ h[i] *= 2. / m * (m / 2. - fabs(i - m / 2.));
+ }
+}
+
+void lsx_apply_blackman(double h[], const int num_points, double alpha /*.16*/)
+{+ int i, m = num_points - 1;
+ for (i = 0; i < num_points; ++i) {+ double x = 2 * M_PI * i / m;
+ h[i] *= (1 - alpha) *.5 - .5 * cos(x) + alpha * .5 * cos(2 * x);
+ }
+}
+
+void lsx_apply_blackman_nutall(double h[], const int num_points)
+{+ int i, m = num_points - 1;
+ for (i = 0; i < num_points; ++i) {+ double x = 2 * M_PI * i / m;
+ h[i] *= .3635819 - .4891775 * cos(x) + .1365995 * cos(2 * x) - .0106411 * cos(3 * x);
+ }
+}
+
+double lsx_kaiser_beta(double att)
+{+ if (att > 100 ) return .1117 * att - 1.11;
+ if (att > 50 ) return .1102 * (att - 8.7);
+ if (att > 20.96) return .58417 * pow(att -20.96, .4) + .07886 * (att - 20.96);
+ return 0;
+}
+
+void lsx_apply_kaiser(double h[], const int num_points, double beta)
+{+ int i, m = num_points - 1;
+ for (i = 0; i <= m; ++i) {+ double x = 2. * i / m - 1;
+ h[i] *= lsx_bessel_I_0(beta * sqrt(1 - x * x)) / lsx_bessel_I_0(beta);
+ }
+}
+
+double * lsx_make_lpf(int num_taps, double Fc, double beta, double scale)
+{+ double * h = malloc(num_taps * sizeof(*h)), sum = 0;
+ int i, m = num_taps - 1;
+ assert(Fc >= 0 && Fc <= 1);
+ 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));
+ if (m - i != i)
+ sum += h[m - i] = h[i];
+ }
+ for (i = 0; i < num_taps; ++i) h[i] *= scale / sum;
+ return h;
+}
+
+double * lsx_design_lpf(
+ double Fp, /* End of pass-band; ~= 0.01dB point */
+ double Fc, /* Start of stop-band */
+ double Fn, /* Nyquist freq; e.g. 0.5, 1, PI */
+ sox_bool allow_aliasing,
+ double att, /* Stop-band attenuation in dB */
+ int * num_taps, /* (Single phase.) 0: value will be estimated */
+ int k) /* Number of phases; 0 for single-phase */
+{+ double tr_bw, beta;
+
+ if (allow_aliasing)
+ Fc += (Fc - Fp) * LSX_TO_3dB;
+ Fp /= Fn, Fc /= Fn; /* Normalise to Fn = 1 */
+ tr_bw = LSX_TO_6dB * (Fc - Fp); /* Transition band-width: 6dB to stop points */
+
+ if (*num_taps == 0) { /* TODO this could be cleaner, esp. for k != 0 */+ double n160 = (.0425* att - 1.4) / tr_bw; /* Half order for att = 160 */
+ int n = n160 * (16.556 / (att - 39.6) + .8625) + .5; /* For att [80,160) */
+ *num_taps = k? 2 * n : 2 * (n + (n & 1)) + 1; /* =1 %4 (0 phase 1/2 band) */
+ }
+ assert(att >= 80);
+ beta = att < 100 ? .1102 * (att - 8.7) : .1117 * att - 1.11;
+ if (k)
+ *num_taps = *num_taps * k - 1;
+ else k = 1;
+ return lsx_make_lpf(*num_taps, (Fc - tr_bw) / k, beta, (double)k);
+}
+
+void lsx_fir_to_phase(double * * h, int * len,
+ int * post_len, double phase0)
+{+ double * work, phase = (phase0 > 50 ? 100 - phase0 : phase0) / 50;
+ int work_len, begin, end, peak = 0, i = *len;
+
+ for (work_len = 32; i > 1; work_len <<= 1, i >>= 1);
+ work = calloc(work_len, sizeof(*work));
+ for (i = 0; i < *len; ++i) work[i] = (*h)[i];
+
+ lsx_safe_rdft(work_len, 1, work); /* Cepstral: */
+ work[0] = log(fabs(work[0])), work[1] = log(fabs(work[1]));
+ for (i = 2; i < work_len; i += 2) {+ work[i] = log(sqrt(sqr(work[i]) + sqr(work[i + 1])));
+ work[i + 1] = 0;
+ }
+ lsx_safe_rdft(work_len, -1, work);
+ for (i = 0; i < work_len; ++i) work[i] *= 2. / work_len;
+ for (i = 1; i < work_len / 2; ++i) { /* Window to reject acausal components */+ work[i] *= 2;
+ work[i + work_len / 2] = 0;
+ }
+ lsx_safe_rdft(work_len, 1, work);
+
+ /* 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. */
+
+ for (i = 2; i < work_len; i += 2) /* Interpolate between linear & min phase */
+ work[i + 1] = phase * M_PI * .5 * i + (1 - phase) * work[i + 1];
+
+ work[0] = exp(work[0]), work[1] = exp(work[1]);
+ for (i = 2; i < work_len; i += 2) {+ double x = exp(work[i]);
+ work[i ] = x * cos(work[i + 1]);
+ work[i + 1] = x * sin(work[i + 1]);
+ }
+ lsx_safe_rdft(work_len, -1, work);
+ for (i = 0; i < work_len; ++i) work[i] *= 2. / work_len;
+
+ for (i = 1; i < work_len; ++i) if (work[i] > work[peak]) /* Find peak pos. */
+ peak = i; /* N.B. peak val. > 0 */
+
+ if (phase == 0)
+ begin = 0;
+ else if (phase == 1)
+ begin = 1 + (work_len - *len) / 2;
+ else {+ if (peak < work_len / 4) { /* Low phases can wrap impulse, so unwrap: */+ memmove(work + work_len / 4, work, work_len / 2 * sizeof(*work));
+ memmove(work, work + work_len * 3 / 4, work_len / 4 * sizeof(*work));
+ peak += work_len / 4;
+ }
+ begin = (.997 - (2 - phase) * .22) * *len + .5;
+ end = (.997 + (0 - phase) * .22) * *len + .5;
+ begin = peak - begin - (begin & 1);
+ end = peak + 1 + end + (end & 1);
+ *len = end - begin;
+ *h = realloc(*h, *len * sizeof(**h));
+ }
+ for (i = 0; i < *len; ++i)
+ (*h)[i] = work[begin + (phase0 > 50 ? *len - 1 - i : i)];
+ *post_len = phase0 > 50 ? peak - begin : begin + *len - (peak + 1);
+ free(work);
+}
--- a/src/loudness.c
+++ b/src/loudness.c
@@ -47,7 +47,7 @@
NUMERIC_PARAMETER(n ,127 ,2047)
} while (0);
p->n = 2 * p->n + 1;
- return argc? lsx_usage(effp) : SOX_SUCCESS;
+ return argc? lsx_usage(effp) : SOX_SUCCESS;
}
static double * make_filter(int n, double start, double delta, double rate)
--- a/src/rate.c
+++ b/src/rate.c
@@ -187,123 +187,6 @@
}
}
-static double * make_lpf(int num_taps, double Fc, double beta, double scale)
-{- double * h = malloc(num_taps * sizeof(*h)), sum = 0;
- int i, m = num_taps - 1;
- assert(Fc >= 0 && Fc <= 1);
- 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));
- if (m - i != i)
- sum += h[m - i] = h[i];
- }
- for (i = 0; i < num_taps; ++i) h[i] *= scale / sum;
- return h;
-}
-
-#define TO_6dB .5869
-#define TO_3dB ((2/3.) * (.5 + TO_6dB))
-#define MAX_TBW0 36.
-#define MAX_TBW0A (MAX_TBW0 / (1 + TO_3dB))
-#define MAX_TBW3 floor(MAX_TBW0 * TO_3dB)
-#define MAX_TBW3A floor(MAX_TBW0A * TO_3dB)
-
-static double * design_lpf(
- double Fp, /* End of pass-band; ~= 0.01dB point */
- double Fc, /* Start of stop-band */
- double Fn, /* Nyquist freq; e.g. 0.5, 1, PI */
- sox_bool allow_aliasing,
- double att, /* Stop-band attenuation in dB */
- int * num_taps, /* (Single phase.) 0: value will be estimated */
- int k) /* Number of phases; 0 for single-phase */
-{- double tr_bw, beta;
-
- if (allow_aliasing)
- Fc += (Fc - Fp) * TO_3dB;
- Fp /= Fn, Fc /= Fn; /* Normalise to Fn = 1 */
- tr_bw = TO_6dB * (Fc - Fp); /* Transition band-width: 6dB to stop points */
-
- if (*num_taps == 0) { /* TODO this could be cleaner, esp. for k != 0 */- double n160 = (.0425* att - 1.4) / tr_bw; /* Half order for att = 160 */
- int n = n160 * (16.556 / (att - 39.6) + .8625) + .5; /* For att [80,160) */
- *num_taps = k? 2 * n : 2 * (n + (n & 1)) + 1; /* =1 %4 (0 phase 1/2 band) */
- }
- assert(att >= 80);
- beta = att < 100 ? .1102 * (att - 8.7) : .1117 * att - 1.11;
- if (k)
- *num_taps = *num_taps * k - 1;
- else k = 1;
- return make_lpf(*num_taps, (Fc - tr_bw) / k, beta, (double)k);
-}
-
-static void fir_to_phase(double * * h, int * len,
- int * post_len, double phase0)
-{- double * work, phase = (phase0 > 50 ? 100 - phase0 : phase0) / 50;
- int work_len, begin, end, peak = 0, i = *len;
-
- for (work_len = 32; i > 1; work_len <<= 1, i >>= 1);
- work = calloc(work_len, sizeof(*work));
- for (i = 0; i < *len; ++i) work[i] = (*h)[i];
-
- lsx_safe_rdft(work_len, 1, work); /* Cepstral: */
- work[0] = log(fabs(work[0])), work[1] = log(fabs(work[1]));
- for (i = 2; i < work_len; i += 2) {- work[i] = log(sqrt(sqr(work[i]) + sqr(work[i + 1])));
- work[i + 1] = 0;
- }
- lsx_safe_rdft(work_len, -1, work);
- for (i = 0; i < work_len; ++i) work[i] *= 2. / work_len;
- for (i = 1; i < work_len / 2; ++i) { /* Window to reject acausal components */- work[i] *= 2;
- work[i + work_len / 2] = 0;
- }
- lsx_safe_rdft(work_len, 1, work);
-
- /* 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. */
-
- for (i = 2; i < work_len; i += 2) /* Interpolate between linear & min phase */
- work[i + 1] = phase * M_PI * .5 * i + (1 - phase) * work[i + 1];
-
- work[0] = exp(work[0]), work[1] = exp(work[1]);
- for (i = 2; i < work_len; i += 2) {- double x = exp(work[i]);
- work[i ] = x * cos(work[i + 1]);
- work[i + 1] = x * sin(work[i + 1]);
- }
- lsx_safe_rdft(work_len, -1, work);
- for (i = 0; i < work_len; ++i) work[i] *= 2. / work_len;
-
- for (i = 1; i < work_len; ++i) if (work[i] > work[peak]) /* Find peak pos. */
- peak = i; /* N.B. peak val. > 0 */
-
- if (phase == 0)
- begin = 0;
- else if (phase == 1)
- begin = 1 + (work_len - *len) / 2;
- else {- if (peak < work_len / 4) { /* Low phases can wrap impulse, so unwrap: */- memmove(work + work_len / 4, work, work_len / 2 * sizeof(*work));
- memmove(work, work + work_len * 3 / 4, work_len / 4 * sizeof(*work));
- peak += work_len / 4;
- }
- begin = (.997 - (2 - phase) * .22) * *len + .5;
- end = (.997 + (0 - phase) * .22) * *len + .5;
- begin = peak - begin - (begin & 1);
- end = peak + 1 + end + (end & 1);
- *len = end - begin;
- *h = realloc(*h, *len * sizeof(**h));
- }
- for (i = 0; i < *len; ++i)
- (*h)[i] = work[begin + (phase0 > 50 ? *len - 1 - i : i)];
- *post_len = phase0 > 50 ? peak - begin : begin + *len - (peak + 1);
- free(work);
-}
-
static void half_band_filter_init(rate_shared_t * p, unsigned which,
int num_taps, sample_t const h[], double Fp, double atten, int multiplier,
double phase, sox_bool allow_aliasing)
@@ -324,10 +207,10 @@
else {/* Adjustment to negate att degradation with intermediate phase */
double att = phase && phase != 50 && phase != 100? atten * (34./33) : atten;
- double * h = design_lpf(Fp, 1., 2., allow_aliasing, att, &num_taps, 0);
+ double * h = lsx_design_lpf(Fp, 1., 2., allow_aliasing, att, &num_taps, 0);
if (phase != 50)
- fir_to_phase(&h, &num_taps, &f->post_peak, phase);
+ lsx_fir_to_phase(&h, &num_taps, &f->post_peak, phase);
else f->post_peak = num_taps / 2;
dft_length = lsx_set_dft_length(num_taps);
@@ -422,8 +305,8 @@
f1 = &f->interp[interp_order];
if (!last_stage.shared->poly_fir_coefs) {int num_taps = 0, phases = divisor == 1? (1 << f1->phase_bits) : divisor;
- raw_coef_t * coefs =
- design_lpf(f->pass, f->stop, 1., sox_false, f->att, &num_taps, phases);
+ raw_coef_t * coefs = lsx_design_lpf(
+ f->pass, f->stop, 1., sox_false, f->att, &num_taps, phases);
assert(num_taps == f->num_coefs * phases - 1);
last_stage.shared->poly_fir_coefs =
prepare_coefs(coefs, f->num_coefs, phases, interp_order, mult);
@@ -445,8 +328,8 @@
{0, NULL, .931, 110}, {0, NULL, .931, 125}, {0, NULL, .931, 170}};filter_t const * f = &filters[quality - Low];
double att = allow_aliasing? (34./33)* f->a : f->a; /* negate att degrade */
- double bw = bandwidth? 1 - (1 - bandwidth / 100) / TO_3dB : f->bw;
- double min = 1 - (allow_aliasing? MAX_TBW0A : MAX_TBW0) / 100;
+ double bw = bandwidth? 1 - (1 - bandwidth / 100) / LSX_TO_3dB : f->bw;
+ double min = 1 - (allow_aliasing? LSX_MAX_TBW0A : LSX_MAX_TBW0) / 100;
assert((size_t)(quality - Low) < array_length(filters));
half_band_filter_init(shared, p->upsample, f->len, f->h, bw, att, mult, phase, allow_aliasing);
if (p->upsample) {@@ -573,7 +456,7 @@
while ((c = getopt(argc, argv, opts)) != -1) switch (c) { GETOPT_NUMERIC('i', coef_interp, 1 , 3) GETOPT_NUMERIC('p', phase, 0 , 100)- GETOPT_NUMERIC('b', bandwidth, 100 - MAX_TBW3, 99.7)+ GETOPT_NUMERIC('b', bandwidth, 100 - LSX_MAX_TBW3, 99.7)case 'M': p->phase = 0; break;
case 'I': p->phase = 25; break;
case 'L': p->phase = 50; break;
@@ -589,8 +472,8 @@
return SOX_EOF;
}
- if (p->bandwidth && p->bandwidth < 100 - MAX_TBW3A && p->allow_aliasing) {- lsx_fail("minimum allowed bandwidth with aliasing is %g%%", 100 - MAX_TBW3A);+ if (p->bandwidth && p->bandwidth < 100 - LSX_MAX_TBW3A && p->allow_aliasing) {+ lsx_fail("minimum allowed bandwidth with aliasing is %g%%", 100 - LSX_MAX_TBW3A);return SOX_EOF;
}
--- a/src/sox.c
+++ b/src/sox.c
@@ -1627,6 +1627,7 @@
"--replay-gain track|album|off Default: off (sox, rec), track (play)",
"-R Use default random numbers (same on each run of SoX)",
"-S, --show-progress Display progress while processing audio data",
+"--single-threaded Disable parallel effects channels processing",
"--temp DIRECTORY Specify the directory to use for temporary files",
"--version Display version number of SoX and exit",
"-V[LEVEL] Increment or set verbosity level (default 2); levels:",
--- a/src/sox_i.h
+++ b/src/sox_i.h
@@ -85,6 +85,8 @@
int lsx_set_dft_length(int num_taps);
extern int * lsx_fft_br;
extern double * lsx_fft_sc;
+void init_fft_cache(void);
+void clear_fft_cache(void);
void lsx_safe_rdft(int len, int type, double * d);
void lsx_safe_cdft(int len, int type, double * d);
void lsx_power_spectrum(int n, double const * in, double * out);
@@ -97,6 +99,23 @@
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_design_lpf(
+ double Fp, /* End of pass-band; ~= 0.01dB point */
+ double Fc, /* Start of stop-band */
+ double Fn, /* Nyquist freq; e.g. 0.5, 1, PI */
+ sox_bool allow_aliasing,
+ double att, /* Stop-band attenuation in dB */
+ int * num_taps, /* (Single phase.) 0: value will be estimated */
+ int k); /* Number of phases; 0 for single-phase */
+void lsx_fir_to_phase(double * * h, int * len,
+ int * post_len, double phase0);
+#define LSX_TO_6dB .5869
+#define LSX_TO_3dB ((2/3.) * (.5 + LSX_TO_6dB))
+#define LSX_MAX_TBW0 36.
+#define LSX_MAX_TBW0A (LSX_MAX_TBW0 / (1 + LSX_TO_3dB))
+#define LSX_MAX_TBW3 floor(LSX_MAX_TBW0 * LSX_TO_3dB)
+#define LSX_MAX_TBW3A floor(LSX_MAX_TBW0A * LSX_TO_3dB)
#ifndef HAVE_STRCASECMP
int strcasecmp(const char *s1, const char *s2);