ref: 4eb9609a1a04a1f72ec472e7362860cb110c5b1a
parent: af2cac8841f5e141028ae556a6b770fd1340b6aa
author: robs <robs>
date: Thu Nov 6 17:23:14 EST 2008
Fix [2223251] mcompand should use linkwitz-riley
--- a/ChangeLog
+++ b/ChangeLog
@@ -67,6 +67,7 @@
o Fix rare crash with `rate' effect. (robs)
o Fix [2190767] `norm' under-amplifying in some cases. (George Yohng)
o Fix [2007062] Earwax effect can overflow; should clip. (robs)
+ o Fix [2223251] mcompand should use linkwitz-riley. (robs)
o Fix `phaser' clicks and overflows. (robs)
o Trim will now skip past 2G point correctly. (cbagwell)
o Improved handling of speed changes in the effects chain. (robs)
--- a/sox.1
+++ b/sox.1
@@ -1777,14 +1777,14 @@
\fBmcompand\fR \(dq\fIattack1\fB,\fIdecay1\fR{\fB,\fIattack2\fB,\fIdecay2\fR}
[\fIsoft-knee-dB\fB:\fR]\fIin-dB1\fR[\fB,\fIout-dB1\fR]{\fB,\fIin-dB2\fB,\fIout-dB2\fR}
.br
-[\fIgain\fR [\fIinitial-volume-dB\fR [\fIdelay\fR]]]\(dq {\fIxover-freq\fR[\fBk\fR] \(dqattack1,...\(dq}
+[\fIgain\fR [\fIinitial-volume-dB\fR [\fIdelay\fR]]]\(dq {\fIcrossover-freq\fR[\fBk\fR] \(dqattack1,...\(dq}
.SP
The multi-band compander is similar to the single-band compander but the
-audio is first divided into bands using Butterworth cross-over filters
+audio is first divided into bands using Linkwitz-Riley cross-over filters
and a separately specifiable compander run on each band. See the
\fBcompand\fR effect for the definition of its parameters. Compand
parameters are specified between double quotes and the crossover
-frequency for that band is given by \fIxover-freq\fR; these can be
+frequency for that band is given by \fIcrossover-freq\fR; these can be
repeated to create multiple bands.
.SP
For example, the following (one long) command shows how multi-band
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -261,15 +261,14 @@
libsox_la_SOURCES += \
band.h bend.c biquad.c biquad.h biquads.c chorus.c compand.c \
compandt.c compandt.h contrast.c dcshift.c delay.c dither.c \
- dither_fir.h dither_iir.h earwax.c \
- echo.c echos.c effects.c effects.h effects_i.c fade.c fft4g.c fft4g.h \
- fifo.h filter.c flanger.c input.c ladspa.c loudness.c mcompand.c \
- mixer.c noiseprof.c noisered.c noisered.h normalise.c output.c pad.c \
- pan.c phaser.c 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
+ dither_fir.h dither_iir.h earwax.c echo.c echos.c effects.c effects.h \
+ effects_i.c fade.c fft4g.c fft4g.h fifo.h filter.c flanger.c input.c \
+ ladspa.c loudness.c mcompand.c mcompand_xover.h mixer.c noiseprof.c \
+ noisered.c noisered.h normalise.c output.c pad.c pan.c phaser.c \
+ 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 shape.c
if HAVE_PNG
libsox_la_SOURCES += spectrogram.c
endif
--- a/src/mcompand.c
+++ b/src/mcompand.c
@@ -3,23 +3,12 @@
*
* Compander code adapted from the SoX compand effect, by Nick Bailey
*
- * Butterworth code adapted from the SoX Butterworth effect family, by
- * Jan Paul Schmidt
+ * SoX is Copyright 1999 Chris Bagwell And Nick Bailey This source code is
+ * freely redistributable and may be used for any purpose. This copyright
+ * notice must be maintained. Chris Bagwell And Nick Bailey are not
+ * responsible for the consequences of using this software.
*
- * SoX is Copyright 1999 Chris Bagwell And Nick Bailey
- * This source code is freely redistributable and may be used for
- * any purpose. This copyright notice must be maintained.
- * Chris Bagwell And Nick Bailey are not responsible for
- * the consequences of using this software.
- */
-
-#include "sox_i.h"
-
-#include <string.h>
-#include <stdlib.h>
-#include "compandt.h"
-
-/*
+ *
* Usage:
* mcompand quoted_compand_args [crossover_frequency
* quoted_compand_args [...]]
@@ -32,164 +21,34 @@
*
* Beware a variety of headroom (clipping) bugaboos.
*
- * Here is an example application, an FM radio sound simulator (or
- * broadcast signal conditioner, if the lowp at the end is skipped -
- * note that the pipeline is set up with US-style 75us preemphasis).
+ * Implementation details:
+ * The input is divided into bands using 4th order Linkwitz-Riley IIRs.
+ * This is akin to the crossover of a loudspeaker, and results in flat
+ * frequency response absent compander action.
*
- * sox -V -t raw -r 44100 -s -w -c 2 - -t raw -r 44100 -s -l -c 2 \
- * - vol -3 db filter 8000- 32 100 mcompand ".005,.1 \
- * -47,-40,-34,-34,-17,-33 0 0 0" 100 ".003,.05 \
- * -47,-40,-34,-34,-17,-33 0 0 0" 400 ".000625,.0125 \
- * -47,-40,-34,-34,-15,-33 0 0 0" 1600 ".0001,.025 \
- * -47,-40,-34,-34,-31,-31,-0,-30 0 0 0" 6400 \
- * "0,.025 -38,-31,-28,-28,-0,-25 0 0 0" vol 27 db vol -12 \
- * db highpass 22 highpass 22 filter -17500 256 vol +12 db \
- * vol -3 db lowp 1780
+ * The outputs of the array of companders is summed, and sample truncation
+ * is done on the final sum.
*
- * implementation details:
- *
- * The input is divided into bands using aligned 2nd order
- * Butterworth IIR's (code adapted from SoX's existing Butterworth
- * effect). This is akin to the crossover of a loudspeaker, and
- * results in flat frequency response (within a db or so) absent
- * compander action.
- *
- * (Crossover design from "Electroacoustic System Design" by
- * Tom Raymann, http://www.traymann.free-online.co.uk/soph.htm)
- *
- * The imported Butterworth code is modified to handle
- * interleaved-channel sample streams, to make the code clean and
- * efficient. The outputs of the array of companders is summed, and
- * sample truncation is done on the final sum.
- *
- * Modifications to the predictive compression code properly
- * maintain alignment of the outputs of the array of companders when
- * the companders have different prediction intervals (volume
- * application delays). Note that the predictive mode of the
- * limiter needs some TLC - in fact, a rewrite - since what's really
- * useful is to assure that a waveform won't be clipped, by slewing
- * the volume in advance so that the peak is at limit (or below, if
- * there's a higher subsequent peak visible in the lookahead window)
- * once it's reached. */
+ * Modifications to the predictive compression code properly maintain
+ * alignment of the outputs of the array of companders when the companders
+ * have different prediction intervals (volume application delays). Note
+ * that the predictive mode of the limiter needs some TLC - in fact, a
+ * rewrite - since what's really useful is to assure that a waveform won't
+ * be clipped, by slewing the volume in advance so that the peak is at
+ * limit (or below, if there's a higher subsequent peak visible in the
+ * lookahead window) once it's reached. */
-struct xy {
- double x [2];
- double y [2];
-} ;
-
-typedef struct {
- struct xy *xy_low, *xy_high;
-
- double a_low[3], a_high[3];
- double b_low[2], b_high[3];
-
- /*
- * Cut off frequencies for respective filters
- */
- double frequency_low, frequency_high;
-
- double bandwidth;
-} butterworth_crossover_t;
-
-static int lowpass_setup (butterworth_crossover_t * butterworth, double frequency, sox_rate_t rate, size_t nchan) {
- double c;
-
- butterworth->xy_low = lsx_calloc(nchan, sizeof(struct xy));
- butterworth->xy_high = lsx_calloc(nchan, sizeof(struct xy));
-
- /* lowpass setup */
- butterworth->frequency_low = frequency/1.3;
-
- c = 1.0 / tan (M_PI * butterworth->frequency_low / rate);
-
- butterworth->a_low[0] = 1.0 / (1.0 + sqrt(2.0) * c + c * c);
- butterworth->a_low[1] = 2.0 * butterworth->a_low [0];
- butterworth->a_low[2] = butterworth->a_low [0];
-
- butterworth->b_low[0] = 2 * (1.0 - c * c) * butterworth->a_low[0];
- butterworth->b_low[1] = (1.0 - sqrt(2.0) * c + c * c) * butterworth->a_low[0];
-
- /* highpass setup */
- butterworth->frequency_high = frequency*1.3;
+#ifdef NDEBUG /* Enable assert always. */
+#undef NDEBUG /* Must undef above assert.h or other that might include it. */
+#endif
- c = tan (M_PI * butterworth->frequency_high / rate);
+#include "sox_i.h"
+#include <assert.h>
+#include <string.h>
+#include <stdlib.h>
+#include "compandt.h"
+#include "mcompand_xover.h"
- butterworth->a_high[0] = 1.0 / (1.0 + sqrt (2.0) * c + c * c);
- butterworth->a_high[1] = -2.0 * butterworth->a_high[0];
- butterworth->a_high[2] = butterworth->a_high[0];
-
- butterworth->b_high[0] = 2 * (c * c - 1.0) * butterworth->a_high[0];
- butterworth->b_high[1] = (1.0 - sqrt(2.0) * c + c * c) * butterworth->a_high[0];
-
- return (SOX_SUCCESS);
-}
-
-static int lowpass_flow(sox_effect_t * effp, butterworth_crossover_t * butterworth, size_t nChan, sox_sample_t *ibuf, sox_sample_t *lowbuf, sox_sample_t *highbuf,
- size_t len) {
- size_t chan;
- double in, out;
-
- size_t done;
-
- sox_sample_t *ibufptr, *lowbufptr, *highbufptr;
-
- for (chan=0;chan<nChan;++chan) {
- ibufptr = ibuf+chan;
- lowbufptr = lowbuf+chan;
- highbufptr = highbuf+chan;
-
- for (done = chan; done < len; done += nChan) {
- in = *ibufptr;
- ibufptr += nChan;
-
- /*
- * Substituting butterworth->a [x] and butterworth->b [x] with
- * variables, which are set outside of the loop, did not increased
- * speed on my AMD Box. GCC seems to do a good job :o)
- */
-
- out =
- butterworth->a_low[0] * in +
- butterworth->a_low [1] * butterworth->xy_low[chan].x [0] +
- butterworth->a_low [2] * butterworth->xy_low[chan].x [1] -
- butterworth->b_low [0] * butterworth->xy_low[chan].y [0] -
- butterworth->b_low [1] * butterworth->xy_low[chan].y [1];
-
- butterworth->xy_low[chan].x [1] = butterworth->xy_low[chan].x [0];
- butterworth->xy_low[chan].x [0] = in;
- butterworth->xy_low[chan].y [1] = butterworth->xy_low[chan].y [0];
- butterworth->xy_low[chan].y [0] = out;
-
- SOX_SAMPLE_CLIP_COUNT(out, effp->clips);
-
- *lowbufptr = out;
-
- out =
- butterworth->a_high[0] * in +
- butterworth->a_high [1] * butterworth->xy_high[chan].x [0] +
- butterworth->a_high [2] * butterworth->xy_high[chan].x [1] -
- butterworth->b_high [0] * butterworth->xy_high[chan].y [0] -
- butterworth->b_high [1] * butterworth->xy_high[chan].y [1];
-
- butterworth->xy_high[chan].x [1] = butterworth->xy_high[chan].x [0];
- butterworth->xy_high[chan].x [0] = in;
- butterworth->xy_high[chan].y [1] = butterworth->xy_high[chan].y [0];
- butterworth->xy_high[chan].y [0] = out;
-
- SOX_SAMPLE_CLIP_COUNT(out, effp->clips);
-
- /* don't forget polarity reversal of high pass! */
-
- *highbufptr = -out;
-
- lowbufptr += nChan;
- highbufptr += nChan;
- }
- }
-
- return (SOX_SUCCESS);
-}
-
typedef struct {
sox_compandt_t transfer_fn;
@@ -200,7 +59,7 @@
double *volume; /* Current "volume" of each channel */
double delay; /* Delay to apply before companding */
double topfreq; /* upper bound crossover frequency */
- butterworth_crossover_t filter;
+ crossover_t filter;
sox_sample_t *delay_buf; /* Old samples, used for delay processing */
size_t delay_size; /* lookahead for this band (in samples) - function of delay, above */
ptrdiff_t delay_buf_ptr; /* Index into delay_buf */
@@ -314,7 +173,7 @@
/* how many bands? */
if (! (n&1)) {
lsx_fail("mcompand accepts only an odd number of arguments:\n"
- " mcompand quoted_compand_args [xover_freq quoted_compand_args [...]");
+ " mcompand quoted_compand_args [crossover_freq quoted_compand_args [...]");
return SOX_EOF;
}
c->nBands = (n+1)>>1;
@@ -387,7 +246,7 @@
l->delay_buf_cnt = 0;
if (l->topfreq != 0)
- lowpass_setup(&l->filter, l->topfreq, effp->out_signal.rate, (size_t) effp->out_signal.channels);
+ crossover_setup(effp, &l->filter, l->topfreq);
}
return (SOX_SUCCESS);
}
@@ -497,7 +356,7 @@
ibuf_copy = lsx_malloc(*isamp * sizeof(sox_sample_t));
memcpy(ibuf_copy, ibuf, *isamp * sizeof(sox_sample_t));
- /* split ibuf into bands using butterworths, pipe each band through sox_mcompand_flow_1, then add back together and write to obuf */
+ /* split ibuf into bands using filters, pipe each band through sox_mcompand_flow_1, then add back together and write to obuf */
memset(obuf,0,len * sizeof *obuf);
for (band=0,abuf=ibuf_copy,bbuf=c->band_buf2,cbuf=c->band_buf1;band<c->nBands;++band) {
@@ -504,7 +363,7 @@
l = &c->bands[band];
if (l->topfreq)
- lowpass_flow(effp, &l->filter, (size_t) effp->out_signal.channels, abuf, bbuf, cbuf, len);
+ crossover_flow(effp, &l->filter, abuf, bbuf, cbuf, len);
else {
bbuf = abuf;
abuf = cbuf;
@@ -595,10 +454,8 @@
for (band = 0; band < c->nBands; band++) {
l = &c->bands[band];
free(l->delay_buf);
- if (l->topfreq != 0) {
- free(l->filter.xy_low);
- free(l->filter.xy_high);
- }
+ if (l->topfreq != 0)
+ free(l->filter.previous);
}
return SOX_SUCCESS;
--- /dev/null
+++ b/src/mcompand_xover.h
@@ -1,0 +1,106 @@
+/* libSoX Compander Crossover Filter (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
+ */
+
+#define N 4 /* 4th order Linkwitz-Riley IIRs */
+#define CONVOLVE _ _ _ _
+
+typedef struct {double in, out_low, out_high;} previous_t[N * 2];
+
+typedef struct {
+ previous_t * previous;
+ size_t pos;
+ double coefs[3 *(N+1)];
+} crossover_t;
+
+static void square_quadratic(char const * name, double const * x, double * y)
+{
+ assert(N == 4);
+ y[0] = x[0] * x[0];
+ y[1] = 2 * x[0] * x[1];
+ y[2] = 2 * x[0] * x[2] + x[1] * x[1];
+ y[3] = 2 * x[1] * x[2];
+ y[4] = x[2] * x[2];
+ lsx_debug("%s=[%.16g %.16g %.16g %.16g %.16g];", name,
+ y[0], y[1], y[2], y[3], y[4]);
+}
+
+static int crossover_setup(sox_effect_t * effp, crossover_t * p, double frequency)
+{
+ double w0 = 2 * M_PI * frequency / effp->in_signal.rate;
+ double Q = sqrt(.5), alpha = sin(w0)/(2*Q);
+ double x[9], norm;
+ int i;
+
+ if (w0 > M_PI) {
+ lsx_fail("frequency must not exceed half the sample-rate (Nyquist rate)");
+ return SOX_EOF;
+ }
+ x[0] = (1 - cos(w0))/2; /* Cf. filter_LPF in biquads.c */
+ x[1] = 1 - cos(w0);
+ x[2] = (1 - cos(w0))/2;
+ x[3] = (1 + cos(w0))/2; /* Cf. filter_HPF in biquads.c */
+ x[4] = -(1 + cos(w0));
+ x[5] = (1 + cos(w0))/2;
+ x[6] = 1 + alpha;
+ x[7] = -2*cos(w0);
+ x[8] = 1 - alpha;
+ for (norm = x[6], i = 0; i < 9; ++i) x[i] /= norm;
+ square_quadratic("lb", x , p->coefs);
+ square_quadratic("hb", x + 3, p->coefs + 5);
+ square_quadratic("a" , x + 6, p->coefs + 10);
+
+ p->previous = lsx_calloc(effp->in_signal.channels, sizeof(*p->previous));
+ return SOX_SUCCESS;
+}
+
+static int crossover_flow(sox_effect_t * effp, crossover_t * p, sox_sample_t
+ *ibuf, sox_sample_t *obuf_low, sox_sample_t *obuf_high, size_t len0)
+{
+ double out_low, out_high;
+ size_t c, len = len0 / effp->in_signal.channels;
+ assert(len * effp->in_signal.channels == len0);
+
+ while (len--) {
+ p->pos = p->pos? p->pos - 1 : N - 1;
+ for (c = 0; c < effp->in_signal.channels; ++c) {
+#define _ out_low += p->coefs[j] * p->previous[c][p->pos + j].in \
+ - p->coefs[2*N+2 + j] * p->previous[c][p->pos + j].out_low, ++j;
+ {
+ int j = 1;
+ out_low = p->coefs[0] * *ibuf;
+ CONVOLVE
+ assert(j == N+1);
+ *obuf_low++ = SOX_ROUND_CLIP_COUNT(out_low, effp->clips);
+ }
+#undef _
+#define _ out_high += p->coefs[j+N+1] * p->previous[c][p->pos + j].in \
+ - p->coefs[2*N+2 + j] * p->previous[c][p->pos + j].out_high, ++j;
+ {
+ int j = 1;
+ out_high = p->coefs[N+1] * *ibuf;
+ CONVOLVE
+ assert(j == N+1);
+ *obuf_high++ = SOX_ROUND_CLIP_COUNT(out_high, effp->clips);
+ }
+ p->previous[c][p->pos + N].in = p->previous[c][p->pos].in = *ibuf++;
+ p->previous[c][p->pos + N].out_low = p->previous[c][p->pos].out_low = out_low;
+ p->previous[c][p->pos + N].out_high = p->previous[c][p->pos].out_high = out_high;
+ }
+ }
+ return SOX_SUCCESS;
+}
+