shithub: sox

Download patch

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;
+}
+