src/noisered.c

   1 /* noisered - Noise Reduction Effect.
   2  *
   3  * Written by Ian Turner (vectro@vectro.org)
   4  *
   5  * Copyright 1999 Ian Turner
   6  * This source code is freely redistributable and may be used for
   7  * any purpose.  This copyright notice must be maintained.
   8  * Authors are not responsible for the consequences of using this software.
   9  */
  10
  11 #include "noisered.h"
  12
  13 #include <stdlib.h>
  14 #include <errno.h>
  15 #include <string.h>
  16 #include <assert.h>
  17
  18 typedef struct {
  19     float *window;
  20     float *lastwindow;
  21     float *noisegate;
  22     float *smoothing;
  23 } chandata_t;
  24
  25 /* Holds profile information */
  26 typedef struct {
  27     char* profile_filename;
  28     float threshold;
  29
  30     chandata_t *chandata;
  31     size_t bufdata;
  32 } priv_t;
  33
  34 static void FFT(unsigned NumSamples,
  35          int InverseTransform,
  36          const float *RealIn, float *ImagIn, float *RealOut, float *ImagOut)
  37 {
  38   unsigned i;
  39   double * work = malloc(2 * NumSamples * sizeof(*work));
  40   for (i = 0; i < 2 * NumSamples; i += 2) {
  41     work[i] = RealIn[i >> 1];
  42     work[i + 1] = ImagIn? ImagIn[i >> 1] : 0;
  43   }
  44   lsx_safe_cdft(2 * (int)NumSamples, InverseTransform? -1 : 1, work);
  45   if (InverseTransform) for (i = 0; i < 2 * NumSamples; i += 2) {
  46     RealOut[i >> 1] = work[i] / NumSamples;
  47     ImagOut[i >> 1] = work[i + 1] / NumSamples;
  48   }
  49   else for (i = 0; i < 2 * NumSamples; i += 2) {
  50     RealOut[i >> 1] = work[i];
  51     ImagOut[i >> 1] = work[i + 1];
  52   }
  53   free(work);
  54 }
  55
  56 /*
  57  * Get the options. Default file is stdin (if the audio
  58  * input file isn't coming from there, of course!)
  59  */
  60 static int sox_noisered_getopts(sox_effect_t * effp, int argc, char **argv)
  61 {
  62   priv_t * p = (priv_t *) effp->priv;
  63   --argc, ++argv;
  64
  65   if (argc > 0) {
  66     p->profile_filename = argv[0];
  67     ++argv;
  68     --argc;
  69   }
  70
  71   p->threshold = 0.5;
  72   do {     /* break-able block */
  73     NUMERIC_PARAMETER(threshold, 0, 1);
  74   } while (0);
  75
  76   return argc? lsx_usage(effp) : SOX_SUCCESS;
  77 }
  78
  79 /*
  80  * Prepare processing.
  81  * Do all initializations.
  82  */
  83 static int sox_noisered_start(sox_effect_t * effp)
  84 {
  85     priv_t * data = (priv_t *) effp->priv;
  86     size_t fchannels = 0;
  87     size_t channels = effp->in_signal.channels;
  88     size_t i;
  89     FILE * ifp = lsx_open_input_file(effp, data->profile_filename, sox_false);
  90
  91     if (!ifp)
  92       return SOX_EOF;
  93
  94     data->chandata = lsx_calloc(channels, sizeof(*(data->chandata)));
  95     data->bufdata = 0;
  96     for (i = 0; i < channels; i ++) {
  97         data->chandata[i].noisegate = lsx_calloc(FREQCOUNT, sizeof(float));
  98         data->chandata[i].smoothing = lsx_calloc(FREQCOUNT, sizeof(float));
  99         data->chandata[i].lastwindow = NULL;
 100     }
 101     while (1) {
 102         unsigned long i1_ul;
 103         size_t i1;
 104         float f1;
 105         if (2 != fscanf(ifp, " Channel %lu: %f", &i1_ul, &f1))
 106             break;
 107         i1 = i1_ul;
 108         if (i1 != fchannels) {
 109             lsx_fail("noisered: Got channel %lu, expected channel %lu.",
 110                     (unsigned long)i1, (unsigned long)fchannels);
 111             return SOX_EOF;
 112         }
 113
 114         data->chandata[fchannels].noisegate[0] = f1;
 115         for (i = 1; i < FREQCOUNT; i ++) {
 116             if (1 != fscanf(ifp, ", %f", &f1)) {
 117                 lsx_fail("noisered: Not enough data for channel %lu "
 118                         "(expected %d, got %lu)", (unsigned long)fchannels, FREQCOUNT, (unsigned long)i);
 119                 return SOX_EOF;
 120             }
 121             data->chandata[fchannels].noisegate[i] = f1;
 122         }
 123         fchannels ++;
 124     }
 125     if (fchannels != channels) {
 126         lsx_fail("noisered: channel mismatch: %lu in input, %lu in profile.",
 127                 (unsigned long)channels, (unsigned long)fchannels);
 128         return SOX_EOF;
 129     }
 130     if (ifp != stdin)
 131       fclose(ifp);
 132
 133   effp->out_signal.length = SOX_UNKNOWN_LEN; /* TODO: calculate actual length */
 134
 135     return (SOX_SUCCESS);
 136 }
 137
 138 /* Mangle a single window. Each output sample (except the first and last
 139  * half-window) is the result of two distinct calls to this function,
 140  * due to overlapping windows. */
 141 static void reduce_noise(chandata_t* chan, float* window, double level)
 142 {
 143     float *inr, *ini, *outr, *outi, *power;
 144     float *smoothing = chan->smoothing;
 145     int i;
 146
 147     inr = lsx_calloc(WINDOWSIZE * 5, sizeof(float));
 148     ini = inr + WINDOWSIZE;
 149     outr = ini + WINDOWSIZE;
 150     outi = outr + WINDOWSIZE;
 151     power = outi + WINDOWSIZE;
 152
 153     for (i = 0; i < FREQCOUNT; i ++)
 154         assert(smoothing[i] >= 0 && smoothing[i] <= 1);
 155
 156     memcpy(inr, window, WINDOWSIZE*sizeof(float));
 157
 158     FFT(WINDOWSIZE, 0, inr, NULL, outr, outi);
 159
 160     memcpy(inr, window, WINDOWSIZE*sizeof(float));
 161     lsx_apply_hann_f(inr, WINDOWSIZE);
 162     lsx_power_spectrum_f(WINDOWSIZE, inr, power);
 163
 164     for (i = 0; i < FREQCOUNT; i ++) {
 165         float smooth;
 166         float plog;
 167         plog = log(power[i]);
 168         if (power[i] != 0 && plog < chan->noisegate[i] + level*8.0)
 169             smooth = 0.0;
 170         else
 171             smooth = 1.0;
 172
 173         smoothing[i] = smooth * 0.5 + smoothing[i] * 0.5;
 174     }
 175
 176     /* Audacity says this code will eliminate tinkle bells.
 177      * I have no idea what that means. */
 178     for (i = 2; i < FREQCOUNT - 2; i ++) {
 179         if (smoothing[i]>=0.5 &&
 180             smoothing[i]<=0.55 &&
 181             smoothing[i-1]<0.1 &&
 182             smoothing[i-2]<0.1 &&
 183             smoothing[i+1]<0.1 &&
 184             smoothing[i+2]<0.1)
 185             smoothing[i] = 0.0;
 186     }
 187
 188     outr[0] *= smoothing[0];
 189     outi[0] *= smoothing[0];
 190     outr[FREQCOUNT-1] *= smoothing[FREQCOUNT-1];
 191     outi[FREQCOUNT-1] *= smoothing[FREQCOUNT-1];
 192
 193     for (i = 1; i < FREQCOUNT-1; i ++) {
 194         int j = WINDOWSIZE - i;
 195         float smooth = smoothing[i];
 196
 197         outr[i] *= smooth;
 198         outi[i] *= smooth;
 199         outr[j] *= smooth;
 200         outi[j] *= smooth;
 201     }
 202
 203     FFT(WINDOWSIZE, 1, outr, outi, inr, ini);
 204     lsx_apply_hann_f(inr, WINDOWSIZE);
 205
 206     memcpy(window, inr, WINDOWSIZE*sizeof(float));
 207
 208     for (i = 0; i < FREQCOUNT; i ++)
 209         assert(smoothing[i] >= 0 && smoothing[i] <= 1);
 210
 211     free(inr);
 212 }
 213
 214 /* Do window management once we have a complete window, including mangling
 215  * the current window. */
 216 static int process_window(sox_effect_t * effp, priv_t * data, unsigned chan_num, unsigned num_chans,
 217                           sox_sample_t *obuf, unsigned len) {
 218     int j;
 219     float* nextwindow;
 220     int use = min(len, WINDOWSIZE)-min(len,(WINDOWSIZE/2));
 221     chandata_t *chan = &(data->chandata[chan_num]);
 222     int first = (chan->lastwindow == NULL);
 223     SOX_SAMPLE_LOCALS;
 224
 225     if ((nextwindow = lsx_calloc(WINDOWSIZE, sizeof(float))) == NULL)
 226         return SOX_EOF;
 227
 228     memcpy(nextwindow, chan->window+WINDOWSIZE/2,
 229            sizeof(float)*(WINDOWSIZE/2));
 230
 231     reduce_noise(chan, chan->window, data->threshold);
 232     if (!first) {
 233         for (j = 0; j < use; j ++) {
 234             float s = chan->window[j] + chan->lastwindow[WINDOWSIZE/2 + j];
 235             obuf[chan_num + num_chans * j] =
 236                 SOX_FLOAT_32BIT_TO_SAMPLE(s, effp->clips);
 237         }
 238         free(chan->lastwindow);
 239     } else {
 240         for (j = 0; j < use; j ++) {
 241             assert(chan->window[j] >= -1 && chan->window[j] <= 1);
 242             obuf[chan_num + num_chans * j] =
 243                 SOX_FLOAT_32BIT_TO_SAMPLE(chan->window[j], effp->clips);
 244         }
 245     }
 246     chan->lastwindow = chan->window;
 247     chan->window = nextwindow;
 248
 249     return use;
 250 }
 251
 252 /*
 253  * Read in windows, and call process_window once we get a whole one.
 254  */
 255 static int sox_noisered_flow(sox_effect_t * effp, const sox_sample_t *ibuf, sox_sample_t *obuf,
 256                     size_t *isamp, size_t *osamp)
 257 {
 258     priv_t * data = (priv_t *) effp->priv;
 259     size_t samp = min(*isamp, *osamp);
 260     size_t tracks = effp->in_signal.channels;
 261     size_t track_samples = samp / tracks;
 262     size_t ncopy = min(track_samples, WINDOWSIZE-data->bufdata);
 263     size_t whole_window = (ncopy + data->bufdata == WINDOWSIZE);
 264     int oldbuf = data->bufdata;
 265     size_t i;
 266
 267     /* FIXME: Make this automatic for all effects */
 268     assert(effp->in_signal.channels == effp->out_signal.channels);
 269
 270     if (whole_window)
 271         data->bufdata = WINDOWSIZE/2;
 272     else
 273         data->bufdata += ncopy;
 274
 275     /* Reduce noise on every channel. */
 276     for (i = 0; i < tracks; i ++) {
 277         SOX_SAMPLE_LOCALS;
 278         chandata_t* chan = &(data->chandata[i]);
 279         size_t j;
 280
 281         if (chan->window == NULL)
 282             chan->window = lsx_calloc(WINDOWSIZE, sizeof(float));
 283
 284         for (j = 0; j < ncopy; j ++)
 285             chan->window[oldbuf + j] =
 286                 SOX_SAMPLE_TO_FLOAT_32BIT(ibuf[i + tracks * j], effp->clips);
 287
 288         if (!whole_window)
 289             continue;
 290         else
 291             process_window(effp, data, (unsigned) i, (unsigned) tracks, obuf, (unsigned) (oldbuf + ncopy));
 292     }
 293
 294     *isamp = tracks*ncopy;
 295     if (whole_window)
 296         *osamp = tracks*(WINDOWSIZE/2);
 297     else
 298         *osamp = 0;
 299
 300     return SOX_SUCCESS;
 301 }
 302
 303 /*
 304  * We have up to half a window left to dump.
 305  */
 306
 307 static int sox_noisered_drain(sox_effect_t * effp, sox_sample_t *obuf, size_t *osamp)
 308 {
 309     priv_t * data = (priv_t *)effp->priv;
 310     unsigned i;
 311     unsigned tracks = effp->in_signal.channels;
 312     for (i = 0; i < tracks; i ++)
 313         *osamp = process_window(effp, data, i, tracks, obuf, (unsigned) data->bufdata);
 314
 315     /* FIXME: This is very picky.  osamp needs to be big enough to get all
 316      * remaining data or it will be discarded.
 317      */
 318     return (SOX_EOF);
 319 }
 320
 321 /*
 322  * Clean up.
 323  */
 324 static int sox_noisered_stop(sox_effect_t * effp)
 325 {
 326     priv_t * data = (priv_t *) effp->priv;
 327     size_t i;
 328
 329     for (i = 0; i < effp->in_signal.channels; i ++) {
 330         chandata_t* chan = &(data->chandata[i]);
 331         free(chan->lastwindow);
 332         free(chan->window);
 333         free(chan->smoothing);
 334         free(chan->noisegate);
 335     }
 336
 337     free(data->chandata);
 338
 339     return (SOX_SUCCESS);
 340 }
 341
 342 static sox_effect_handler_t sox_noisered_effect = {
 343   "noisered",
 344   "[profile-file [amount]]",
 345   SOX_EFF_MCHAN|SOX_EFF_LENGTH,
 346   sox_noisered_getopts,
 347   sox_noisered_start,
 348   sox_noisered_flow,
 349   sox_noisered_drain,
 350   sox_noisered_stop,
 351   NULL, sizeof(priv_t)
 352 };
 353
 354 const sox_effect_handler_t *lsx_noisered_effect_fn(void)
 355 {
 356     return &sox_noisered_effect;
 357 }