ref: 105e2bc5e35b1b747b8a93da398ae34cac0a5dde
parent: 791cff8c4282427ea3e8c527af4ce40bf6056b7c
author: cbagwell <cbagwell>
date: Mon Dec 20 16:13:02 EST 2004
Adding new noise removal effect
--- /dev/null
+++ b/src/FFT.c
@@ -1,0 +1,401 @@
+/*
+ *
+ * FFT.c
+ *
+ * Based on FFT.cpp from Audacity, which contained the following text:
+ *
+ * This file contains a few FFT routines, including a real-FFT
+ * routine that is almost twice as fast as a normal complex FFT,
+ * and a power spectrum routine when you know you don't care
+ * about phase information.
+ *
+ * Some of this code was based on a free implementation of an FFT
+ * by Don Cross, available on the web at:
+ *
+ * http://www.intersrv.com/~dcross/fft.html
+ *
+ * The basic algorithm for his code was based on Numerican Recipes
+ * in Fortran. I optimized his code further by reducing array
+ * accesses, caching the bit reversal table, and eliminating
+ * float-to-double conversions, and I added the routines to
+ * calculate a real FFT and a real power spectrum.
+ *
+ * This file is now part of SoX, and is copyright Ian Turner and others.
+ *
+ * SoX is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * Foobar 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 General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with SoX; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+
+#include "FFT.h"
+
+int **gFFTBitTable = NULL;
+const int MaxFastBits = 16;
+
+/* Declare Static functions */
+static int IsPowerOfTwo(int x);
+static int NumberOfBitsNeeded(int PowerOfTwo);
+static int ReverseBits(int index, int NumBits);
+static void InitFFT();
+
+int IsPowerOfTwo(int x)
+{+ if (x < 2)
+ return 0;
+
+ return !(x & (x-1)); /* Thanks to 'byang' for this cute trick! */
+}
+
+int NumberOfBitsNeeded(int PowerOfTwo)
+{+ int i;
+
+ if (PowerOfTwo < 2) {+ fprintf(stderr, "Error: FFT called with size %d\n", PowerOfTwo);
+ exit(1);
+ }
+
+ for (i = 0;; i++)
+ if (PowerOfTwo & (1 << i))
+ return i;
+}
+
+int ReverseBits(int index, int NumBits)
+{+ int i, rev;
+
+ for (i = rev = 0; i < NumBits; i++) {+ rev = (rev << 1) | (index & 1);
+ index >>= 1;
+ }
+
+ return rev;
+}
+
+
+/* This function allocates about 250k (actually (2**16)-2 ints) which is never
+ * freed, to use as a lookup table for bit-reversal. The good news is that
+ * we bascially need this until the very end, so the fact that it's not freed
+ * is OK. */
+void InitFFT()
+{+ int len, b;
+
+ gFFTBitTable = (int**)calloc(MaxFastBits, sizeof(*gFFTBitTable));
+
+ len = 2;
+ for (b = 1; b <= MaxFastBits; b++) {+ int i;
+
+ gFFTBitTable[b - 1] = (int*)calloc(len, sizeof(**gFFTBitTable));
+
+ for (i = 0; i < len; i++) {+ gFFTBitTable[b - 1][i] = ReverseBits(i, b);
+ }
+
+ len <<= 1;
+ }
+}
+
+inline int FastReverseBits(int i, int NumBits)
+{+ if (NumBits <= MaxFastBits)
+ return gFFTBitTable[NumBits - 1][i];
+ else
+ return ReverseBits(i, NumBits);
+}
+
+/*
+ * Complex Fast Fourier Transform
+ */
+
+void FFT(int NumSamples,
+ int InverseTransform,
+ float *RealIn, float *ImagIn, float *RealOut, float *ImagOut)
+{+ int NumBits; /* Number of bits needed to store indices */
+ int i, j, k, n;
+ int BlockSize, BlockEnd;
+
+ double angle_numerator = 2.0 * M_PI;
+ float tr, ti; /* temp real, temp imaginary */
+
+ if (!IsPowerOfTwo(NumSamples)) {+ fprintf(stderr, "%d is not a power of two\n", NumSamples);
+ exit(1);
+ }
+
+ if (!gFFTBitTable)
+ InitFFT();
+
+ if (InverseTransform)
+ angle_numerator = -angle_numerator;
+
+ NumBits = NumberOfBitsNeeded(NumSamples);
+
+ /*
+ ** Do simultaneous data copy and bit-reversal ordering into outputs...
+ */
+
+ for (i = 0; i < NumSamples; i++) {+ j = FastReverseBits(i, NumBits);
+ RealOut[j] = RealIn[i];
+ ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i];
+ }
+
+ /*
+ ** Do the FFT itself...
+ */
+
+ BlockEnd = 1;
+ for (BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) {+
+ double delta_angle = angle_numerator / (double) BlockSize;
+
+ float sm2 = sin(-2 * delta_angle);
+ float sm1 = sin(-delta_angle);
+ float cm2 = cos(-2 * delta_angle);
+ float cm1 = cos(-delta_angle);
+ float w = 2 * cm1;
+ float ar0, ar1, ar2, ai0, ai1, ai2;
+
+ for (i = 0; i < NumSamples; i += BlockSize) {+ ar2 = cm2;
+ ar1 = cm1;
+
+ ai2 = sm2;
+ ai1 = sm1;
+
+ for (j = i, n = 0; n < BlockEnd; j++, n++) {+ ar0 = w * ar1 - ar2;
+ ar2 = ar1;
+ ar1 = ar0;
+
+ ai0 = w * ai1 - ai2;
+ ai2 = ai1;
+ ai1 = ai0;
+
+ k = j + BlockEnd;
+ tr = ar0 * RealOut[k] - ai0 * ImagOut[k];
+ ti = ar0 * ImagOut[k] + ai0 * RealOut[k];
+
+ RealOut[k] = RealOut[j] - tr;
+ ImagOut[k] = ImagOut[j] - ti;
+
+ RealOut[j] += tr;
+ ImagOut[j] += ti;
+ }
+ }
+
+ BlockEnd = BlockSize;
+ }
+
+ /*
+ ** Need to normalize if inverse transform...
+ */
+
+ if (InverseTransform) {+ float denom = (float) NumSamples;
+
+ for (i = 0; i < NumSamples; i++) {+ RealOut[i] /= denom;
+ ImagOut[i] /= denom;
+ }
+ }
+}
+
+/*
+ * Real Fast Fourier Transform
+ *
+ * This function was based on the code in Numerical Recipes for C.
+ * In Num. Rec., the inner loop is based on a single 1-based array
+ * of interleaved real and imaginary numbers. Because we have two
+ * separate zero-based arrays, our indices are quite different.
+ * Here is the correspondence between Num. Rec. indices and our indices:
+ *
+ * i1 <-> real[i]
+ * i2 <-> imag[i]
+ * i3 <-> real[n/2-i]
+ * i4 <-> imag[n/2-i]
+ */
+
+void RealFFT(int NumSamples, float *RealIn, float *RealOut, float *ImagOut)
+{+ int Half = NumSamples / 2;
+ int i;
+
+ float theta = M_PI / Half;
+ float wtemp = (float) sin(0.5 * theta);
+ float wpr = -2.0 * wtemp * wtemp;
+ float wpi = (float) sin(theta);
+ float wr = 1.0 + wpr;
+ float wi = wpi;
+
+ int i3;
+
+ float h1r, h1i, h2r, h2i;
+
+
+ float *tmpReal = (float*)calloc(Half, sizeof(float));
+ float *tmpImag = (float*)calloc(Half, sizeof(float));
+
+ for (i = 0; i < Half; i++) {+ tmpReal[i] = RealIn[2 * i];
+ tmpImag[i] = RealIn[2 * i + 1];
+ }
+
+ FFT(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);
+
+
+ for (i = 1; i < Half / 2; i++) {+
+ i3 = Half - i;
+
+ h1r = 0.5 * (RealOut[i] + RealOut[i3]);
+ h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
+ h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
+ h2i = -0.5 * (RealOut[i] - RealOut[i3]);
+
+ RealOut[i] = h1r + wr * h2r - wi * h2i;
+ ImagOut[i] = h1i + wr * h2i + wi * h2r;
+ RealOut[i3] = h1r - wr * h2r + wi * h2i;
+ ImagOut[i3] = -h1i + wr * h2i + wi * h2r;
+
+ wr = (wtemp = wr) * wpr - wi * wpi + wr;
+ wi = wi * wpr + wtemp * wpi + wi;
+ }
+
+ RealOut[0] = (h1r = RealOut[0]) + ImagOut[0];
+ ImagOut[0] = h1r - ImagOut[0];
+
+ free(tmpReal);
+ free(tmpImag);
+}
+
+/*
+ * PowerSpectrum
+ *
+ * This function computes the same as RealFFT, above, but
+ * adds the squares of the real and imaginary part of each
+ * coefficient, extracting the power and throwing away the
+ * phase.
+ *
+ * For speed, it does not call RealFFT, but duplicates some
+ * of its code.
+ */
+
+void PowerSpectrum(int NumSamples, float *In, float *Out)
+{+ int Half, i, i3;
+ float theta, wtemp, wpr, wpi, wr, wi;
+ float h1r, h1i, h2r, h2i, rt, it;
+ float *tmpReal;
+ float *tmpImag;
+ float *RealOut;
+ float *ImagOut;
+
+ Half = NumSamples / 2;
+
+ theta = M_PI / Half;
+
+ tmpReal = (float*)calloc(Half, sizeof(float));
+ tmpImag = (float*)calloc(Half, sizeof(float));
+ RealOut = (float*)calloc(Half, sizeof(float));
+ ImagOut = (float*)calloc(Half, sizeof(float));
+
+ for (i = 0; i < Half; i++) {+ tmpReal[i] = In[2 * i];
+ tmpImag[i] = In[2 * i + 1];
+ }
+
+ FFT(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);
+
+ wtemp = (float) sin(0.5 * theta);
+
+ wpr = -2.0 * wtemp * wtemp;
+ wpi = (float) sin(theta);
+ wr = 1.0 + wpr;
+ wi = wpi;
+
+ for (i = 1; i < Half / 2; i++) {+
+ i3 = Half - i;
+
+ h1r = 0.5 * (RealOut[i] + RealOut[i3]);
+ h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
+ h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
+ h2i = -0.5 * (RealOut[i] - RealOut[i3]);
+
+ rt = h1r + wr * h2r - wi * h2i;
+ it = h1i + wr * h2i + wi * h2r;
+
+ Out[i] = rt * rt + it * it;
+
+ rt = h1r - wr * h2r + wi * h2i;
+ it = -h1i + wr * h2i + wi * h2r;
+
+ Out[i3] = rt * rt + it * it;
+
+ wr = (wtemp = wr) * wpr - wi * wpi + wr;
+ wi = wi * wpr + wtemp * wpi + wi;
+ }
+
+ rt = (h1r = RealOut[0]) + ImagOut[0];
+ it = h1r - ImagOut[0];
+ Out[0] = rt * rt + it * it;
+
+ rt = RealOut[Half / 2];
+ it = ImagOut[Half / 2];
+ Out[Half / 2] = rt * rt + it * it;
+
+ free(tmpReal);
+ free(tmpImag);
+ free(RealOut);
+ free(ImagOut);
+}
+
+/*
+ * Windowing Functions
+ */
+
+void WindowFunc(windowfunc_t whichFunction, int NumSamples, float *in)
+{+ int i;
+
+ switch (whichFunction) {+ case BARTLETT:
+ for (i = 0; i < NumSamples / 2; i++) {+ in[i] *= (i / (float) (NumSamples / 2));
+ in[i + (NumSamples / 2)] *=
+ (1.0 - (i / (float) (NumSamples / 2)));
+ }
+ break;
+
+ case HAMMING:
+ for (i = 0; i < NumSamples; i++)
+ in[i] *= 0.54 - 0.46 * cos(2 * M_PI * i / (NumSamples - 1));
+ break;
+
+ case HANNING:
+ for (i = 0; i < NumSamples; i++)
+ in[i] *= 0.50 - 0.50 * cos(2 * M_PI * i / (NumSamples - 1));
+ break;
+ case RECTANGULAR:
+ break;
+ }
+}
--- /dev/null
+++ b/src/FFT.h
@@ -1,0 +1,96 @@
+/*
+ *
+ * FFT.h
+ *
+ * Based on FFT.h from Audacity, which contained the following text:
+ *
+ * This file contains a few FFT routines, including a real-FFT
+ * routine that is almost twice as fast as a normal complex FFT,
+ * and a power spectrum routine when you know you don't care
+ * about phase information.
+ *
+ * Some of this code was based on a free implementation of an FFT
+ * by Don Cross, available on the web at:
+ *
+ * http://www.intersrv.com/~dcross/fft.html
+ *
+ * The basic algorithm for his code was based on Numerican Recipes
+ * in Fortran. I optimized his code further by reducing array
+ * accesses, caching the bit reversal table, and eliminating
+ * float-to-double conversions, and I added the routines to
+ * calculate a real FFT and a real power spectrum.
+ *
+ * Note: all of these routines use single-precision floats.
+ * I have found that in practice, floats work well until you
+ * get above 8192 samples. If you need to do a larger FFT,
+ * you need to use doubles.
+ *
+ * This file is now part of SoX, and is copyright Ian Turner and others.
+ *
+ * SoX is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * Foobar 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 General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with SoX; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ */
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846 /* pi */
+#endif
+
+#ifdef __cplusplus
+extern "C" {+#endif
+
+/*
+ * This is the function you will use the most often.
+ * Given an array of floats, this will compute the power
+ * spectrum by doing a Real FFT and then computing the
+ * sum of the squares of the real and imaginary parts.
+ * Note that the output array is half the length of the
+ * input array, and that NumSamples must be a power of two.
+ */
+
+void PowerSpectrum(int NumSamples, float *In, float *Out);
+
+/*
+ * Computes an FFT when the input data is real but you still
+ * want complex data as output. The output arrays are half
+ * the length of the input, and NumSamples must be a power of
+ * two.
+ */
+
+void RealFFT(int NumSamples,
+ float *RealIn, float *RealOut, float *ImagOut);
+
+/*
+ * Computes a FFT of complex input and returns complex output.
+ * Currently this is the only function here that supports the
+ * inverse transform as well.
+ */
+
+void FFT(int NumSamples,
+ int InverseTransform,
+ float *RealIn, float *ImagIn, float *RealOut, float *ImagOut);
+
+/*
+ * Applies a windowing function to the data in place
+ */
+typedef enum {RECTANGULAR = 0, /* no window */+ BARTLETT = 1, /* triangular */
+ HAMMING = 2,
+ HANNING = 3} windowfunc_t;
+
+void WindowFunc(windowfunc_t whichFunction, int NumSamples, float *data);
+
+#ifdef __cplusplus
+}
+#endif
--- /dev/null
+++ b/src/noiseprof.c
@@ -1,0 +1,212 @@
+/*
+ * noiseprof - Noise Profiling Effect.
+ *
+ * Written by Ian Turner (vectro@vectro.org)
+ *
+ * Copyright 1999 Ian Turner and others
+ * This file is part of SoX.
+
+ * SoX 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 of the License, or
+ * (at your option) any later version.
+ *
+ */
+
+#include "noisered.h"
+
+#include <assert.h>
+#include <string.h>
+#include <errno.h>
+
+typedef struct chandata {+ float *sum;
+ int *profilecount;
+
+ float *window;
+} chandata_t;
+
+typedef struct profdata {+ char* output_filename;
+ FILE* output_file;
+
+ chandata_t *chandata;
+ st_size_t bufdata;
+} * profdata_t;
+
+/*
+ * Get the filename, if any. We don't open it until st_noiseprof_start.
+ */
+int st_noiseprof_getopts(eff_t effp, int n, char **argv)
+{+ profdata_t data = (profdata_t) effp->priv;
+
+ if (n == 1) {+ data->output_filename = argv[0];
+ } else if (n > 1) {+ st_fail("Usage: noiseprof [filename]");+ return (ST_EOF);
+ }
+
+ return (ST_SUCCESS);
+}
+
+/*
+ * Prepare processing.
+ * Do all initializations.
+ */
+int st_noiseprof_start(eff_t effp)
+{+ profdata_t data = (profdata_t) effp->priv;
+ int channels = effp->ininfo.channels;
+ int i;
+
+ if (data->output_filename != NULL) {+ data->output_file = fopen(data->output_filename, "w");
+ if (data->output_file == NULL) {+ st_fail("Couldn't open output file %s: %s",+ data->output_filename, strerror(errno));
+ }
+ } else {+ /* FIXME: We should detect output to stdout, and redirect to stderr. */
+ data->output_file = stderr;
+ }
+
+ data->chandata = (chandata_t*)calloc(channels, sizeof(*(data->chandata)));
+ for (i = 0; i < channels; i ++) {+ data->chandata[i].sum =
+ (float*)calloc(FREQCOUNT, sizeof(float));
+ data->chandata[i].profilecount =
+ (int*)calloc(FREQCOUNT, sizeof(int));
+ data->chandata[i].window =
+ (float*)calloc(WINDOWSIZE, sizeof(float));
+ }
+ data->bufdata = 0;
+ return (ST_SUCCESS);
+}
+
+/* Collect statistics from the complete window on channel chan. */
+static void collect_data(profdata_t data, chandata_t* chan) {+ float *out = (float*)calloc(FREQCOUNT, sizeof(float));
+
+ int i;
+
+ PowerSpectrum(WINDOWSIZE, chan->window, out);
+
+ for (i = 0; i < FREQCOUNT; i ++) {+ float value = log(out[i]);
+ if (finite(value)) {+ chan->sum[i] += value;
+ chan->profilecount[i] ++;
+ }
+ }
+
+ free(out);
+}
+
+/*
+ * Grab what we can from ibuf, and process if we have a whole window.
+ */
+int st_noiseprof_flow(eff_t effp, st_sample_t *ibuf, st_sample_t *obuf,
+ st_size_t *isamp, st_size_t *osamp)
+{+ profdata_t data = (profdata_t) effp->priv;
+ int samp = min(*isamp, *osamp);
+ int tracks = effp->ininfo.channels;
+ int track_samples = samp / tracks;
+ int ncopy = 0;
+ int i;
+
+ assert(effp->ininfo.channels == effp->outinfo.channels);
+
+ /* How many samples per track to analyze? */
+ ncopy = min(track_samples, WINDOWSIZE-data->bufdata);
+
+ /* Collect data for every channel. */
+ for (i = 0; i < tracks; i ++) {+ chandata_t* chan = &(data->chandata[i]);
+ int j;
+ for (j = 0; j < ncopy; j ++) {+ chan->window[j+data->bufdata] =
+ ST_SAMPLE_TO_FLOAT_DWORD(ibuf[i+j*tracks]);
+ }
+ if (ncopy + data->bufdata == WINDOWSIZE) {+ collect_data(data, chan);
+ }
+ }
+
+ data->bufdata += ncopy;
+ assert(data->bufdata <= WINDOWSIZE);
+ if (data->bufdata == WINDOWSIZE)
+ data->bufdata = 0;
+
+ memcpy(obuf, ibuf, ncopy*tracks);
+ *isamp = *osamp = ncopy*tracks;
+
+ return (ST_SUCCESS);
+}
+
+/*
+ * Finish off the last window.
+ */
+
+int st_noiseprof_drain(eff_t effp, st_sample_t *obuf, st_size_t *osamp)
+{+ profdata_t data = (profdata_t) effp->priv;
+ int tracks = effp->ininfo.channels;
+ int i;
+
+ *osamp = 0;
+
+ if (data->bufdata == 0) {+ return ST_SUCCESS;
+ }
+
+ for (i = 0; i < tracks; i ++) {+ int j;
+ for (j = data->bufdata+1; j < WINDOWSIZE; j ++) {+ data->chandata[i].window[j] = 0;
+ }
+ collect_data(data, &(data->chandata[i]));
+ }
+
+ return (ST_SUCCESS);
+}
+
+/*
+ * Print profile and clean up.
+ */
+int st_noiseprof_stop(eff_t effp)
+{+ profdata_t data = (profdata_t) effp->priv;
+ int i;
+
+ for (i = 0; i < effp->ininfo.channels; i ++) {+ chandata_t* chan = &(data->chandata[i]);
+ int j;
+
+ fprintf(data->output_file, "Channel %d: ", i);
+
+ for (j = 0; j < FREQCOUNT; j ++) {+ fprintf(data->output_file, "%s%f", j == 0 ? "" : ", ",
+ chan->sum[j] / chan->profilecount[j]);
+ }
+ fprintf(data->output_file, "\n");
+
+ free(chan->sum);
+ free(chan->profilecount);
+ }
+
+ free(data->chandata);
+
+ if (data->output_file != stderr && data->output_file != stdout) {+ fclose(data->output_file);
+ }
+
+ return (ST_SUCCESS);
+}
+
+/* For Emacs:
+ Local Variables:
+ c-basic-offset: 4
+*/
--- /dev/null
+++ b/src/noisered.c
@@ -1,0 +1,340 @@
+/*
+ * noiseprof - Noise Profiling Effect.
+ *
+ * Written by Ian Turner (vectro@vectro.org)
+ *
+ * Copyright 1999 Ian Turner
+ * This source code is freely redistributable and may be used for
+ * any purpose. This copyright notice must be maintained.
+ * Authors are not responsible for the consequences of using this software.
+ */
+
+#include "noisered.h"
+
+#include <stdlib.h>
+#include <errno.h>
+#include <string.h>
+#include <assert.h>
+
+typedef struct chandata {+ float *window;
+ float *lastwindow;
+ float *noisegate;
+ float *smoothing;
+} chandata_t;
+
+/* Holds profile information */
+typedef struct reddata {+ char* profile_filename;
+ float threshold;
+
+ chandata_t *chandata;
+ st_size_t bufdata;
+} * reddata_t;
+
+/*
+ * Get the options. Filename is mandatory, though a reasonable default would
+ * be stdin (if the input file isn't coming from there, of course!)
+ */
+int st_noisered_getopts(eff_t effp, int n, char **argv)
+{+ reddata_t data = (reddata_t) effp->priv;
+
+ if (n > 2 || n < 1) {+ st_fail("Usage: noiseprof profile-file [threshold]");+ return (ST_EOF);
+ }
+ data->threshold = 0.5;
+ data->profile_filename = argv[0];
+ if (n == 2)
+ {+ data->threshold = atof(argv[1]);
+
+ if (data->threshold > 1)
+ {+ data->threshold = 1;
+ } else if (data->threshold < 0)
+ {+ data->threshold = 0;
+ }
+ }
+ return (ST_SUCCESS);
+}
+
+/*
+ * Prepare processing.
+ * Do all initializations.
+ */
+int st_noisered_start(eff_t effp)
+{+ reddata_t data = (reddata_t) effp->priv;
+ int fchannels = 0;
+ int channels = effp->ininfo.channels;
+ int i;
+ FILE* ifd;
+
+ data->chandata = (chandata_t*)calloc(channels, sizeof(*(data->chandata)));
+ for (i = 0; i < channels; i ++) {+ data->chandata[i].noisegate = (float*)calloc(FREQCOUNT, sizeof(float));
+ data->chandata[i].smoothing = (float*)calloc(FREQCOUNT, sizeof(float));
+ data->chandata[i].lastwindow = NULL;
+ }
+ data->bufdata = 0;
+
+ /* Here we actually open the input file. */
+ ifd = fopen(data->profile_filename, "r");
+ if (ifd == NULL) {+ st_fail("Couldn't open profile file %s: %s",+ data->profile_filename, strerror(errno));
+ return ST_EOF;
+ }
+
+ while (1) {+ int i1;
+ float f1;
+ if (2 != fscanf(ifd, " Channel %d: %f", &i1, &f1))
+ break;
+ if (i1 != fchannels) {+ st_fail("noisered: Got channel %d, expected channel %d.",+ i1, fchannels);
+ return ST_EOF;
+ }
+
+ data->chandata[fchannels].noisegate[0] = f1;
+ for (i = 1; i < FREQCOUNT; i ++) {+ if (1 != fscanf(ifd, ", %f", &f1)) {+ st_fail("noisered: Not enough datums for channel %d "+ "(expected %d, got %d)", fchannels, FREQCOUNT, i);
+ return ST_EOF;
+ }
+ data->chandata[fchannels].noisegate[i] = f1;
+ }
+ fchannels ++;
+ }
+ if (fchannels != channels) {+ st_fail("noisered: channel mismatch: %d in input, %d in profile.\n",+ channels, fchannels);
+ return ST_EOF;
+ }
+ fclose(ifd);
+
+ return (ST_SUCCESS);
+}
+
+/* Mangle a single window. Each output sample (except the first and last
+ * half-window) is the result of two distinct calls to this function,
+ * due to overlapping windows. */
+static void reduce_noise(chandata_t* chan, float* window, float level)
+{+ float *inr = (float*)calloc(WINDOWSIZE, sizeof(float));
+ float *ini = (float*)calloc(WINDOWSIZE, sizeof(float));
+ float *outr = (float*)calloc(WINDOWSIZE, sizeof(float));
+ float *outi = (float*)calloc(WINDOWSIZE, sizeof(float));
+ float *power = (float*)calloc(WINDOWSIZE, sizeof(float));
+ float *smoothing = chan->smoothing;
+ static int callnum = 0;
+
+ callnum ++;
+
+ int i;
+
+ for (i = 0; i < FREQCOUNT; i ++) {+ assert(smoothing[i] >= 0 && smoothing[i] <= 1);
+ }
+
+ memcpy(inr, window, WINDOWSIZE*sizeof(float));
+
+ FFT(WINDOWSIZE, 0, inr, NULL, outr, outi);
+
+ memcpy(inr, window, WINDOWSIZE*sizeof(float));
+ WindowFunc(HANNING, WINDOWSIZE, inr);
+ PowerSpectrum(WINDOWSIZE, inr, power);
+
+ for(i = 0; i < FREQCOUNT; i ++) {+ float smooth;
+ float plog;
+ plog = log(power[i]);
+ if (power[i] != 0 && plog < chan->noisegate[i] + level*8.0)
+ smooth = 0.0;
+ else
+ smooth = 1.0;
+
+ smoothing[i] = smooth * 0.5 + smoothing[i] * 0.5;
+ }
+
+ /* Audacity says this code will eliminate tinkle bells.
+ * I have no idea what that means. */
+ for (i = 2; i < FREQCOUNT - 2; i ++) {+ if (smoothing[i]>=0.5 &&
+ smoothing[i]<=0.55 &&
+ smoothing[i-1]<0.1 &&
+ smoothing[i-2]<0.1 &&
+ smoothing[i+1]<0.1 &&
+ smoothing[i+2]<0.1)
+ smoothing[i] = 0.0;
+ }
+
+ outr[0] *= smoothing[0];
+ outi[0] *= smoothing[0];
+ outr[FREQCOUNT-1] *= smoothing[FREQCOUNT-1];
+ outi[FREQCOUNT-1] *= smoothing[FREQCOUNT-1];
+
+ for (i = 1; i < FREQCOUNT-1; i ++) {+ int j = WINDOWSIZE - i;
+ float smooth = smoothing[i];
+
+ outr[i] *= smooth;
+ outi[i] *= smooth;
+ outr[j] *= smooth;
+ outi[j] *= smooth;
+ }
+
+ FFT(WINDOWSIZE, 1, outr, outi, inr, ini);
+ WindowFunc(HANNING, WINDOWSIZE, inr);
+
+ memcpy(window, inr, WINDOWSIZE*sizeof(float));
+
+ free(inr);
+ free(ini);
+ free(outr);
+ free(outi);
+ free(power);
+
+ for (i = 0; i < FREQCOUNT; i ++) {+ assert(smoothing[i] >= 0 && smoothing[i] <= 1);
+ }
+}
+
+/* Do window management once we have a complete window, including mangling
+ * the current window. */
+static int process_window(reddata_t data, int chan_num, int num_chans,
+ st_sample_t *obuf, int len) {+ int j;
+ float* nextwindow;
+ int use = min(len, WINDOWSIZE)-(WINDOWSIZE/2);
+ chandata_t *chan = &(data->chandata[chan_num]);
+ int first = (chan->lastwindow == NULL);
+
+ nextwindow = (float*)calloc(WINDOWSIZE, sizeof(float));
+ memcpy(nextwindow, chan->window+WINDOWSIZE/2,
+ sizeof(float)*(WINDOWSIZE/2));
+
+ reduce_noise(chan, chan->window, data->threshold);
+
+ if (!first) {+ for (j = 0; j < use; j ++) {+ float s = chan->window[j] + chan->lastwindow[WINDOWSIZE/2 + j];
+ assert(s >= -1 && s <= 1);
+ obuf[chan_num + num_chans * j] =
+ ST_FLOAT_DWORD_TO_SAMPLE(s);
+ }
+ free(chan->lastwindow);
+ } else {+ for (j = 0; j < use; j ++) {+ assert(chan->window[j] >= -1 && chan->window[j] <= 1);
+ obuf[chan_num + num_chans * j] =
+ ST_FLOAT_DWORD_TO_SAMPLE(chan->window[j]);
+ }
+ }
+ chan->lastwindow = chan->window;
+ chan->window = nextwindow;
+
+ return use;
+}
+
+/*
+ * Read in windows, and call process_window once we get a whole one.
+ */
+int st_noisered_flow(eff_t effp, st_sample_t *ibuf, st_sample_t *obuf,
+ st_size_t *isamp, st_size_t *osamp)
+{+ reddata_t data = (reddata_t) effp->priv;
+ int samp = min(*isamp, *osamp);
+ int tracks = effp->ininfo.channels;
+ int track_samples = samp / tracks;
+ int ncopy = min(track_samples, WINDOWSIZE-data->bufdata);
+ int whole_window = (ncopy + data->bufdata == WINDOWSIZE);
+ int oldbuf = data->bufdata;
+ int i;
+ assert(effp->ininfo.channels == effp->outinfo.channels);
+
+ if (whole_window) {+ data->bufdata = WINDOWSIZE/2;
+ } else {+ data->bufdata += ncopy;
+ }
+
+ /* Reduce noise on every channel. */
+ for (i = 0; i < tracks; i ++) {+ chandata_t* chan = &(data->chandata[i]);
+ int j;
+ if (chan->window == NULL) {+ chan->window = (float*)calloc(WINDOWSIZE, sizeof(float));
+ }
+
+ for (j = 0; j < ncopy; j ++) {+ chan->window[oldbuf + j] =
+ ST_SAMPLE_TO_FLOAT_DWORD(ibuf[i + tracks * j]);
+ }
+
+ if (!whole_window)
+ continue;
+ else {+ process_window(data, i, tracks, obuf, oldbuf + ncopy);
+ }
+ }
+
+ *isamp = tracks*ncopy;
+ if (whole_window) {+ *osamp = tracks*(WINDOWSIZE/2);
+ } else {+ *osamp = 0;
+ }
+
+ return (ST_SUCCESS);
+}
+
+/*
+ * We have up to half a window left to dump.
+ */
+
+int st_noisered_drain(eff_t effp, st_sample_t *obuf, st_size_t *osamp)
+{+ reddata_t data = (reddata_t)effp->priv;
+ int i;
+ int tracks = effp->ininfo.channels;
+ for (i = 0; i < tracks; i ++) {+ *osamp = process_window(data, i, tracks, obuf, data->bufdata);
+ }
+ return (ST_SUCCESS);
+}
+
+/*
+ * Clean up.
+ */
+int st_noisered_stop(eff_t effp)
+{+ reddata_t data = (reddata_t) effp->priv;
+ int i;
+
+ for (i = 0; i < effp->ininfo.channels; i ++) {+ chandata_t* chan = &(data->chandata[i]);
+ if (chan->lastwindow != NULL) {+ free(chan->lastwindow);
+ }
+ if (chan->window != NULL) {+ free(chan->window);
+ }
+ free(chan->smoothing);
+ free(chan->noisegate);
+ }
+
+ free(data->chandata);
+
+ return (ST_SUCCESS);
+}
+
+/* For Emacs:
+ Local Variables:
+ c-basic-offset: 4
+*/
--- /dev/null
+++ b/src/noisered.h
@@ -1,0 +1,26 @@
+/*
+ * noiseprof.h - Headers for Noise Profiling Effect.
+ *
+ * Written by Ian Turner (vectro@vectro.org)
+ *
+ * Copyright 1999 Ian Turner and others
+ * This file is part of SoX.
+
+ * SoX 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 of the License, or
+ * (at your option) any later version.
+ */
+
+#include "st_i.h"
+#include "FFT.h"
+
+#include <math.h>
+
+#define WINDOWSIZE 2048
+#define FREQCOUNT (WINDOWSIZE/2+1)
+
+static int min(int a, int b) {+ if (a < b) return a;
+ else return b;
+}