ref: 2ce9ae35fc417bb9fb66ce7ca3dec12135bfc7be
parent: fe37415768718214125ee0a4cf9ae8d3237755ea
author: robs <robs>
date: Sun Sep 21 16:46:07 EDT 2008
new loudness effect
--- a/AUTHORS
+++ b/AUTHORS
@@ -17,7 +17,7 @@
EFFECTS: tempo, pad, bass, treble, delay, new reverb, new
flanger, soft-knee companding, speed via resampling, filters
makeover inc. gnuplot & octave plotting, splice, remix, norm,
- contrast, new rate, spectrogram, new pitch, riaa;
+ contrast, new rate, spectrogram, new pitch, riaa, loudness;
new effects chain with buffering and any # channels.
OTHERS: open input files via URL, file merge, play
multiple files with mixed format types, play with replay-gain,
--- a/ChangeLog
+++ b/ChangeLog
@@ -44,7 +44,9 @@
Effects:
- o New `riaa' effect; RIAA vinyl playback EQ. (robs)
+ o New `riaa' effect: RIAA vinyl playback EQ. (robs)
+ o New `loudness' effect: gain control with ISO 226 loudness
+ compensation. (robs)
o New -b option for the norm effect; can be used to fix stereo
imbalance. (robs)
o Change default settings for `rate' effect; add -s -L options to
--- a/README
+++ b/README
@@ -127,6 +127,7 @@
o dcshift Apply or remove DC offset
o fade Apply a fade-in and/or fade-out to the audio
o gain Apply gain or attenuation
+ o loudness Gain control with ISO 226 loudness compensation
o mcompand Multi-band compression/expansion/limiting
o norm Normalise to 0dB (or other), or fix imbalance
o vol Adjust audio volume
--- a/soxeffect.7
+++ b/soxeffect.7
@@ -652,6 +652,23 @@
and one audio output port can be used. If found, the environment varible
LADSPA_PATH will be used as search path for plugins.
.TP
+\fBloudness [\fIgain\fR [\fIreference\fR]]
+Loudness control\*msimilar to the
+.B gain
+effect, but provides equalisation for the human auditory system. See
+http://en.wikipedia.org/wiki/Loudness for a detailed description of
+loudness. The gain is adjusted by the given
+.I gain
+parameter (usually negative) and the signal equalised according to ISO
+226 w.r.t. a reference level of 65dB, though an alternative
+.I reference
+level may be given if the original audio has been equalised for some
+other optimal level.
+.SP
+See also the
+.B gain
+effect.
+.TP
\fBlowpass\fR [\fB\-1\fR|\fB\-2\fR] \fIfrequency\fR[\fBk\fR]\fR [\fRwidth\fR[\fBq\fR\^|\^\fBo\fR\^|\^\fBh\fR\^|\^\fBk\fR]]
Apply a low-pass filter.
See the description of the \fBhighpass\fR effect for details.
@@ -967,7 +984,7 @@
audio. If no quality option is given, the quality level used is `high'.
.SP
The `quick' algorithm uses cubic interpolation; all others use
-band-width limited interpolation. The `quick' and `low' quality
+band-limited interpolation. The `quick' and `low' quality
algorithms have a `linear' phase response; for `medium', `high' and
`very high', the phase response is configurable (see below), but
defaults to `intermediate'.
@@ -1010,7 +1027,7 @@
Any phase response (0 = minimum, 25 = intermediate, 50 = linear, 100 = maximum)
T}
\-s Steep filter (band-width = 99%)
-\-b\ 74\-99.7 Any band-width %
+\-b\ 74\-99\*d7 Any band-width %
\-a Allow aliasing above the pass-band
.TE
.DT
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -28,15 +28,16 @@
# Format with: !xargs echo|tr ' ' '\n'|sort|column|expand|sed 's/^/ /'
set(effects_srcs
- biquad earwax mixer polyphas speed
- biquads echo noiseprof rate splice
- chorus echos noisered remix stat
- compand fade normalise repeat swap
- compandt fft4g output resample synth
- contrast filter pad reverb tempo
- dcshift flanger pan reverse tremolo
- delay input phaser silence trim
- dither mcompand pitch skeleff vol
+ biquad echo noiseprof remix swap
+ biquads echos noisered repeat synth
+ chorus fade normalise resample tempo
+ compand fft4g output reverb tremolo
+ compandt filter pad reverse trim
+ contrast flanger pan silence vol
+ dcshift input phaser skeleff
+ delay loudness pitch speed
+ dither mcompand polyphas splice
+ earwax mixer rate stat
)
set(formats_srcs
8svx dat ima-fmt s3-fmt u4-fmt
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -252,7 +252,7 @@
band.h biquad.c biquad.h biquads.c chorus.c compand.c compandt.c \
compandt.h contrast.c dcshift.c delay.c dither.c 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 mcompand.c \
+ 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 pitch.c polyphas.c rabbit.c rate.c \
rate_filters.h rate_half_fir.h rate_poly_fir0.h rate_poly_fir.h \
--- a/src/effects.h
+++ b/src/effects.h
@@ -41,6 +41,7 @@
#ifdef HAVE_LADSPA_H
EFFECT(ladspa)
#endif
+ EFFECT(loudness)
EFFECT(lowpass)
EFFECT(mcompand)
EFFECT(mixer)
--- a/src/effects_i.c
+++ b/src/effects_i.c
@@ -68,6 +68,56 @@
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;
+}
+
enum_item const lsx_wave_enum[] = {ENUM_ITEM(SOX_WAVE_,SINE)
ENUM_ITEM(SOX_WAVE_,TRIANGLE)
@@ -279,6 +329,15 @@
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;
@@ -381,6 +440,24 @@
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)
+{+ 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);
}
}
--- /dev/null
+++ b/src/loudness.c
@@ -1,0 +1,221 @@
+/* Effect: loudness filter 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
+ */
+
+#include "sox_i.h"
+#include "fft4g.h"
+#define FIFO_SIZE_T int
+#include "fifo.h"
+#include <string.h>
+
+typedef struct {+ int dft_length, num_taps;
+ double * coefs;
+} filter_t;
+
+typedef struct {+ double delta, start;
+ int n;
+ size_t samples_in, samples_out;
+ fifo_t input_fifo, output_fifo;
+ filter_t filter, * filter_ptr;
+} priv_t;
+
+static int create(sox_effect_t * effp, int argc, char **argv)
+{+ priv_t * p = (priv_t *)effp->priv;
+ p->filter_ptr = &p->filter;
+ p->delta = -10;
+ p->start = 65;
+ p->n = 1023;
+ do { /* break-able block */+ NUMERIC_PARAMETER(delta,-50 , 15)
+ NUMERIC_PARAMETER(start, 50 , 75)
+ NUMERIC_PARAMETER(n ,127 ,2047)
+ } while (0);
+ p->n = 2 * p->n + 1;
+ return argc? lsx_usage(effp) : SOX_SUCCESS;
+}
+
+static double * make_filter(int n, double start, double delta, double rate)
+{+ static const struct {double f, af, lu, tf;} iso226_table[] = {+ { 20,0.532,-31.6,78.5},{ 25,0.506,-27.2,68.7},{ 31.5,0.480,-23.0,59.5},+ { 40,0.455,-19.1,51.1},{ 50,0.432,-15.9,44.0},{ 63,0.409,-13.0,37.5},+ { 80,0.387,-10.3,31.5},{ 100,0.367, -8.1,26.5},{ 125,0.349, -6.2,22.1},+ { 160,0.330, -4.5,17.9},{ 200,0.315, -3.1,14.4},{ 250,0.301, -2.0,11.4},+ { 315,0.288, -1.1, 8.6},{ 400,0.276, -0.4, 6.2},{ 500,0.267, 0.0, 4.4},+ { 630,0.259, 0.3, 3.0},{ 800,0.253, 0.5, 2.2},{ 1000,0.250, 0.0, 2.4},+ { 1250,0.246, -2.7, 3.5},{ 1600,0.244, -4.1, 1.7},{ 2000,0.243, -1.0,-1.3},+ { 2500,0.243, 1.7,-4.2},{ 3150,0.243, 2.5,-6.0},{ 4000,0.242, 1.2,-5.4},+ { 5000,0.242, -2.1,-1.5},{ 6300,0.245, -7.1, 6.0},{ 8000,0.254,-11.2,12.6},+ {10000,0.271,-10.7,13.9},{12500,0.301, -3.1,12.3},+ };
+ #define LEN (array_length(iso226_table) + 2)
+ #define SPL(phon, t) (10 / t.af * log10(4.47e-3 * (pow(10., .025 * (phon)) - \
+ 1.15) + pow(.4 * pow(10., (t.tf + t.lu) / 10 - 9), t.af)) - t.lu + 94)
+ double fs[LEN], spl[LEN], d[LEN], * work, * h;
+ int i, work_len;
+
+ fs[0] = log(1.);
+ spl[0] = delta * .2;
+ for (i = 0; i < (int)LEN - 2; ++i) {+ spl[i + 1] = SPL(start + delta, iso226_table[i]) -
+ SPL(start , iso226_table[i]);
+ fs[i + 1] = log(iso226_table[i].f);
+ }
+ fs[i + 1] = log(100000.);
+ spl[i + 1] = spl[0];
+ lsx_prepare_spline3(fs, spl, LEN, HUGE_VAL, HUGE_VAL, d);
+
+ for (work_len = 8192; work_len < rate / 2; work_len <<= 1);
+ work = lsx_calloc(work_len, sizeof(*work));
+ h = lsx_calloc(n, sizeof(*h));
+
+ for (i = 1; i <= work_len / 2; ++i) {+ double f = rate * i / work_len;
+ double spl1 = f < 1? spl[0] : lsx_spline3(fs, spl, d, LEN, log(f));
+ work[i < work_len / 2 ? 2 * i : 1] = dB_to_linear(spl1);
+ }
+ lsx_safe_rdft(work_len, -1, work);
+ for (i = 0; i < n; ++i)
+ h[i] = work[(work_len - n / 2 + i) % work_len] * 2. / work_len;
+ lsx_apply_kaiser(h, n, lsx_kaiser_beta(40 + 2./3 * fabs(delta)));
+
+ free(work);
+ return h;
+ #undef SPL
+ #undef LEN
+}
+
+static int start(sox_effect_t * effp)
+{+ priv_t * p = (priv_t *) effp->priv;
+ int i, half = p->n / 2, dft_length = lsx_set_dft_length(p->n);
+ double * h;
+
+ if (p->delta == 0)
+ return SOX_EFF_NULL;
+
+ if (!p->filter_ptr->num_taps) {+ h = make_filter(p->n, p->start, p->delta, effp->in_signal.rate);
+ p->filter_ptr->coefs = lsx_calloc(dft_length, sizeof(*p->filter_ptr->coefs));
+ for (i = 0; i < p->n; ++i)
+ p->filter_ptr->coefs[(i + dft_length - p->n + 1) & (dft_length - 1)]
+ = h[i] / dft_length * 2;
+ free(h);
+ p->filter_ptr->num_taps = p->n;
+ p->filter_ptr->dft_length = dft_length;
+ lsx_safe_rdft(dft_length, 1, p->filter_ptr->coefs);
+ }
+ fifo_create(&p->input_fifo, (int)sizeof(double));
+ memset(fifo_reserve(&p->input_fifo, half), 0, sizeof(double) * half);
+ fifo_create(&p->output_fifo, (int)sizeof(double));
+ return SOX_SUCCESS;
+}
+
+static void filter(priv_t * p)
+{+ int i, num_in = max(0, fifo_occupancy(&p->input_fifo));
+ filter_t const * f = p->filter_ptr;
+ int const overlap = f->num_taps - 1;
+ double * output;
+
+ while (num_in >= f->dft_length) {+ double const * input = fifo_read_ptr(&p->input_fifo);
+ fifo_read(&p->input_fifo, f->dft_length - overlap, NULL);
+ num_in -= f->dft_length - overlap;
+
+ output = fifo_reserve(&p->output_fifo, f->dft_length);
+ fifo_trim_by(&p->output_fifo, overlap);
+ memcpy(output, input, f->dft_length * sizeof(*output));
+
+ lsx_rdft(f->dft_length, 1, output, lsx_fft_br, lsx_fft_sc);
+ output[0] *= f->coefs[0];
+ output[1] *= f->coefs[1];
+ for (i = 2; i < f->dft_length; i += 2) {+ double tmp = output[i];
+ output[i ] = f->coefs[i ] * tmp - f->coefs[i+1] * output[i+1];
+ output[i+1] = f->coefs[i+1] * tmp + f->coefs[i ] * output[i+1];
+ }
+ lsx_rdft(f->dft_length, -1, output, lsx_fft_br, lsx_fft_sc);
+ }
+}
+
+static int flow(sox_effect_t * effp, const sox_sample_t * ibuf,
+ sox_sample_t * obuf, size_t * isamp, size_t * osamp)
+{+ priv_t * p = (priv_t *)effp->priv;
+ size_t i, odone = min(*osamp, (size_t)fifo_occupancy(&p->output_fifo));
+ double const * s = fifo_read(&p->output_fifo, (int)odone, NULL);
+
+ for (i = 0; i < odone; ++i)
+ *obuf++ = SOX_FLOAT_64BIT_TO_SAMPLE(*s++, effp->clips);
+ p->samples_out += odone;
+
+ if (*isamp && odone < *osamp) {+ double * t = fifo_write(&p->input_fifo, (int)*isamp, NULL);
+ p->samples_in += (int)*isamp;
+
+ for (i = *isamp; i; --i)
+ *t++ = SOX_SAMPLE_TO_FLOAT_64BIT(*ibuf++, effp->clips);
+ filter(p);
+ }
+ else *isamp = 0;
+ *osamp = odone;
+ return SOX_SUCCESS;
+}
+
+static int drain(sox_effect_t * effp, sox_sample_t * obuf, size_t * osamp)
+{+ priv_t * p = (priv_t *)effp->priv;
+ static size_t isamp = 0;
+ size_t samples_out = p->samples_in;
+ size_t remaining = samples_out - p->samples_out;
+ double * buff = lsx_calloc(1024, sizeof(*buff));
+
+ if ((int)remaining > 0) {+ while ((size_t)fifo_occupancy(&p->output_fifo) < remaining) {+ fifo_write(&p->input_fifo, 1024, buff);
+ p->samples_in += 1024;
+ filter(p);
+ }
+ fifo_trim_to(&p->output_fifo, (int)remaining);
+ p->samples_in = 0;
+ }
+ free(buff);
+ return flow(effp, 0, obuf, &isamp, osamp);
+}
+
+static int stop(sox_effect_t * effp)
+{+ priv_t * p = (priv_t *) effp->priv;
+
+ fifo_delete(&p->input_fifo);
+ fifo_delete(&p->output_fifo);
+ free(p->filter_ptr->coefs);
+ memset(p->filter_ptr, 0, sizeof(*p->filter_ptr));
+ return SOX_SUCCESS;
+}
+
+sox_effect_handler_t const * sox_loudness_effect_fn(void)
+{+ static sox_effect_handler_t handler = {+ "loudness", "[gain [ref]]", 0,
+ create, start, flow, drain, stop, NULL, sizeof(priv_t)
+ };
+ return &handler;
+}
--- a/src/rate.c
+++ b/src/rate.c
@@ -301,15 +301,6 @@
free(work);
}
-static int 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;
-}
-
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)
@@ -320,7 +311,7 @@
if (f->num_taps)
return;
if (h) {- dft_length = set_dft_length(num_taps);
+ dft_length = lsx_set_dft_length(num_taps);
f->coefs = calloc(dft_length, sizeof(*f->coefs));
for (i = 0; i < num_taps; ++i)
f->coefs[(i + dft_length - num_taps + 1) & (dft_length - 1)]
@@ -337,7 +328,7 @@
fir_to_phase(&h, &num_taps, &f->post_peak, phase);
else f->post_peak = num_taps / 2;
- dft_length = set_dft_length(num_taps);
+ dft_length = lsx_set_dft_length(num_taps);
f->coefs = calloc(dft_length, sizeof(*f->coefs));
for (i = 0; i < num_taps; ++i)
f->coefs[(i + dft_length - num_taps + 1) & (dft_length - 1)]
--- a/src/sox_i.h
+++ b/src/sox_i.h
@@ -61,7 +61,13 @@
unsigned lsx_gcd(unsigned a, unsigned b);
unsigned lsx_lcm(unsigned a, unsigned b);
+void lsx_prepare_spline3(double const * x, double const * y, int n,
+ double start_1d, double end_1d, double * y_2d);
+double lsx_spline3(double const * x, double const * y, double const * y_2d,
+ int n, double x1);
+
double lsx_bessel_I_0(double x);
+int lsx_set_dft_length(int num_taps);
extern int * lsx_fft_br;
extern double * lsx_fft_sc;
void lsx_safe_rdft(int len, int type, double * d);
@@ -72,6 +78,8 @@
void lsx_apply_hann(double h[], const int num_points);
void lsx_apply_hamming(double h[], const int num_points);
void lsx_apply_bartlett(double h[], const int num_points);
+void lsx_apply_blackman(double h[], const int num_points, double alpha);
+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);
--- a/src/test-srcs
+++ b/src/test-srcs
@@ -20,15 +20,7 @@
# Create input file
$SOX -r $1 -n -twavpcm $IN synth 20 sin 0:`expr $1 / 2` gain -2
for src in \
- "polyphase" \
- "rabbit -c0" \
- "rabbit -c1" \
- "rate -h" \
- "rate -v" \
- "resample -q" \
- "resample -ql" \
- "ssrc" \
- "ssrc_hp" \
+ "rate -vsL" \
; do
echo $src
case "$src" in