fftdconvEfficient.jsfx

desc:Convolution Ratio Denoiser

// loosely based off github.com/nbickford/REAPERDenoiser

// This defines a combo box that allows the user to select "Denoise Input" or
// "Record Noise Sample". The default value is 0 (Denoise Input). The maximum
// value is 1 (Record Noise Sample), and it increases in steps of 1.
slider1:0<0,1,1{Denoise Input, Record Noise Sample}>Noise Collection Mode

// This defines a slider that can be varied between 0.0 and 10.0 in steps of
// 0.001, with default value 1.0. (If slider2 is equal to 0.0, this plugin
// shouldn't really do anything to the input audio.)
slider2:1.0<1.0,100.0,0.001>Noise Scale

// This defines a slider that specifies the max level (in dB)
// of the noise that is recorded
slider3:-120<-120,0,0.1>Noise detect min (dB)

// This defines a slider that specifies the minimum level (in dB)
// that a noise is recorded with (useful to filter out silence)
slider4:0<-120,0,0.1>Noise detext max (dB)

// This defines a slider to mix between no noise reduction or
// any amount of noise reduction, up to 100%
slider5:100<0,100,0.5>Wet/dry mix (0 = all dry, 100 = all wet)

slider6:0<0,10,0.01>Decay factor

slider7:1<0.01,1,0.01>Tail Size

// Here we can label our input and output pins. This also tells REAPER how many
// channels we can handle. In this case, the plugin is stereo (a monophonic
// plugin would be simpler, but I almost always use this to denoise stereo
// audio), so we define two input and output pins.
in_pin:Noisy Audio 1
in_pin:Noisy Audio 2
out_pin:Denoised Audio 1
out_pin:Denoised Audio 2


import cookdsp.jsfx-inc

@gfx 1044 400
gfx_set(0, 1, 0);
i = 0;
gi = 10;
idiff = max(1, DSIZE/1024);
loop(min(1024, DSIZE),
  gfx_line(gi, 200, gi, 200-50*visBuf[i]);
  i += idiff;
  gi += 1;
);


@init
// don't add noise to an already silent signal
ext_nodenorm = 1;
ext_noinit = 1;

// function to convert decibel measurement to decimal value
function dBtoDec(value)
(
  2^(value/6);
);
// On initialization, initialize all of our variables.

function memalloc_(sz) (
    // compensate for cookdsp adding some extra space
    // (avoids problems with memory layout when using
    //    FFT/MDCT and trying to manually keep track of the layout)
    __memory_next -= 8;
    memalloc(sz);
);


function memalloc_(sz, reinit) (
    // compensate for cookdsp adding some extra space
    // (avoids problems with memory layout when using
    //    FFT/MDCT and trying to manually keep track of the layout)
    __memory_next -= 8;
    memalloc(sz, reinit);
);

// The FFT size will always be constant.
SIZE = 8192;
DSIZE = 2*SIZE;
invSIZE = 1/SIZE;
invDSIZE = 0.5/SIZE;
HALFSIZE = SIZE*0.5;

convBufferL = memalloc_(DSIZE);
convBufferR = memalloc_(DSIZE);

fftOutL = memalloc_(DSIZE);
fftOutR = memalloc_(DSIZE);

visBuf = memalloc_(DSIZE);

window_buf = memalloc_(SIZE);

outBufL = memalloc_(SIZE+1);
outBufR = memalloc_(SIZE+1);

delL = memalloc_(HALFSIZE);
delR = memalloc_(HALFSIZE);

noiseBufferL = memalloc_(HALFSIZE+1, 1);
noiseBufferR = memalloc_(HALFSIZE+1, 1);

// fill window_buf with a hahn window
i = 0;
loop(SIZE,
  window_buf[i] = 0.5 - 0.5*cos(2*$pi*i/SIZE);
  i += 1;
);

runningSumL = memalloc_(1);
runningSumR = memalloc_(1);

freembuf(__memory_next);

@slider
// A simple function to zero out the noise buffers when switching mode to "Record Noise Sample"
// previousMode should default to 0 on first initialization, but setting it to 0 in @init will cause
// this code to get run again, and the noise profile lost even when switching to "Denoise Input"
slider1 > 0.5 ? (
  previousMode > 0.5 ? (
    memset(noiseBufferL, 0, SIZE+2);
    previousMode = 0;
  )
) : previousMode = 1;
noiseLow = dBtoDec(slider3);
noiseHigh = dBtoDec(slider4);
//dryfrac = slider5*0.01;
//wetfrac = 1 - slider5*0.02;
dryfrac = 1 - slider5/100;
wetfrac = slider5/100;
strength = sqr(slider2);
df = exp(slider6);
tail = 0|max(4, min(SIZE, SIZE*slider7));

ss = 1 + sqrt(strength);
invss = 1 / ss;

pdc_delay = HALFSIZE;
pdc_bot_ch = 0;
pdc_top_ch = 2;


@sample

// We'll write a function to denoise a single channel, and then we'll call this
// for each of the channels.
// In this case, we'll pass in the channel number, the four input and output
// tiles, and the current sample.
// We also need to specify which variables will be local to the function (i.e.
// which variables have local instead of global scope).
// Note that channels are zero-indexed (so the left channel is channel 0, and
// the right channel is channel 1).
// Functions can return values, but this one won't return anything.
// Swapping tiles and resetting samplesCollected will be managed by the caller.
function denoiseChannel(channel fftBuffer inBuf outBuf noiseBuffer runningSum convBuffer samplesCollected)
(
  // Read out input audio
  sample = spl(channel); // You can also use spl0 or spl1.
  runningSum[] += sample * sample;
  
  dry = inBuf[samplesCollected];

  wet = outBuf[samplesCollected];

  inBuf[samplesCollected] = sample;

  outBuf[samplesCollected] = outBuf[samplesCollected + HALFSIZE];
  
  fftBuffer[samplesCollected*2] = sample*invSIZE;
  fftBuffer[samplesCollected*2+1] = 0;
  fftBuffer[samplesCollected*2+SIZE] = 0;
  fftBuffer[samplesCollected*2+SIZE+1] = 0;
  
  /*
  bandIndex = 0;
  loop(HALFSIZE,
    outBuf[bandIndex] = outBuf[bandIndex+1] + convBuffer[2*bandIndex] * sample;
    bandIndex += 1;
  );*/

  outVal = wet * wetfrac + dry * dryfrac;

  // Reached end of a tile?
  samplesCollected >= HALFSIZE-1 ? (
    fft(fftBuffer, SIZE);
    
    // If slider1 is greater than 0.5 (i.e. the user selected "Record Noise
    // Sample", we store the FFT mangitudes of each of these buffers, but only when the level is correct
    runningSumVal = sqrt(runningSum[]*invSIZE);
    runningSum[] = 0;
    (slider1 > 0.5 && runningSumVal >= noiseLow
      && runningSumVal <= noiseHigh) ?
    (
      // for each band, compare the norm of the noise in this frame.
      // If it is greater than what's already there for this band, then copy
      // it into the noiseBuffer
      // treat dc and ny bands seperately:
      // band 0 (dc)
      index = 0;
      loop(HALFSIZE+1,
        squareMagnitudeNew = sqr(fftBuffer[2*index]) + sqr(fftBuffer[2*index+1]);
        squareMagnitudeOld = noiseBuffer[index];
        noiseBuffer[index] = 0.25 * squareMagnitudeNew + 0.75 * squareMagnitudeOld;
        index += 1;
      );
    );

    // Apply Norbert Weiner's filtering algorithm,
    //   X(f) = Y(f) * (|Y(f)|^2)/(|Y(f)|^2 + k^2 |N(f)|^2)
    // sqr() computes the square of a number, and abs() computes the absolute
    // value of a number. We also include a factor of 1/SIZE, to normalize the
    // FFT (so that if we don't do any denoising, the input signal is equal to
    // the output signal).

    // Loop over each band, from bandIndex = 1 to HALFSIZE - 1.
    // caluclate attenuationfactor and store in convBuffer as a real valued complex number.
    bandIndex = 0;
    loop(SIZE,
      // Compute |Y(f)|^2 = real(Y(f))^2 + imaginary(Y(f))^2
      yNorm = sqr(fftBuffer[bandIndex*2]) + sqr(fftBuffer[2*bandIndex + 1]);
      // The same for the noise component:
      nNorm = noiseBuffer[min(bandIndex, SIZE-bandIndex)];
      
      attenuationFactor = yNorm / (strength * nNorm + yNorm);
      
      ((yNorm + nNorm)*100000000 == 0) ? attenuationFactor = 1;

      convBuffer[2*bandIndex] = attenuationFactor;
      convBuffer[2*bandIndex+1] = 0;
      
      bandIndex += 1;
    );
    
    /*
    bandIndex = 0;
    loop(DSIZE,
      visBuf[bandIndex] = convBuffer[bandIndex];
      bandIndex += 1;
    );
    //*/
    
    // find min of convBuffer
    minVal = 2;
    bandIndex = 0;
    loop(SIZE,
      convBuffer[bandIndex] > 0 ? minVal = min(convBuffer[bandIndex], minVal);
      bandIndex += 2;
    );
    
    // go a bit past the min
    minVal *= 0.0000001;
    
    // add a bias to all bands, to prevent taking log of 0,
    // and take log
    bandIndex = 0;
    loop(SIZE,
      // calculate log(|H(k)|^2)/4 == log(|H(k)|)/2
      convBuffer[bandIndex] = log(convBuffer[bandIndex]+minVal)*invSIZE;
      bandIndex += 2;
    );
    
    /*
    bandIndex = 0;
    loop(DSIZE,
      visBuf[bandIndex] = convBuffer[bandIndex]*SIZE;
      bandIndex += 1;
    );
    //*/

    // convert to cepstrum
    ifft(convBuffer, SIZE);

    /*
    bandIndex = 0;
    loop(DSIZE,
      visBuf[bandIndex] = convBuffer[bandIndex];
      bandIndex += 1;
    );
    //*/
    
    convBuffer[0] *= 0.5;
    /*bandIndex = 1;
    loop(SIZE-1,
      convBuffer[bandIndex] *= 1;// * window_buf[bandIndex];
      bandIndex += 1;
    );*/

    memset(convBuffer+SIZE, 0, SIZE);
    
    /*
    bandIndex = 0;
    loop(SIZE,
      visBuf[bandIndex] = convBuffer[bandIndex]*SIZE;
      bandIndex += 1;
    );
    //*/
    
    // reverse cepstrum into fft realm
    fft(convBuffer, SIZE);

    /*
    bandIndex = 0;
    loop(DSIZE,
      visBuf[bandIndex] = convBuffer[bandIndex];
      bandIndex += 1;
    );
    //*/

    // apply complex exp to reverse log from before
    bandIndex = 0;
    loop(SIZE,
      // e^(a+bj) = e*a(cos(b) + j*sin(b))
      mult = exp(convBuffer[bandIndex])*invSIZE;
      ph = convBuffer[bandIndex+1];
      convBuffer[bandIndex+1] = mult * sin(ph);
      convBuffer[bandIndex] = mult * cos(ph);
      bandIndex += 2;
    );

    // clean convBuffer
    ifft(convBuffer, SIZE);
    
    
    //bandIndex = 0;
    //loop(HALFSIZE,
    //  convBuffer[bandIndex] *= window_buf[bandIndex+HALFSIZE];
    //);
    
    /*
    bandIndex = 0;
    loop(DSIZE,
      visBuf[bandIndex] = convBuffer[bandIndex]*SIZE;
      bandIndex += 1;
    );
    //*/
    
    bandIndex = 0;
    loop(tail,
      w = (0.5+0.5*cos($pi*bandIndex/tail))^df;
      convBuffer[bandIndex] *= w;
      convBuffer[bandIndex+1] *= w;
      bandIndex += 2;
    );
    
    // discard second half of the fir
    memset(convBuffer + tail, 0, DSIZE - tail);
  
    //bandIndex = 0;
    //loop(SIZE,
    //  convBuffer[bandIndex] *= window_buf[1];
    //  bandIndex += 2;
    //);
    //convBuffer[0] = 1;
    
    
    convBuffer[0] *= invss;
    //*
    bandIndex = 0;
    loop(DSIZE,
      visBuf[bandIndex] = convBuffer[bandIndex*0.5]*44*ss;
      bandIndex += 1;
    );
     //*/

    fft(convBuffer, SIZE);
    
    /*
    bandIndex = 0;
    loop(SIZE,
      visBuf[bandIndex+SIZE] = convBuffer[bandIndex*2]*SIZE;
      bandIndex += 1;
    );
    //*/

    convolve_c(fftBuffer, convBuffer, SIZE);
    
    ifft(fftBuffer, SIZE);
    
    // copy to outBuf
    bandIndex = 0;
    loop(HALFSIZE,
      outBuf[bandIndex] += fftBuffer[bandIndex*2];
      bandIndex += 1;
    );
    loop(HALFSIZE,
      outBuf[bandIndex] = fftBuffer[bandIndex*2];
      bandIndex += 1;
    );
  ); // samplesCollected >= HALFSIZE-1
  outVal;
); // function denoiseChannel(...

// Now, call denoiseChannel for each of the channels.
spl0 = denoiseChannel(0, fftOutL, delL, outBufL, noiseBufferL, runningSumL, convBufferL, samplesCollected);
spl1 = denoiseChannel(1, fftOutR, delR, outBufR, noiseBufferR, runningSumR, convBufferR, samplesCollected);
// Go to the next sample
samplesCollected += 1;
samplesCollected %= HALFSIZE;


@serialize
// Sliders are serialized automatically, so all we have to serialize is the two
// noise buffers. JSFX's serialization works in a clever way: when reading the
// state of the plugin from a serialized version, these functions copy data into
// noiseBufferL and noiseBufferR. But when writing out the state of the plugin,
// they work the other way, copying data out of noiseBufferL and noiseBufferR.
file_var(0, previousMode);
file_mem(0, noiseBufferL, SIZE+2);