REAPERDenoiser.eel

desc:Ratio Denoiser (github.com/nbickford/REAPERDenoiser)
// That's the description of the plugin. This is how it'll show up in the effect
// search dialog, as well as the text at the start of its user interface. We use
// it as the first line of the script per the JSFX documentation's
// recommendation (https://www.reaper.fm/sdk/js/js.php#js_file)

// Define our user interface.
// Our FFT size will always be the same, so we only need controls for
// the noise collection mode and the noise scale (k).


///IFTEST* //this section not used by testing code
// 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:0<0.0,10.0,0.001>Noise Scale

// This defines a combo box for selecting the FFT size used in the STFT process. 
// The options range from 256 to 32768, allowing the user to choose the resolution 
// of the Fourier transform. The default value is set to the third option, 1024.
// This setting affects the granularity and latency of the audio processing, 
// with larger sizes providing finer frequency resolution but increased processing time.
slider3:5<0,7,1{256,512,1024,2048,4096,8192,16384,32768}>FFT size

// This defines a combo box for selecting the overlap factor used in the STFT process. 
// The options are 2, 4, 8, and 16, representing the division of the FFT size 
// to determine the overlap between successive FFT windows. The default value is 1, 
// corresponding to an overlap factor of 4. This setting affects the smoothness 
// and quality of the audio processing, with higher overlaps providing smoother results 
// at the cost of increased computational load.
slider4:1<0,3,1{2,4,8,16}>FFT overlap

//IFTEST*/ //this section not used by testing code

/*IFJSFX{
// 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
}IFJSFX*/



///IFTEST*
@init
//IFTEST*/
///IFJSFX*
slider1=0;
slider2=0;
slider3=5;
slider4=1;
//IFJSFX*/

//IFTEST srate = 24000;

// Memory allocation helper function
_free_memory = 0; // Global variable to track allocated memory
function simple_alloc(amount) local(_free_memory_old) global(_free_memory) (
  _free_memory_old = _free_memory;
  _free_memory += amount;
  _free_memory_old; // Return the starting address of the allocated block
);

function reset_state() (
	// clear buffer state
	memset(outputBuffer, 0, MAX_FFT_SIZE*2);
	memset(inputBuffer, 0, MAX_FFT_SIZE*2);
	samplesCollected = 0;
	fftCounter = 0;
	silence = overlaps;
);

function window(r) (
	sin(r*$pi)*sqrt(2);
);

function gen_window() (
	i = 0.5;
	loop(fftSize,
		windowBuffer[i] = window(i/fftSize);
		i += 1;
	);
);

function slider_code() (
	// 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(noiseBuffer, 0, MAX_FFT_SIZE+2);
			previousMode = 1;
		)
	) : previousMode = 0;
	changed = 0;
	fftSize = 2^(slider3+8);
	changed |= fftSize != oldFftSize;
	oldFftSize = fftSize;
	overlaps = 2^(slider4+1);
	changed |= overlaps != oldOverlaps;
	oldOverlaps = overlaps;
	fftInterval = fftSize/overlaps;
	fftScalingFactor = 1/overlaps/fftSize;
	kSquared = sqr(slider2); // slider2 is the Noise Scale from above.
	changed ? (
		reset_state();
		gen_window();
	);
	pdc_delay = fftSize; 
	pdc_bot_ch=0;
	pdc_top_ch=2;
);

// On initialization, initialize all of our variables.
MAX_FFT_SIZE = 32768; // Maximum FFT size supported

// Allocate memory for the complex FFT buffer, inputs, outputs, and noise spectrum
// Note: Each "complex" sample requires two slots: one for the real part and one for the imaginary part
fftBuffer = simple_alloc(MAX_FFT_SIZE*2); // Complex FFT buffer for processing
inputBuffer = simple_alloc(MAX_FFT_SIZE*2); // Complex input buffer
outputBuffer = simple_alloc(MAX_FFT_SIZE*2); // Complex output buffer
windowBuffer = simple_alloc(MAX_FFT_SIZE);
noiseBuffer = simple_alloc(MAX_FFT_SIZE+2); // Real buffer for noise spectrum

slider_code();

/*IFJSFX{
freembuf(_free_memory);

@slider
slider_code();
}IFJSFX*/


///IFTEST*
@sample
//IFTEST*/

function get_weight(yNorm, nNorm) (
	yNorm > 0 ? yNorm / (yNorm + kSquared*nNorm) : 0;
);

function process_stft_segment(fftBuffer, fftSize) (
	// If slider1 is greater than 0.5 (i.e. the user selected "Record Noise
	// Sample", we store the FFTs of each of these buffers.
	slider1 > 0.5? (
		// 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
		fft_bin = 0;
		loop(fftSize/2+1,
			fft_bin2 = fft_bin ? (fftSize - fft_bin) : 0;
			left_real = fftBuffer[2*fft_bin] + fftBuffer[2*fft_bin2];
			left_imag = fftBuffer[2*fft_bin + 1] - fftBuffer[2*fft_bin2 + 1];
			right_real = fftBuffer[2*fft_bin + 1] + fftBuffer[2*fft_bin2 + 1];
			right_imag = -fftBuffer[2*fft_bin] + fftBuffer[2*fft_bin2];

			normSquareNew_left = sqr(left_real) + sqr(left_imag);
			normSquareNew_right = sqr(right_real) + sqr(right_imag);

			normSquareOld_left = noiseBuffer[2*fft_bin];
			normSquareNew_left >= normSquareOld_left ? (
				noiseBuffer[2*fft_bin] = normSquareNew_left;
			);

			normSquareOld_right = noiseBuffer[2*fft_bin+1];
			normSquareNew_right >= normSquareOld_right ? (
				noiseBuffer[2*fft_bin+1] = normSquareNew_right;
			);
			fft_bin += 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 = 0 to SIZE - 1.
	fft_bin = 0;
	loop(fftSize/2+1,
		fft_bin2 = fft_bin ? (fftSize - fft_bin) : 0;
		// unfold complex spectra
		left_real = fftBuffer[2*fft_bin] + fftBuffer[2*fft_bin2];
		left_imag = fftBuffer[2*fft_bin + 1] - fftBuffer[2*fft_bin2 + 1];
		right_real = fftBuffer[2*fft_bin + 1] + fftBuffer[2*fft_bin2 + 1];
		right_imag = -fftBuffer[2*fft_bin] + fftBuffer[2*fft_bin2];

		// apply denoising algorithm to find a multiplier for these spectra.
		weight_left = get_weight(sqr(left_real) + sqr(left_imag), noiseBuffer[2*fft_bin]);
		weight_right = get_weight(sqr(right_real) + sqr(right_imag), noiseBuffer[2*fft_bin+1]);

		// multiply spectra
		left_real *= weight_left;
		left_imag *= weight_left;
		right_real *= weight_right;
		right_imag *= weight_right;

		// re fold complex spectra
		fftBuffer[2*fft_bin] = (left_real - right_imag)*0.5;
		fftBuffer[2*fft_bin + 1] = (left_imag + right_real)*0.5;
		fftBuffer[2*fft_bin2] = (left_real + right_imag)*0.5;
		fftBuffer[2*fft_bin2 + 1] = (-left_imag + right_real)*0.5;
		fft_bin += 1;
	);
);

//IFTEST function sample_code() (

// input 1 stereo sample
inputBuffer[samplesCollected*2] = spl0;
inputBuffer[samplesCollected*2+1] = spl1;

// output 1 stereo sample
silence > 0 ? (
	spl0 = spl1 = 0;
	// silence for fft init
) : (
	spl0 = outputBuffer[samplesCollected*2];
	spl1 = outputBuffer[samplesCollected*2+1];
);

// clear 1 output buffer sample
outputBuffer[samplesCollected*2] = outputBuffer[samplesCollected*2+1] = 0;


fftCounter += 1;

fftCounter >= fftInterval ? (
	fftCounter = 0;

	// copy input to fft
	bufferIndex = samplesCollected + 1;
	i = 0;
	loop(fftSize,
		bufferIndex >= fftSize ? bufferIndex -= fftSize;
		fftBuffer[i*2] = inputBuffer[bufferIndex*2]*windowBuffer[i];
		fftBuffer[i*2+1] = inputBuffer[bufferIndex*2+1]*windowBuffer[i];
		bufferIndex += 1;
		i += 1;
	);

	// process
	fft(fftBuffer, fftSize);
	fft_permute(fftBuffer, fftSize);

	process_stft_segment(fftBuffer, fftSize);

	fft_ipermute(fftBuffer, fftSize);
	ifft(fftBuffer, fftSize);

	// copy fft back

	i = 0;
	bufferIndex = samplesCollected + 1;
	loop(fftSize,
		bufferIndex >= fftSize ? bufferIndex -= fftSize;
		outputBuffer[2*bufferIndex] += fftBuffer[2*i]*fftScalingFactor*windowBuffer[i];
		outputBuffer[2*bufferIndex+1] += fftBuffer[2*i+1]*fftScalingFactor*windowBuffer[i];
		bufferIndex += 1;
		i += 1;
	);
	silence > 0 ? silence -= 1;
);

samplesCollected += 1;
samplesCollected %= fftSize;
//IFTEST ); // function sample_code()


/*IFJSFX{
@serialize
}IFJSFX*/
//IFTEST function file_var(file, val) (printf("FILE_VAR, FILE: %d, VAL: %g\n", file, val));
//IFTEST function file_mem(file, val, count) (printf("FILE_MEM, FILE: %d, MEM: %g[0..%d]\n", file, val, count-1));
/*
// 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, noiseBuffer, MAX_FFT_SIZE+2);
*/
/*IFTEST{ // main test block

// helpers
function sum_first_pdc_samples(s0val, s1val) (
	reset_state();
	spl0sum = 0;
	spl1sum = 0;
	loop(pdc_delay,
		spl0=s0val;
		spl1=s1val;
		sample_code();
		spl0sum += abs(spl0);
		spl1sum += abs(spl1);
	);
	assert_near_equal(0, 0.00000001, spl0sum, "spl0sum for init");
	assert_near_equal(0, 0.00000001, spl1sum, "spl1sum for init");
	spl0sum = 0;
	spl1sum = 0;
);

printf("SAMPLE_0_0_TEST\n");
// spl0=0 spl1=0 for 100 output samples
sum_first_pdc_samples(0, 0);
loop(100,
	spl0=0;
	spl1=0;
	sample_code();
	spl0sum += abs(spl0);
	spl1sum += abs(spl1);
);
assert_near_equal(0, 0.001, spl0sum, "SAMPLE_0_0_TEST: spl0sum was not as expected");
assert_near_equal(0, 0.001, spl1sum, "SAMPLE_0_0_TEST: spl1sum was not as expected");
printf("spl0sum=%g, spl1sum=%g\n", spl0sum, spl1sum);

printf("SAMPLE_1_0_TEST\n");
// spl0=1 spl1=0 for for 100 output samples
sum_first_pdc_samples(1, 0);
loop(100,
	spl0=1;
	spl1=0;
	sample_code();
	spl0sum += abs(spl0);
	spl1sum += abs(spl1);
);
assert_near_equal(100, 0.001, spl0sum, "SAMPLE_1_0_TEST: spl0sum was not as expected");
assert_near_equal(0, 0.001, spl1sum, "SAMPLE_1_0_TEST: spl1sum was not as expected");
printf("spl0sum=%g, spl1sum=%g\n", spl0sum, spl1sum);

printf("SAMPLE_0_1_TEST\n");
// spl0=0 spl1=1 for for 100 output samples
sum_first_pdc_samples(0, 1);
loop(100,
	spl0=0;
	spl1=1;
	sample_code();
	spl0sum += abs(spl0);
	spl1sum += abs(spl1);
);
assert_near_equal(0, 0.001, spl0sum, "SAMPLE_0_1_TEST: spl0sum was not as expected");
assert_near_equal(100, 0.001, spl1sum, "SAMPLE_0_1_TEST: spl1sum was not as expected");
printf("spl0sum=%g, spl1sum=%g\n", spl0sum, spl1sum);

printf("SAMPLE_1_1_TEST\n");
// spl0=1 spl1=1 for for 100 output samples
sum_first_pdc_samples(1, 1);
loop(100,
	spl0=1;
	spl1=1;
	sample_code();
	spl0sum += abs(spl0);
	spl1sum += abs(spl1);
);
assert_near_equal(100, 0.001, spl0sum, "SAMPLE_1_1_TEST: spl0sum was not as expected");
assert_near_equal(100, 0.001, spl1sum, "SAMPLE_1_1_TEST: spl1sum was not as expected");
printf("spl0sum=%g, spl1sum=%g\n", spl0sum, spl1sum);

test_summary();
}IFTEST*/ // main test block