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);
  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     return (SOX_SUCCESS);
 134 }
 135
 136 /* Mangle a single window. Each output sample (except the first and last
 137  * half-window) is the result of two distinct calls to this function,
 138  * due to overlapping windows. */
 139 static void reduce_noise(chandata_t* chan, float* window, double level)
 140 {
 141     float *inr, *ini, *outr, *outi, *power;
 142     float *smoothing = chan->smoothing;
 143     int i;
 144
 145     inr = lsx_calloc(WINDOWSIZE * 5, sizeof(float));
 146     ini = inr + WINDOWSIZE;
 147     outr = ini + WINDOWSIZE;
 148     outi = outr + WINDOWSIZE;
 149     power = outi + WINDOWSIZE;
 150
 151     for (i = 0; i < FREQCOUNT; i ++)
 152         assert(smoothing[i] >= 0 && smoothing[i] <= 1);
 153
 154     memcpy(inr, window, WINDOWSIZE*sizeof(float));
 155
 156     FFT(WINDOWSIZE, 0, inr, NULL, outr, outi);
 157
 158     memcpy(inr, window, WINDOWSIZE*sizeof(float));
 159     lsx_apply_hann_f(inr, WINDOWSIZE);
 160     lsx_power_spectrum_f(WINDOWSIZE, inr, power);
 161
 162     for (i = 0; i < FREQCOUNT; i ++) {
 163         float smooth;
 164         float plog;
 165         plog = log(power[i]);
 166         if (power[i] != 0 && plog < chan->noisegate[i] + level*8.0)
 167             smooth = 0.0;
 168         else
 169             smooth = 1.0;
 170
 171         smoothing[i] = smooth * 0.5 + smoothing[i] * 0.5;
 172     }
 173
 174     /* Audacity says this code will eliminate tinkle bells.
 175      * I have no idea what that means. */
 176     for (i = 2; i < FREQCOUNT - 2; i ++) {
 177         if (smoothing[i]>=0.5 &&
 178             smoothing[i]<=0.55 &&
 179             smoothing[i-1]<0.1 &&
 180             smoothing[i-2]<0.1 &&
 181             smoothing[i+1]<0.1 &&
 182             smoothing[i+2]<0.1)
 183             smoothing[i] = 0.0;
 184     }
 185
 186     outr[0] *= smoothing[0];
 187     outi[0] *= smoothing[0];
 188     outr[FREQCOUNT-1] *= smoothing[FREQCOUNT-1];
 189     outi[FREQCOUNT-1] *= smoothing[FREQCOUNT-1];
 190
 191     for (i = 1; i < FREQCOUNT-1; i ++) {
 192         int j = WINDOWSIZE - i;
 193         float smooth = smoothing[i];
 194
 195         outr[i] *= smooth;
 196         outi[i] *= smooth;
 197         outr[j] *= smooth;
 198         outi[j] *= smooth;
 199     }
 200
 201     FFT(WINDOWSIZE, 1, outr, outi, inr, ini);
 202     lsx_apply_hann_f(inr, WINDOWSIZE);
 203
 204     memcpy(window, inr, WINDOWSIZE*sizeof(float));
 205
 206     for (i = 0; i < FREQCOUNT; i ++)
 207         assert(smoothing[i] >= 0 && smoothing[i] <= 1);
 208
 209     free(inr);
 210 }
 211
 212 /* Do window management once we have a complete window, including mangling
 213  * the current window. */
 214 static int process_window(sox_effect_t * effp, priv_t * data, unsigned chan_num, unsigned num_chans,
 215                           sox_sample_t *obuf, unsigned len) {
 216     int j;
 217     float* nextwindow;
 218     int use = min(len, WINDOWSIZE)-min(len,(WINDOWSIZE/2));
 219     chandata_t *chan = &(data->chandata[chan_num]);
 220     int first = (chan->lastwindow == NULL);
 221     SOX_SAMPLE_LOCALS;
 222
 223     if ((nextwindow = lsx_calloc(WINDOWSIZE, sizeof(float))) == NULL)
 224         return SOX_EOF;
 225
 226     memcpy(nextwindow, chan->window+WINDOWSIZE/2,
 227            sizeof(float)*(WINDOWSIZE/2));
 228
 229     reduce_noise(chan, chan->window, data->threshold);
 230     if (!first) {
 231         for (j = 0; j < use; j ++) {
 232             float s = chan->window[j] + chan->lastwindow[WINDOWSIZE/2 + j];
 233             obuf[chan_num + num_chans * j] =
 234                 SOX_FLOAT_32BIT_TO_SAMPLE(s, effp->clips);
 235         }
 236         free(chan->lastwindow);
 237     } else {
 238         for (j = 0; j < use; j ++) {
 239             assert(chan->window[j] >= -1 && chan->window[j] <= 1);
 240             obuf[chan_num + num_chans * j] =
 241                 SOX_FLOAT_32BIT_TO_SAMPLE(chan->window[j], effp->clips);
 242         }
 243     }
 244     chan->lastwindow = chan->window;
 245     chan->window = nextwindow;
 246
 247     return use;
 248 }
 249
 250 /*
 251  * Read in windows, and call process_window once we get a whole one.
 252  */
 253 static int sox_noisered_flow(sox_effect_t * effp, const sox_sample_t *ibuf, sox_sample_t *obuf,
 254                     size_t *isamp, size_t *osamp)
 255 {
 256     priv_t * data = (priv_t *) effp->priv;
 257     size_t samp = min(*isamp, *osamp);
 258     size_t tracks = effp->in_signal.channels;
 259     size_t track_samples = samp / tracks;
 260     size_t ncopy = min(track_samples, WINDOWSIZE-data->bufdata);
 261     size_t whole_window = (ncopy + data->bufdata == WINDOWSIZE);
 262     int oldbuf = data->bufdata;
 263     size_t i;
 264
 265     /* FIXME: Make this automatic for all effects */
 266     assert(effp->in_signal.channels == effp->out_signal.channels);
 267
 268     if (whole_window)
 269         data->bufdata = WINDOWSIZE/2;
 270     else
 271         data->bufdata += ncopy;
 272
 273     /* Reduce noise on every channel. */
 274     for (i = 0; i < tracks; i ++) {
 275         SOX_SAMPLE_LOCALS;
 276         chandata_t* chan = &(data->chandata[i]);
 277         size_t j;
 278
 279         if (chan->window == NULL)
 280             chan->window = lsx_calloc(WINDOWSIZE, sizeof(float));
 281
 282         for (j = 0; j < ncopy; j ++)
 283             chan->window[oldbuf + j] =
 284                 SOX_SAMPLE_TO_FLOAT_32BIT(ibuf[i + tracks * j], effp->clips);
 285
 286         if (!whole_window)
 287             continue;
 288         else
 289             process_window(effp, data, (unsigned) i, (unsigned) tracks, obuf, (unsigned) (oldbuf + ncopy));
 290     }
 291
 292     *isamp = tracks*ncopy;
 293     if (whole_window)
 294         *osamp = tracks*(WINDOWSIZE/2);
 295     else
 296         *osamp = 0;
 297
 298     return SOX_SUCCESS;
 299 }
 300
 301 /*
 302  * We have up to half a window left to dump.
 303  */
 304
 305 static int sox_noisered_drain(sox_effect_t * effp, sox_sample_t *obuf, size_t *osamp)
 306 {
 307     priv_t * data = (priv_t *)effp->priv;
 308     unsigned i;
 309     unsigned tracks = effp->in_signal.channels;
 310     for (i = 0; i < tracks; i ++)
 311         *osamp = process_window(effp, data, i, tracks, obuf, (unsigned) data->bufdata);
 312
 313     /* FIXME: This is very picky.  osamp needs to be big enough to get all
 314      * remaining data or it will be discarded.
 315      */
 316     return (SOX_EOF);
 317 }
 318
 319 /*
 320  * Clean up.
 321  */
 322 static int sox_noisered_stop(sox_effect_t * effp)
 323 {
 324     priv_t * data = (priv_t *) effp->priv;
 325     size_t i;
 326
 327     for (i = 0; i < effp->in_signal.channels; i ++) {
 328         chandata_t* chan = &(data->chandata[i]);
 329         free(chan->lastwindow);
 330         free(chan->window);
 331         free(chan->smoothing);
 332         free(chan->noisegate);
 333     }
 334
 335     free(data->chandata);
 336
 337     return (SOX_SUCCESS);
 338 }
 339
 340 static sox_effect_handler_t sox_noisered_effect = {
 341   "noisered",
 342   "[profile-file [amount]]",
 343   SOX_EFF_MCHAN|SOX_EFF_LENGTH,
 344   sox_noisered_getopts,
 345   sox_noisered_start,
 346   sox_noisered_flow,
 347   sox_noisered_drain,
 348   sox_noisered_stop,
 349   NULL, sizeof(priv_t)
 350 };
 351
 352 const sox_effect_handler_t *lsx_noisered_effect_fn(void)
 353 {
 354     return &sox_noisered_effect;
 355 }