denconv-alt.jsfx

desc:Alternate algorithm convolution ratio denoising

// 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.0,10.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.
slider5:100<0,100,0.5>Wet/dry mix (0 = all dry, 100 = all wet)

// checkbox? to enable/disable difference monitoring
slider6:0<0,1,1{Off, On}>Difference monitoring

// 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 512 200
gfx_set(1, 0, 0);
bi = 0;
loop(min(512, SIZE/2),
  gfx_line(bi, 0, bi, 200*tileBufL[bi]);
  bi += 1;
);

@init
// 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 ceps_permute(buf, size)
local(i, temp)
(
  // normally swap i and i + SIZE/2, but instead
  // swap each adjacent pair, starting from 0, 2, ...
  // because the permutation is bit-reversed after fft.
  i = 0;
  loop(size/2,
    // swap Re for these two bands
    temp = buf[i];
    buf[i] = buf[size + 2];
    buf[size + i] = temp;
    i += 1;
    // swap Im for these two bands
    temp = buf[i];
    buf[i] = buf[size + 2];
    buf[size + i] = temp;
    i += 2;
  );
);

// The FFT size will always be constant.
SIZE = 2048;
DSIZE = 2*SIZE;

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

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

tempfft = memalloc_(2*DSIZE); // 

//fftOutL1 = fftOutL;
//fftOutL2 = fftOutL + DSIZE;

//fftOutR1 = fftOutR;
//fftOutR2 = fftOutR + DSIZE;

outLbuf = memalloc_(SIZE);
outRbuf = memalloc_(SIZE);

noiseBufferL = memalloc_(SIZE);
noiseBufferR = memalloc_(SIZE);

tileBufL = memalloc_(SIZE);
tileBufR = memalloc_(SIZE);

delBufL = memalloc_(SIZE/2);
delBufR = memalloc_(SIZE/2);

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


// freembuf is only needed if the memory ever gets assigned to, and then can be freed.

// samplesCollected will be our position in *_history
samplesCollected = 0;

// Finally, the algorithm we use outputs modified audio SIZE samples after we
// input it. If we tell REAPER that the plugin has a delay of SIZE samples,
// REAPER can automatically compensate for this and make it appear as if there's
// no delay at all.
pdc_delay=SIZE/2; 
pdc_bot_ch=0;
pdc_top_ch=2;


@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, DSIZE);
    previousMode = 1;
  )
) : previousMode = 0;
noiseLow = dBtoDec(slider3);
noiseHigh = dBtoDec(slider4);
slider6 == 0 ?
(
  wetfrac = slider5/100;
  dryfrac = (100-slider5)/100;
) : (
  dryfrac = 1;
  wetfrac = - slider5/100;
);


@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 tilebuf inDelBuf fftBuffer outBuf convBuffer noiseBuffer runningSum samplesCollected)
(
  // Read out input audio
  sample = spl(channel); // You can also use spl0 or spl1.
  runningSum[] += sample * sample;
  
  dry = inDelBuf[samplesCollected];
  
  inDelBuf[samplesCollected] = sample;
  
  fftBuffer[samplesCollected*2] = sample;
  fftBuffer[samplesCollected*2+1] = 0;

  // Reached end of a tile? - update convBuffer
  samplesCollected >= SIZE/2-1 ? (
    memset(fftBuffer + SIZE, 0, SIZE); 
    fft(fftBuffer, SIZE);
    fft_permute(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(2*runningSum[]/SIZE);
    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
      index = 0;
      loop(SIZE,
        squareMagnitudeNew = sqr(fftBuffer[2*index]) + sqr(fftBuffer[2*index+1]);
        squareMagnitudeOld = noiseBuffer[index];
        noiseBuffer[index] = (squareMagnitudeNew + 3 * squareMagnitudeOld) / 4;
        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).
    kSquared = sqr(slider2); // slider2 is the Noise Scale from above.

    // Loop over each band, from bandIndex = 0 to SIZE - 1.
    bandIndex = 0;
    loop(SIZE,
      // Compute |Y(f)|^2 = real(Y(f))^2 + imaginary(Y(f))^2
      yNorm = sqr(fftBuffer[bandIndex]) + sqr(fftBuffer[bandIndex + 1]);
      // The same for the noise component:
      nNorm = noiseBuffer[0.5*bandIndex];
      
      attenuationFactor = yNorm / (yNorm + kSquared * nNorm);
      
      ((yNorm + nNorm)*100000000 == 0) ? attenuationFactor = 1;

      tileBuf[bandIndex*0.5] = attenuationFactor;
      
      convBuffer[bandIndex] = attenuationFactor / SIZE;
      bandIndex += 2;
    );
    
    ZZZ1 = convBuffer[0];
    
    /*bandIndex = 0;
    loop(SIZE,
      // calculate log(|H(k)|^2)/4
      ov = log(sqr(convBuffer[bandIndex]) + sqr(convBuffer[bandIndex+1]));
      convBuffer[bandIndex] = max(-1099511627776, ov)*0.25/SIZE;
      bandIndex += 1;
    );

    // prepare to take cepstrum
    ceps_ipermute(convBuffer, SIZE);
    ZZZ1a = convbuffer[2];

    // take ifft, converting to cepstrum
    ifft(convBuffer, SIZE);
    
    //fft_ipermute(convBuffer, SIZE);
    
    // homomorphic filter
    //memset(convBuffer+SIZE, 0, SIZE-2);
    convBuffer[0] *= -1;
    bandIndex = 2;
    loop(SIZE-1,
      convBuffer[bandIndex] *= 2;
      bandIndex += 2;
    );
    
    //fft_permute(convBuffer, SIZE);
    
    // take FFT to reverse cepstrum into fft realm
    fft(convBuffer, SIZE);

    // re-permute
    ceps_permute(convBuffer, SIZE);

    ZZZ2 = convBuffer[SIZE] / SIZE;

    // 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] / SIZE) / SIZE;
      convBuffer[bandIndex]   = mult * cos(convBuffer[bandIndex+1]);
      convBuffer[bandIndex+1] = mult * sin(convBuffer[bandIndex+1]);
      bandIndex += 2;
    );

    ceps_ipermute(convBuffer, SIZE);
    
    
    
    // take convBuffer back to time domain to extract FIR
    /*/
    fft_ipermute(convBuffer, SIZE);
    ifft(convBuffer, SIZE);
    
    
    ZZZ4 = convBuffer[0];
    //ZZZ3 = convBuffer[0];
    
    // discard second half of the fir
    // memset(convBuffer + SIZE, 0, SIZE);
    
    //fft(convBuffer, SIZE);
    
    //ZZZ4 = convBuffer[0]; // expect this to be 1 when no denoising 

    // perform covolution
    //convolve_c(fftBuffer, convBuffer, SIZE);

    // convert fft to real convolved signal
    //ifft(fftBuffer, SIZE);
    // add first half of fftBuffer
    /*bandIndex = 0;
    loop(SIZE/2,
      outBuf[bandIndex] += fftBuffer[bandIndex*2];
      bandIndex += 1;
    );
    // replace second half
    loop(SIZE/2,
      outBuf[bandIndex] = fftBuffer[bandIndex*2];
      bandIndex += 1;
    );*/
  ); // samplescollected >= SIZE/2 - 1
  
  
  bandIndex = 0;
  loop(SIZE/4,
    outBuf[bandIndex] = outBuf[bandIndex + 1] + dry * convBuffer[bandIndex*2];
    bandIndex += 1;
  );
  wet = outBuf[0];
  
  spl(channel) = wet * wetfrac + dry * dryfrac;
  
  //// copy from second half
  //outBuf[samplesCollected] = outBuf[samplesCollected+SIZE/2];
);

// Now, call denoiseChannel for each of the channels.
denoiseChannel(0, tileBufL, delBufL, fftOutL, outLbuf, convBufferL, noiseBufferL, runningSumL, samplesCollected);
//denoiseChannel(1, tileBufR, delBufR, fftOutR, outRbuf, convBufferR, noiseBufferR, runningSumR, samplesCollected);
// Go to the next sample
samplesCollected += 1;
samplesCollected %= SIZE/2;

@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_mem(0, noiseBufferL, DSIZE);