From 6e0c4f9d7142212a90d40a934feecc284524e922 Mon Sep 17 00:00:00 2001 From: micsthepick <11528421+micsthepick@users.noreply.github.com> Date: Thu, 29 Feb 2024 13:21:49 +1100 Subject: [PATCH] fix Blurry, better testing should fix some glitches --- .gitignore | 3 +- process_scripts.sh | 16 + test_eel.sh | 76 +- vocalrediso-jamesdsp.eel => vocalrediso.eel | 897 +++++++++++--------- vocalrediso.jsfx | 233 ++++- vocalrediso.jsfx-template | 478 +++++++++++ vocalredisoBlurry.eel | 517 +++++++++++ vocalredisoBlurry.jsfx | 324 +++++-- vocalredisoBlurry.jsfx-template | 517 +++++++++++ 9 files changed, 2502 insertions(+), 559 deletions(-) create mode 100755 process_scripts.sh rename vocalrediso-jamesdsp.eel => vocalrediso.eel (69%) create mode 100644 vocalrediso.jsfx-template create mode 100644 vocalredisoBlurry.eel create mode 100644 vocalredisoBlurry.jsfx-template diff --git a/.gitignore b/.gitignore index 43d5d8a..8b71034 100644 --- a/.gitignore +++ b/.gitignore @@ -361,4 +361,5 @@ MigrationBackup/ # Fody - auto-generated XML schema FodyWeavers.xsd -*.test.eel2 \ No newline at end of file +*.test.eel2 +*.test \ No newline at end of file diff --git a/process_scripts.sh b/process_scripts.sh new file mode 100755 index 0000000..cab7ff4 --- /dev/null +++ b/process_scripts.sh @@ -0,0 +1,16 @@ +function process_eel_script { + local input_script="$1" + local output_script="$2" + local replace_str="$3" # New argument for the replacement + + # Process the second script and append its contents to the output script + sed -r -e "s/\/\/IF${replace_str}|\/\*IF${replace_str}\{|\}IF${replace_str}\*\///" \ + "$input_script" > "$output_script" +} + +process_eel_script vocalrediso.jsfx-template vocalrediso.test TEST +process_eel_script vocalrediso.jsfx-template vocalrediso.eel EEL +process_eel_script vocalrediso.jsfx-template vocalrediso.jsfx JSFX +process_eel_script vocalredisoBlurry.jsfx-template vocalredisoBlurry.test TEST +process_eel_script vocalredisoBlurry.jsfx-template vocalredisoBlurry.eel EEL +process_eel_script vocalredisoBlurry.jsfx-template vocalredisoBlurry.jsfx JSFX \ No newline at end of file diff --git a/test_eel.sh b/test_eel.sh index dde23af..dbb5e31 100755 --- a/test_eel.sh +++ b/test_eel.sh @@ -1,44 +1,32 @@ -cat testing_defines.eel2 > vocalrediso.test.eel2 -sed -r -e 's/^(desc|slider[0-9]+):.*|^@(init|slider|block|serialize|sample)//' \ - -e 's/\/\/DEBUGPRINT/printf/' \ - -e 's/\/\/IFTEST|\/\*IFTEST\{\*|\*\}IFTEST\*\///' \ - -e 's/\*IFNTEST\*|IFNTEST\{\*|\*\}IFNTEST//' \ - vocalrediso-jamesdsp.eel >> vocalrediso.test.eel2 - -# Initialize a flag to indicate stderr output -stderr_output=0 - -output=$(./WDL/WDL/eel2/loose_eel ./vocalrediso.test.eel2 2>&1) - -echo "$output" - -# Get the last line of the output -last_line=$(echo "$output" | tail -n 1) - -# Check if the last line starts with 'FAILURE' -if echo "$last_line" | grep -q "^FAILURE"; then - echo "Failed Test Cases, will return -1!" - exit -1 -fi - -##//DEBUGPRINT("HI"); -##//IFTEST code_here(); -##/*IFTEST{* -## more_code(); -##*}IFTEST*/ -##/*IFNTEST*/called_when_not_testing(); -##/*IFNTEST{*/ -## also_called_when_not_testing(); -##/*}IFNTEST*/ - -# will transform to - -##printf("HI"); -## code_here(); -## -## more_code(); -## -##//called_when_not_testing(); -##/*IFNTEST{*/ -## also_called_when_not_testing(); -##/*IFNTEST}*/ \ No newline at end of file +# ensure latest testing script is available +./process_scripts.sh + +# add defines to head of test script +cat testing_defines.eel2 vocalrediso.test > vocalrediso.test.eel2 + +function run_test { + local test_script="$1" + local test_output + + echo "TESTING $1" + + # Ensure latest testing script is available + ./process_scripts.sh + + # Add defines to head of test script + cat testing_defines.eel2 "$test_script" > "${test_script}.eel2" + + # Run the test and capture output + test_output=$(./WDL/WDL/eel2/loose_eel "./${test_script}.eel2" 2>&1) + echo "$test_output" + + # Get the last line of the output and check for failure + if echo "$test_output" | tail -n 1 | grep -q "^FAILURE"; then + echo "Failed Test Cases, will return -1!" + exit -1 + fi +} + +# Run tests for each script +run_test vocalrediso.test +run_test vocalredisoBlurry.test diff --git a/vocalrediso-jamesdsp.eel b/vocalrediso.eel similarity index 69% rename from vocalrediso-jamesdsp.eel rename to vocalrediso.eel index 7719321..ccd3ffd 100644 --- a/vocalrediso-jamesdsp.eel +++ b/vocalrediso.eel @@ -1,419 +1,478 @@ -// Based on Vocalrediso.ny, a nyquist filter for audacity, currently released under the GNU GPL V3 -// see https://www.gnu.org/licenses/gpl-3.0.en.html, https://www.audacityteam.org/ -// credits to Neil Bickford for his denoiser, https://github.com/Nbickford/REAPERDenoiser, -// which was used as a starting point for this script -// which now uses the STFT template code by geraintluff: https://forum.cockos.com/showthread.php?t=225955 - -desc: vocal removal/isolation -//tags: processing vocals stereo -//author: Michael Pannekoek - -//slider1:fft_size_index=2<0,7,1{256,512,1024,2048,4096,8192,16384,32768}>FFT size -slider1:0<-100,100,0.1>dry mix -slider2:0<-100,100,0.1>C mix (Vocals) -slider3:0<-5,5,0.001>strength at Low Cut -slider4:0<-5,5,0.001>strength at High Cut -slider5:80<0,24000,10>Low Cut (Vocals) -slider6:24000<0,24000,10>High Cut (Vocals) -slider7:0<-90,90,0.1>Phase (Degrees) -slider8:90<1,180,0.1>Phase width at Low Cut (Degrees) -slider9:90<1,180,0.1>Phase width at High Cut (Degrees) -slider10:1<0,1,0.05>Attenuate if different volume -slider11:1<0,1,1{No,Yes}>undo input corrections -slider12:0<-180,180,0.05>Phase2 (Degrees) - - -@init -slider1=0; -slider2=0; -slider3=0; -slider4=0; -slider5=80; -slider6=24000; -slider7=0; -slider8=90; -slider9=90; -slider10=1; -slider11=1; -slider12=0; - -//IFTEST srate = 24000; - -_free_memory = 0; -function simple_alloc(amount) -local(_free_memory_old) -global(_free_memory) -( - _free_memory_old = _free_memory; - _free_memory += amount; - _free_memory_old; -); - -prev_fft_size = 0; - -function setup_stft_state(fft_size, first_time) ( - //////////////////////// Setup block - // This is called when playback starts, or when the FFT size is changed - //memset(fft_buffer, 0, MAX_FFT_SIZE*2); - memset(output_buffer, 0, buffer_length*2); - memset(input_buffer, 0, buffer_length*2); - // reset indexes and counters - buffer_index = 0; - output_index = 0; - fft_counter = 0; - // force silence for first overlaps where fft not yet run - silence = overlap_factor; - //////////////////////// -); - -setup_stft_state(fft_size, prev_fft_size == 0); - -function process_stft_segment(fft_buffer, fft_size) local(fft_bin, left_real, left_imag, right_real, right_imag) ( - fft_bin = 0; // FFT bin number - loop(fft_size/2+1, - fft_bin2 = fft_bin ? (fft_size - fft_bin) : 0; - - // Unfold complex spectrum into two real spectra - left_real = fft_buffer[2*fft_bin] + fft_buffer[2*fft_bin2]; - left_imag = fft_buffer[2*fft_bin + 1] - fft_buffer[2*fft_bin2 + 1]; - right_real = fft_buffer[2*fft_bin + 1] + fft_buffer[2*fft_bin2 + 1]; - right_imag = -fft_buffer[2*fft_bin] + fft_buffer[2*fft_bin2]; - - // input corrections - // first, change the phase of L based on phase2: - l_real_twisted = left_real*cosine2 + left_imag*sine2; - l_imag_twisted = left_imag*cosine2 - left_real*sine2; - - // now mix L&R together based on phase - r_real_rotated = right_real*cosine + l_real_twisted*sine; - r_imag_rotated = right_imag*cosine + l_imag_twisted*sine; - l_real_rotated = l_real_twisted*cosine - right_real*sine; - l_imag_rotated = l_imag_twisted*cosine - right_imag*sine; - - //////////////////////// Main STFT block - // The 'meat' of the algorithm, the code in this block will most resemble the code from vocalrediso.ny - - // first, apply vocal reduction algorithm only in the right bands - fft_bin >= low_bin && fft_bin <= high_bin ? ( - // get constants for this bin - strength = strength_buffer[fft_bin]; - phase_width = phase_width_buffer[fft_bin]; - - // cacluate energy for this bin - norm_left = sqrt(sqr(l_real_rotated) + sqr(l_imag_rotated)); - norm_right = sqrt(sqr(r_real_rotated) + sqr(r_imag_rotated)); - - // calculate phase difference between left & right, divide by phase_width - w1 = acos((l_real_rotated * r_real_rotated + l_imag_rotated * r_imag_rotated) / (norm_left * norm_right)) / phase_width; - // calculate weight: truncate w1 to [0, 1] and apply strength, then take 1 - the result, and multiply - // by 1 - the square of the difference between the two norm values divided by the sum of the two, moderated by strength * slider10 - weight = (1 - min(1, max(0, w1)) ^ strength) * ( - 1 - sqr(norm_left - norm_right)/sqr(norm_left + norm_right) - ) ^ (slider10 / strength) * 0.5; - - // find the c channel, the sum of the two complex numbers - c_real = l_real_rotated + r_real_rotated; - c_imag = l_imag_rotated + r_imag_rotated; - - // apply weight to c channel - c_real *= weight; - c_imag *= weight; - ) : - ( - c_real = 0; - c_imag = 0; - ); - //////////////////////// END MAIN STFT block - - // apply wet dry mix - out_l_real = l_real_rotated * dry_mix + c_real * wet_mix; - out_l_imag = l_imag_rotated * dry_mix + c_imag * wet_mix; - out_r_real = r_real_rotated * dry_mix + c_real * wet_mix; - out_r_imag = r_imag_rotated * dry_mix + c_imag * wet_mix; - - // output corrections - slider11 > 0.5 ? ( - // if requested, undo input corrections - - // unmix by phase - l_real_out_twisted = out_l_real*cosine - out_r_real*sine; - l_imag_out_twisted = out_l_imag*cosine - out_r_imag*sine; - right_real = out_r_real*cosine + out_l_real*sine; - right_imag = out_r_imag*cosine + out_l_imag*sine; - - left_real = l_real_out_twisted * cosine2 - l_imag_out_twisted * sine2; - left_imag = l_imag_out_twisted * cosine2 + l_real_out_twisted * sine2; - ) : - ( - // else, just copy the values - left_real = out_l_real; - left_imag = out_l_imag; - right_real = out_r_real; - right_imag = out_r_imag; - ); - - // Re-fold back into complex spectrum - fft_buffer[2*fft_bin] = (left_real - right_imag)*0.5; - fft_buffer[2*fft_bin + 1] = (left_imag + right_real)*0.5; - fft_buffer[2*fft_bin2] = (left_real + right_imag)*0.5; - fft_buffer[2*fft_bin2 + 1] = (-left_imag + right_real)*0.5; - - fft_bin += 1; - ); -); - -MAX_FFT_SIZE = 32768; -fft_size = 8192; - -freemem = 0; -fft_buffer = simple_alloc(MAX_FFT_SIZE*2); -window_buffer = simple_alloc(MAX_FFT_SIZE); - -strength_buffer = simple_alloc(MAX_FFT_SIZE); -phase_width_buffer = simple_alloc(MAX_FFT_SIZE); - -buffer_length = srate; -buffer_index = 0; -input_buffer = simple_alloc(buffer_length*2); -output_buffer = simple_alloc(buffer_length*2); - -function window(r) ( - // When squared, the Hann window adds up perfectly for overlap >= 4, so it's suitable for perfect reconstruction - //(0.5 - 0.5*cos(r*2*$pi))/sqrt(0.375); - // the MLT sine window also appears to add up correctly, with sigma = sqrt(2). - sin(r*$pi)*sqrt(2); -); - -overlap_factor = 4; -fft_interval = fft_size/overlap_factor; -fft_scaling_factor = 1/overlap_factor/fft_size; - -fft_size != prev_fft_size ? ( - setup_stft_state(fft_size, prev_fft_size == 0); - prev_fft_size = fft_size; - // Fill window buffer - i = 0; - loop(fft_size, - r = i/fft_size; - window_buffer[i] = window(r); - i += 1; - ); -); - -pdc_delay = fft_size; -pdc_bot_ch = 0; -pdc_top_ch = 2; - -freembuf(freemem); - - -@slider -//IFTEST function slider_code() ( -// convert low cut and high cut to bins every time a slider is changed -low_bin = min(slider5, slider6) / srate * fft_size; -high_bin = max(slider6, slider5) / srate * fft_size; -// convert to radians -rotation = slider7*$pi/180; -// convert percentage to raw scale factor -dry_mix = slider1/100; -wet_mix = slider2/100; -low_strength = slider3; -high_strength = slider4; -phase_width_low = slider8*$pi/180; -phase_width_high = slider9*$pi/180; -cosine = cos(rotation); -sine = sin(rotation); -cosine2 = cos(slider12*$pi/180); -sine2 = sin(slider12*$pi/180); -// fill strength_buffer and phase_width_buffer -band_index = 0; -loop(fft_size, - band_index >= low_bin && band_index <= high_bin ? - ( - // only set values for the appropriate frequency range - frac = (band_index - low_bin)/(high_bin - low_bin - 1); - frac = max(0, min(1, frac)); - // fraction of progress through range [low_bin, high_bin) - strength = low_strength* (1 - frac) + high_strength * frac; - strength_buffer[band_index] = 10^strength; - // precaculate strength (actual value should be positive, so it makes - // sense to take the power of ten, but only after the - // linear mapping over the spectrum is done. - // precalculate phase width - phase_width_buffer[band_index] = phase_width_low * (1 - frac) + phase_width_high * frac; - ); - - band_index += 1; - // next index -); -//IFTEST ); // slider_code() ( - - -@sample -//IFTEST function sample_code() ( -input_buffer[buffer_index*2] = spl0; -input_buffer[buffer_index*2 + 1] = spl1; - -output_index = buffer_index - fft_size; -output_index < 0 ? output_index += buffer_length; -silence > 0 ? ( - spl0 = spl1 = 0; - // silence for fft init -) : ( - spl0 = output_buffer[output_index*2]; - spl1 = output_buffer[output_index*2 + 1]; -); -output_buffer[output_index*2] = 0; // clear the sample we just read -output_buffer[output_index*2 + 1] = 0; - - -fft_counter += 1; -fft_counter >= fft_interval ? ( - fft_counter = 0; - - // Copy input to buffer - bi = buffer_index - fft_size + 1; - i = 0; - loop(fft_size, - i2 = bi + i; - i2 < 0 ? i2 += buffer_length; - - fft_buffer[2*i] = input_buffer[2*i2]*window_buffer[i]; - fft_buffer[2*i + 1] = input_buffer[2*i2 + 1]*window_buffer[i]; - - i += 1; - ); - - // Process buffer - fft(fft_buffer, fft_size); - fft_permute(fft_buffer, fft_size); - - process_stft_segment(fft_buffer, fft_size); - - fft_ipermute(fft_buffer, fft_size); - ifft(fft_buffer, fft_size); - - // Add to output - bi = buffer_index - fft_size + 1; - i = 0; - loop(fft_size, - i2 = bi + i; - (i2 < 0) ? i2 += buffer_length; - - output_buffer[2*i2] += fft_buffer[2*i]*fft_scaling_factor*window_buffer[i]; - output_buffer[2*i2 + 1] += fft_buffer[2*i + 1]*fft_scaling_factor*window_buffer[i]; - - i += 1; - ); - silence > 0 ? silence -= 1; -); - -buffer_index = (buffer_index + 1)%buffer_length; - -//IFTEST ); // function sample_code() ( - - -/*IFTEST{* // main test block - -// helpers -function sum_first_pdc_samples(s0val, s1val) ( - setup_stft_state(fft_size, 1); - spl0sum = 0; - spl1sum = 0; - loop(pdc_delay, - spl0=s0val; - spl1=s1val; - sample_code(); - spl0sum += abs(spl0); - spl1sum += abs(spl1); - ); -); - -//DEBUGPRINT("SETUP: fft range is [0, srate/2] dry=100 wet=100\n"); -slider1 = slider2 = 100; -slider5 = 0; -slider6 = srate/2; -slider_code(); - -//DEBUGPRINT("SAMPLE_0_0_TEST\n"); -// spl0=0 spl1=0 for 100 output samples -sum_first_pdc_samples(0, 0); -assert_equal_exact(0, spl0sum, "spl0sum for init"); -assert_equal_exact(0, spl1sum, "spl1sum for init"); - -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"); -//DEBUGPRINT("spl0sum=%g, spl1sum=%g\n", spl0sum, spl1sum); - -//DEBUGPRINT("SAMPLE_1_0_TEST\n"); -// spl0=1 spl1=0 for for 100 output samples -sum_first_pdc_samples(1, 0); -assert_equal_exact(0, spl0sum, "spl0sum for init"); -assert_equal_exact(0, spl1sum, "spl1sum for init"); - -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"); -//DEBUGPRINT("spl0sum=%g, spl1sum=%g\n", spl0sum, spl1sum); - -//DEBUGPRINT("SAMPLE_0_1_TEST\n"); -// spl0=0 spl1=1 for for 100 output samples -sum_first_pdc_samples(0, 1); -assert_equal_exact(0, spl0sum, "spl0sum for init"); -assert_equal_exact(0, spl1sum, "spl1sum for init"); - -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"); -//DEBUGPRINT("spl0sum=%g, spl1sum=%g\n", spl0sum, spl1sum); - -//DEBUGPRINT("SAMPLE_1_1_TEST\n"); -// spl0=1 spl1=1 for for 100 output samples -sum_first_pdc_samples(1, 1); -assert_equal_exact(0, spl0sum, "spl0sum for init"); -assert_equal_exact(0, spl1sum, "spl1sum for init"); - -loop(100, - spl0=1; - spl1=1; - sample_code(); - spl0sum += abs(spl0); - spl1sum += abs(spl1); -); -assert_near_equal(200, 0.001, spl0sum, "SAMPLE_1_1_TEST: spl0sum was not as expected"); -assert_near_equal(200, 0.001, spl1sum, "SAMPLE_1_1_TEST: spl1sum was not as expected"); -//DEBUGPRINT("spl0sum=%g, spl1sum=%g\n", spl0sum, spl1sum); - - -//DEBUGPRINT("SETUP: fft range is [srate/4, srate/2] dry=100 wet=100\n"); -// check no values before first output sample -slider1 = slider2 = 100; -slider5 = srate/4; -slider6 = srate/2; -//DEBUGPRINT("PDC_SILENCE_ON_PARTIAL_RANGE_TEST\n"); -slider_code(); -sum_first_pdc_samples(1, 1); -assert_equal_exact(0, spl0sum, "PDC_SILENCE_ON_PARTIAL_RANGE_TEST failed on spl0"); -assert_equal_exact(0, spl1sum, "PDC_SILENCE_ON_PARTIAL_RANGE_TEST failed on spl1"); - -test_summary(); -//*}IFTEST*/ // main test block +// Based on Vocalrediso.ny, a nyquist filter for audacity, currently released under the GNU GPL V3 +// see https://www.gnu.org/licenses/gpl-3.0.en.html, https://www.audacityteam.org/ +// credits to Neil Bickford for his denoiser, https://github.com/Nbickford/REAPERDenoiser, +// which was used as a starting point for this script +// which now uses the STFT template code by geraintluff: https://forum.cockos.com/showthread.php?t=225955 + +//IFTEST /* //this section not used by testing code +desc: vocal removal/isolation +//tags: processing vocals stereo +//author: Michael Pannekoek + +//slider1:fft_size_index=2<0,7,1{256,512,1024,2048,4096,8192,16384,32768}>FFT size +slider1:0<-100,100,0.1>dry mix +slider2:0<-100,100,0.1>C mix (Vocals) +slider3:0<-5,5,0.001>strength at Low Cut +slider4:0<-5,5,0.001>strength at High Cut +slider5:80<0,24000,10>Low Cut (Vocals) +slider6:24000<0,24000,10>High Cut (Vocals) +slider7:0<-90,90,0.1>Phase (Degrees) +slider8:0<-5,5,0.001>Phase width at Low Cut +slider9:0<-5,5,0.001>Phase width at High Cut +slider10:1<0,1,0.05>Attenuate if different volume +slider11:1<0,1,1{No,Yes}>undo input corrections +slider12:0<-180,180,0.05>Phase2 (Degrees) +//IFTEST */ //this section not used by testing code + +/*IFJSFX{ +in_pin:left input +in_pin:right input +out_pin:left output +out_pin:right output +}IFJSFX*/ + + + +//IFTEST /* +@init +//IFTEST */ +//IFJSFX /* +slider1=0; +slider2=0; +slider3=0; +slider4=0; +slider5=80; +slider6=24000; +slider7=0; +slider8=0; +slider9=0; +slider10=1; +slider11=1; +slider12=0; +//IFJSFX */ + +//IFTEST srate = 24000; + +_free_memory = 0; +function simple_alloc(amount) +local(_free_memory_old) +global(_free_memory) +( + _free_memory_old = _free_memory; + _free_memory += amount; + _free_memory_old; +); + +prev_fft_size = 0; + +function setup_stft_state(fft_size, first_time) ( + //////////////////////// Setup block + // This is called when playback starts, or when the FFT size is changed + memset(output_buffer, 0, buffer_length*2); + memset(input_buffer, 0, buffer_length*2); + // reset indexes and counters + buffer_index = 0; + output_index = 0; + fft_counter = 0; + // force silence for first overlaps where fft not yet run + silence = overlap_factor; + //////////////////////// +); + +function get_weight(strength, phase_width, l_real, l_imag, r_real, r_imag) local(norm_left, norm_right, w1, weight) ( + // cacluate energy for this bin + norm_left = sqrt(sqr(l_real) + sqr(l_imag)); + norm_right = sqrt(sqr(r_real) + sqr(r_imag)); + + // calculate phase difference between left & right, divide by phase_width + w1 = (l_real * r_real + l_imag * r_imag) / (norm_left * norm_right * phase_width * 2); + weight = ((max(0, min(1, w1+0.5))) ^ strength); + // calculate weight: truncate w1 to [0, 1] and apply strength, then take 1 - the result, and multiply + // by 1 - the square of the difference between the two norm values divided by the sum of the two, moderated by strength * slider10 + slider10 > 0 ? weight *= ( + 1 - sqr(norm_left - norm_right)/sqr(norm_left + norm_right) + ) ^ (slider10 / strength) * 0.5; + weight; +); + +setup_stft_state(fft_size, prev_fft_size == 0); + +function process_stft_segment(fft_buffer, fft_size) local(fft_bin, left_real, left_imag, right_real, right_imag) ( + fft_bin = 0; // FFT bin number + loop(fft_size/2+1, + fft_bin2 = fft_bin ? (fft_size - fft_bin) : 0; + + // Unfold complex spectrum into two real spectra + left_real = fft_buffer[2*fft_bin] + fft_buffer[2*fft_bin2]; + left_imag = fft_buffer[2*fft_bin + 1] - fft_buffer[2*fft_bin2 + 1]; + right_real = fft_buffer[2*fft_bin + 1] + fft_buffer[2*fft_bin2 + 1]; + right_imag = -fft_buffer[2*fft_bin] + fft_buffer[2*fft_bin2]; + + // input corrections + // first, change the phase of L based on phase2: + l_real_twisted = left_real*cosine2 + left_imag*sine2; + l_imag_twisted = left_imag*cosine2 - left_real*sine2; + + // now mix L&R together based on phase + r_real_rotated = right_real*cosine + l_real_twisted*sine; + r_imag_rotated = right_imag*cosine + l_imag_twisted*sine; + l_real_rotated = l_real_twisted*cosine - right_real*sine; + l_imag_rotated = l_imag_twisted*cosine - right_imag*sine; + + //////////////////////// Main STFT block + // The 'meat' of the algorithm, the code in this block will most resemble the code from vocalrediso.ny + + weight = 0; + // first, apply vocal reduction algorithm only in the right bands + fft_bin >= low_bin && fft_bin <= high_bin ? ( + weight = get_weight(strength_buffer[fft_bin], phase_width_buffer[fft_bin], l_real_rotated, l_imag_rotated, r_real_rotated, r_imag_rotated); + // find the c channel, the sum of the two complex numbers + c_real = l_real_rotated + r_real_rotated; + c_imag = l_imag_rotated + r_imag_rotated; + + // apply weight to c channel + c_real *= weight; + c_imag *= weight; + ) : + ( + // let wet signal have 0 for fft coefficients when out of bounds + c_real = 0; + c_imag = 0; + ); + //////////////////////// END MAIN STFT block + + // apply wet dry mix + out_l_real = l_real_rotated * dry_mix + c_real * wet_mix; + out_l_imag = l_imag_rotated * dry_mix + c_imag * wet_mix; + out_r_real = r_real_rotated * dry_mix + c_real * wet_mix; + out_r_imag = r_imag_rotated * dry_mix + c_imag * wet_mix; + + // output corrections + slider11 > 0.5 ? ( + // if requested, undo input corrections + + // unmix by phase + l_real_out_twisted = out_l_real*cosine - out_r_real*sine; + l_imag_out_twisted = out_l_imag*cosine - out_r_imag*sine; + right_real = out_r_real*cosine + out_l_real*sine; + right_imag = out_r_imag*cosine + out_l_imag*sine; + + left_real = l_real_out_twisted * cosine2 - l_imag_out_twisted * sine2; + left_imag = l_imag_out_twisted * cosine2 + l_real_out_twisted * sine2; + ) : + ( + // else, just copy the values + left_real = out_l_real; + left_imag = out_l_imag; + right_real = out_r_real; + right_imag = out_r_imag; + ); + + // Re-fold back into complex spectrum + fft_buffer[2*fft_bin] = (left_real - right_imag)*0.5; + fft_buffer[2*fft_bin + 1] = (left_imag + right_real)*0.5; + fft_buffer[2*fft_bin2] = (left_real + right_imag)*0.5; + fft_buffer[2*fft_bin2 + 1] = (-left_imag + right_real)*0.5; + + fft_bin += 1; + ); +); + +MAX_FFT_SIZE = 32768; +fft_size = 8192; + +freemem = 0; +fft_buffer = simple_alloc(MAX_FFT_SIZE*2); +window_buffer = simple_alloc(MAX_FFT_SIZE); + +strength_buffer = simple_alloc(MAX_FFT_SIZE); +phase_width_buffer = simple_alloc(MAX_FFT_SIZE); + +buffer_length = srate; +buffer_index = 0; +input_buffer = simple_alloc(buffer_length*2); +output_buffer = simple_alloc(buffer_length*2); + +function window(r) ( + // When squared, the Hann window adds up perfectly for overlap >= 4, so it's suitable for perfect reconstruction + //(0.5 - 0.5*cos(r*2*$pi))/sqrt(0.375); + // the MLT sine window also appears to add up correctly, with sigma = sqrt(2). + sin(r*$pi)*sqrt(2); +); + +overlap_factor = 4; +fft_interval = fft_size/overlap_factor; +fft_scaling_factor = 1/overlap_factor/fft_size; + +fft_size != prev_fft_size ? ( + setup_stft_state(fft_size, prev_fft_size == 0); + prev_fft_size = fft_size; + // Fill window buffer + i = 0; + loop(fft_size, + r = i/fft_size; + window_buffer[i] = window(r); + i += 1; + ); +); + +pdc_delay = fft_size; +pdc_bot_ch = 0; +pdc_top_ch = 2; + +/*IFJSFX{ +freembuf(freemem); + +@slider +}IFJSFX*/ +//IFTEST function slider_code() ( +// convert low cut and high cut to bins every time a slider is changed +low_bin = min(slider5, slider6) / srate * fft_size; +high_bin = max(slider6, slider5) / srate * fft_size; +// convert to radians +rotation = slider7*$pi/180; +// convert percentage to raw scale factor +dry_mix = slider1/100; +wet_mix = slider2/100; +low_strength = slider3; +high_strength = slider4; +phase_width_low = slider8; +phase_width_high = slider9; +cosine = cos(rotation); +sine = sin(rotation); +cosine2 = cos(slider12*$pi/180); +sine2 = sin(slider12*$pi/180); +// fill strength_buffer and phase_width_buffer +band_index = 0; +loop(fft_size/2+1, + band_index >= low_bin && band_index <= high_bin ? + ( + // only set values for the appropriate frequency range + frac = (band_index - low_bin)/(high_bin - low_bin + 1); + // fraction of progress through range [low_bin, high_bin) + strength = low_strength* (1 - frac) + high_strength * frac; + strength_buffer[band_index] = exp(strength); + // precaculate strength (actual value should be positive, so it makes + // sense to take the power of ten, but only after the + // linear mapping over the spectrum is done. + // precalculate phase width + phase_width = phase_width_low * (1 - frac) + phase_width_high * frac; + phase_width_buffer[band_index] = exp(phase_width); + ); + + band_index += 1; + // next index +); +//IFTEST ); // slider_code() ( + +//IFTEST /* +@sample +//IFTEST */ +//IFTEST function sample_code() ( +input_buffer[buffer_index*2] = spl0; +input_buffer[buffer_index*2 + 1] = spl1; + +output_index = buffer_index - fft_size; +output_index < 0 ? output_index += buffer_length; +silence > 0 ? ( + spl0 = spl1 = 0; + // silence for fft init +) : ( + spl0 = output_buffer[output_index*2]; + spl1 = output_buffer[output_index*2 + 1]; +); +output_buffer[output_index*2] = 0; // clear the sample we just read +output_buffer[output_index*2 + 1] = 0; + +fft_counter += 1; +fft_counter >= fft_interval ? ( + fft_counter = 0; + + // Copy input to buffer + bi = buffer_index - fft_size + 1; + i = 0; + loop(fft_size, + i2 = bi + i; + i2 < 0 ? i2 += buffer_length; + + fft_buffer[2*i] = input_buffer[2*i2]*window_buffer[i]; + fft_buffer[2*i + 1] = input_buffer[2*i2 + 1]*window_buffer[i]; + + i += 1; + ); + + // Process buffer + fft(fft_buffer, fft_size); + fft_permute(fft_buffer, fft_size); + + process_stft_segment(fft_buffer, fft_size); + + fft_ipermute(fft_buffer, fft_size); + ifft(fft_buffer, fft_size); + + // Add to output + bi = buffer_index - fft_size + 1; + i = 0; + loop(fft_size, + i2 = bi + i; + (i2 < 0) ? i2 += buffer_length; + + output_buffer[2*i2] += fft_buffer[2*i]*fft_scaling_factor*window_buffer[i]; + output_buffer[2*i2 + 1] += fft_buffer[2*i + 1]*fft_scaling_factor*window_buffer[i]; + + i += 1; + ); + silence > 0 ? silence -= 1; +); + +buffer_index = (buffer_index + 1)%buffer_length; + +//IFTEST ); // function sample_code() + + +/*IFJSFX{ +@serialize +}IFJSFX*/ + /* +//IFJSFX /* +function file_var(file, val) ( + printf("FILE_VAR, FILE: %d, VAL: %g\n", file, val) +); +//IFJSFX */ +serial_version = 1.00; +file_var(0, serial_version); + */ +// nothing serialized yet, but keep track of the serial_version +// for the preset state of the original plugin, serial_version would now be euqal to 0. + +/*IFTEST{ // main test block + +// helpers +function sum_first_pdc_samples(s0val, s1val) ( + setup_stft_state(fft_size, 1); + spl0sum = 0; + spl1sum = 0; + loop(pdc_delay, + spl0=s0val; + spl1=s1val; + sample_code(); + spl0sum += abs(spl0); + spl1sum += abs(spl1); + ); +); + +printf("SETUP: fft range is [0, srate/2] dry=100 wet=100 attenuate_diff_vol=1\n"); +slider1 = slider2 = 100; +slider5 = 0; +slider6 = srate/2; +slider10 = 1; +slider_code(); +printf("SETUP: got low_bin=%g, high_bin=%g\n", low_bin, high_bin); + +printf("get_weight_ALL_ZERO\n"); +assert_equal_exact(0, get_weight(1, 1, 0, 0, 0, 0)); + +printf("get_weight_PERPENDICULAR\n"); +assert_equal_exact(0.25, get_weight(1, 1, 0, 1, 1, 0)); + +printf("get_weight_EQUIVALENT\n"); +assert_equal_exact(0.5, get_weight(1, 1, 1, 0, 1, 0)); + +printf("get_weight_INVERSE\n"); +assert_equal_exact(0, get_weight(1, 1, 1, 0, -1, 0)); + +printf("get_weight_ONE_ONLY\n"); +assert_equal_exact(0, get_weight(1, 1, 1, 0, 0, 0)); + +printf("get_weight_HALF\n"); +assert_equal_exact(0.5, get_weight(1, 1, 0.5, 0, 0.5, 0)); + +printf("get_weight_STRANGE\n"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -33.28894848631978, -662.24573314730117, -33.28894848631978, -662.24573314730117), "STRANGE01"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, 452.28715419177155, 53.67874245364300, 452.28715419177155, 53.67874245364300), "STRANGE02"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, 267.79446776740104, 17.14616043377475, 267.79446776740104, 17.14616043377475), "STRANGE03"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -32.55549374859478, 261.41171292127513, -32.55549374859478, 261.41171292127513), "STRANGE04"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -2.75689966923304, -217.67592724424455, -2.75689966923304, -217.67592724424455), "STRANGE05"); + +printf("SAMPLE_0_0_TEST\n"); +// spl0=0 spl1=0 for 100 output samples +sum_first_pdc_samples(0, 0); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +loop(100, + spl0=1; + spl1=1; + sample_code(); + spl0sum += abs(spl0); + spl1sum += abs(spl1); +); +assert_near_equal(200, 0.001, spl0sum, "SAMPLE_1_1_TEST: spl0sum was not as expected"); +assert_near_equal(200, 0.001, spl1sum, "SAMPLE_1_1_TEST: spl1sum was not as expected"); +printf("spl0sum=%g, spl1sum=%g\n", spl0sum, spl1sum); + + +printf("SETUP: fft range is [srate/4, srate/2] dry=100 wet=100\n"); +// check no values before first output sample +slider1 = slider2 = 100; +slider5 = srate/4; +slider6 = srate/2; +printf("PDC_SILENCE_ON_PARTIAL_RANGE_TEST\n"); +slider_code(); +sum_first_pdc_samples(1, 1); +assert_equal_exact(0, spl0sum, "PDC_SILENCE_ON_PARTIAL_RANGE_TEST failed on spl0"); +assert_equal_exact(0, spl1sum, "PDC_SILENCE_ON_PARTIAL_RANGE_TEST failed on spl1"); + +test_summary(); +}IFTEST*/ // main test block diff --git a/vocalrediso.jsfx b/vocalrediso.jsfx index 1d2e118..ab619a6 100644 --- a/vocalrediso.jsfx +++ b/vocalrediso.jsfx @@ -4,6 +4,7 @@ // which was used as a starting point for this script // which now uses the STFT template code by geraintluff: https://forum.cockos.com/showthread.php?t=225955 +//IFTEST /* //this section not used by testing code desc: vocal removal/isolation //tags: processing vocals stereo //author: Michael Pannekoek @@ -16,11 +17,13 @@ slider4:0<-5,5,0.001>strength at High Cut slider5:80<0,24000,10>Low Cut (Vocals) slider6:24000<0,24000,10>High Cut (Vocals) slider7:0<-90,90,0.1>Phase (Degrees) -slider8:90<1,180,0.1>Phase width at Low Cut (Degrees) -slider9:90<1,180,0.1>Phase width at High Cut (Degrees) +slider8:0<-5,5,0.001>Phase width at Low Cut +slider9:0<-5,5,0.001>Phase width at High Cut slider10:1<0,1,0.05>Attenuate if different volume slider11:1<0,1,1{No,Yes}>undo input corrections slider12:0<-180,180,0.05>Phase2 (Degrees) +//IFTEST */ //this section not used by testing code + in_pin:left input in_pin:right input @@ -28,7 +31,28 @@ out_pin:left output out_pin:right output + + +//IFTEST /* @init +//IFTEST */ + /* +slider1=0; +slider2=0; +slider3=0; +slider4=0; +slider5=80; +slider6=24000; +slider7=0; +slider8=0; +slider9=0; +slider10=1; +slider11=1; +slider12=0; + */ + +//IFTEST srate = 24000; + _free_memory = 0; function simple_alloc(amount) local(_free_memory_old) @@ -44,7 +68,6 @@ prev_fft_size = 0; function setup_stft_state(fft_size, first_time) ( //////////////////////// Setup block // This is called when playback starts, or when the FFT size is changed - //memset(fft_buffer, 0, MAX_FFT_SIZE*2); memset(output_buffer, 0, buffer_length*2); memset(input_buffer, 0, buffer_length*2); // reset indexes and counters @@ -56,6 +79,22 @@ function setup_stft_state(fft_size, first_time) ( //////////////////////// ); +function get_weight(strength, phase_width, l_real, l_imag, r_real, r_imag) local(norm_left, norm_right, w1, weight) ( + // cacluate energy for this bin + norm_left = sqrt(sqr(l_real) + sqr(l_imag)); + norm_right = sqrt(sqr(r_real) + sqr(r_imag)); + + // calculate phase difference between left & right, divide by phase_width + w1 = (l_real * r_real + l_imag * r_imag) / (norm_left * norm_right * phase_width * 2); + weight = ((max(0, min(1, w1+0.5))) ^ strength); + // calculate weight: truncate w1 to [0, 1] and apply strength, then take 1 - the result, and multiply + // by 1 - the square of the difference between the two norm values divided by the sum of the two, moderated by strength * slider10 + slider10 > 0 ? weight *= ( + 1 - sqr(norm_left - norm_right)/sqr(norm_left + norm_right) + ) ^ (slider10 / strength) * 0.5; + weight; +); + setup_stft_state(fft_size, prev_fft_size == 0); function process_stft_segment(fft_buffer, fft_size) local(fft_bin, left_real, left_imag, right_real, right_imag) ( @@ -83,24 +122,10 @@ function process_stft_segment(fft_buffer, fft_size) local(fft_bin, left_real, le //////////////////////// Main STFT block // The 'meat' of the algorithm, the code in this block will most resemble the code from vocalrediso.ny + weight = 0; // first, apply vocal reduction algorithm only in the right bands fft_bin >= low_bin && fft_bin <= high_bin ? ( - // get constants for this bin - strength = strength_buffer[fft_bin]; - phase_width = phase_width_buffer[fft_bin]; - - // cacluate energy for this bin - norm_left = sqrt(sqr(l_real_rotated) + sqr(l_imag_rotated)); - norm_right = sqrt(sqr(r_real_rotated) + sqr(r_imag_rotated)); - - // calculate phase difference between left & right, divide by phase_width - w1 = acos((l_real_rotated * r_real_rotated + l_imag_rotated * r_imag_rotated) / (norm_left * norm_right)) / phase_width; - // calculate weight: truncate w1 to [0, 1] and apply strength, then take 1 - the result, and multiply - // by 1 - the square of the difference between the two norm values divided by the sum of the two, moderated by strength * slider10 - weight = (1 - min(1, max(0, w1)) ^ strength) * ( - 1 - sqr(norm_left - norm_right)/sqr(norm_left + norm_right) - ) ^ (slider10 / strength) * 0.5; - + weight = get_weight(strength_buffer[fft_bin], phase_width_buffer[fft_bin], l_real_rotated, l_imag_rotated, r_real_rotated, r_imag_rotated); // find the c channel, the sum of the two complex numbers c_real = l_real_rotated + r_real_rotated; c_imag = l_imag_rotated + r_imag_rotated; @@ -110,6 +135,7 @@ function process_stft_segment(fft_buffer, fft_size) local(fft_bin, left_real, le c_imag *= weight; ) : ( + // let wet signal have 0 for fft coefficients when out of bounds c_real = 0; c_imag = 0; ); @@ -194,10 +220,12 @@ pdc_delay = fft_size; pdc_bot_ch = 0; pdc_top_ch = 2; -freembuf(freemem); +freembuf(freemem); @slider + +//IFTEST function slider_code() ( // convert low cut and high cut to bins every time a slider is changed low_bin = min(slider5, slider6) / srate * fft_size; high_bin = max(slider6, slider5) / srate * fft_size; @@ -208,36 +236,39 @@ dry_mix = slider1/100; wet_mix = slider2/100; low_strength = slider3; high_strength = slider4; -phase_width_low = slider8*$pi/180; -phase_width_high = slider9*$pi/180; +phase_width_low = slider8; +phase_width_high = slider9; cosine = cos(rotation); sine = sin(rotation); cosine2 = cos(slider12*$pi/180); sine2 = sin(slider12*$pi/180); // fill strength_buffer and phase_width_buffer band_index = 0; -loop(fft_size, +loop(fft_size/2+1, band_index >= low_bin && band_index <= high_bin ? ( // only set values for the appropriate frequency range - frac = (band_index - low_bin)/(high_bin - low_bin - 1); - frac = max(0, min(1, frac)); + frac = (band_index - low_bin)/(high_bin - low_bin + 1); // fraction of progress through range [low_bin, high_bin) strength = low_strength* (1 - frac) + high_strength * frac; - strength_buffer[band_index] = 10^strength; + strength_buffer[band_index] = exp(strength); // precaculate strength (actual value should be positive, so it makes // sense to take the power of ten, but only after the // linear mapping over the spectrum is done. // precalculate phase width - phase_width_buffer[band_index] = phase_width_low * (1 - frac) + phase_width_high * frac; + phase_width = phase_width_low * (1 - frac) + phase_width_high * frac; + phase_width_buffer[band_index] = exp(phase_width); ); band_index += 1; // next index ); +//IFTEST ); // slider_code() ( - +//IFTEST /* @sample +//IFTEST */ +//IFTEST function sample_code() ( input_buffer[buffer_index*2] = spl0; input_buffer[buffer_index*2 + 1] = spl1; @@ -253,7 +284,6 @@ silence > 0 ? ( output_buffer[output_index*2] = 0; // clear the sample we just read output_buffer[output_index*2 + 1] = 0; - fft_counter += 1; fft_counter >= fft_interval ? ( fft_counter = 0; @@ -292,12 +322,157 @@ fft_counter >= fft_interval ? ( i += 1; ); + silence > 0 ? silence -= 1; ); buffer_index = (buffer_index + 1)%buffer_length; +//IFTEST ); // function sample_code() + + + @serialize + +//IFEEL /* + /* +function file_var(file, val) ( + printf("FILE_VAR, FILE: %d, VAL: %g\n", file, val) +); + */ serial_version = 1.00; file_var(0, serial_version); +//IFEEL */ // nothing serialized yet, but keep track of the serial_version // for the preset state of the original plugin, serial_version would now be euqal to 0. + +/*IFTEST{ // main test block + +// helpers +function sum_first_pdc_samples(s0val, s1val) ( + setup_stft_state(fft_size, 1); + spl0sum = 0; + spl1sum = 0; + loop(pdc_delay, + spl0=s0val; + spl1=s1val; + sample_code(); + spl0sum += abs(spl0); + spl1sum += abs(spl1); + ); +); + +printf("SETUP: fft range is [0, srate/2] dry=100 wet=100 attenuate_diff_vol=1\n"); +slider1 = slider2 = 100; +slider5 = 0; +slider6 = srate/2; +slider10 = 1; +slider_code(); +printf("SETUP: got low_bin=%g, high_bin=%g\n", low_bin, high_bin); + +printf("get_weight_ALL_ZERO\n"); +assert_equal_exact(0, get_weight(1, 1, 0, 0, 0, 0)); + +printf("get_weight_PERPENDICULAR\n"); +assert_equal_exact(0.25, get_weight(1, 1, 0, 1, 1, 0)); + +printf("get_weight_EQUIVALENT\n"); +assert_equal_exact(0.5, get_weight(1, 1, 1, 0, 1, 0)); + +printf("get_weight_INVERSE\n"); +assert_equal_exact(0, get_weight(1, 1, 1, 0, -1, 0)); + +printf("get_weight_ONE_ONLY\n"); +assert_equal_exact(0, get_weight(1, 1, 1, 0, 0, 0)); + +printf("get_weight_HALF\n"); +assert_equal_exact(0.5, get_weight(1, 1, 0.5, 0, 0.5, 0)); + +printf("get_weight_STRANGE\n"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -33.28894848631978, -662.24573314730117, -33.28894848631978, -662.24573314730117), "STRANGE01"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, 452.28715419177155, 53.67874245364300, 452.28715419177155, 53.67874245364300), "STRANGE02"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, 267.79446776740104, 17.14616043377475, 267.79446776740104, 17.14616043377475), "STRANGE03"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -32.55549374859478, 261.41171292127513, -32.55549374859478, 261.41171292127513), "STRANGE04"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -2.75689966923304, -217.67592724424455, -2.75689966923304, -217.67592724424455), "STRANGE05"); + +printf("SAMPLE_0_0_TEST\n"); +// spl0=0 spl1=0 for 100 output samples +sum_first_pdc_samples(0, 0); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +loop(100, + spl0=1; + spl1=1; + sample_code(); + spl0sum += abs(spl0); + spl1sum += abs(spl1); +); +assert_near_equal(200, 0.001, spl0sum, "SAMPLE_1_1_TEST: spl0sum was not as expected"); +assert_near_equal(200, 0.001, spl1sum, "SAMPLE_1_1_TEST: spl1sum was not as expected"); +printf("spl0sum=%g, spl1sum=%g\n", spl0sum, spl1sum); + + +printf("SETUP: fft range is [srate/4, srate/2] dry=100 wet=100\n"); +// check no values before first output sample +slider1 = slider2 = 100; +slider5 = srate/4; +slider6 = srate/2; +printf("PDC_SILENCE_ON_PARTIAL_RANGE_TEST\n"); +slider_code(); +sum_first_pdc_samples(1, 1); +assert_equal_exact(0, spl0sum, "PDC_SILENCE_ON_PARTIAL_RANGE_TEST failed on spl0"); +assert_equal_exact(0, spl1sum, "PDC_SILENCE_ON_PARTIAL_RANGE_TEST failed on spl1"); + +test_summary(); +}IFTEST*/ // main test block diff --git a/vocalrediso.jsfx-template b/vocalrediso.jsfx-template new file mode 100644 index 0000000..5d32519 --- /dev/null +++ b/vocalrediso.jsfx-template @@ -0,0 +1,478 @@ +// Based on Vocalrediso.ny, a nyquist filter for audacity, currently released under the GNU GPL V3 +// see https://www.gnu.org/licenses/gpl-3.0.en.html, https://www.audacityteam.org/ +// credits to Neil Bickford for his denoiser, https://github.com/Nbickford/REAPERDenoiser, +// which was used as a starting point for this script +// which now uses the STFT template code by geraintluff: https://forum.cockos.com/showthread.php?t=225955 + +//IFTEST /* //this section not used by testing code +desc: vocal removal/isolation +//tags: processing vocals stereo +//author: Michael Pannekoek + +//slider1:fft_size_index=2<0,7,1{256,512,1024,2048,4096,8192,16384,32768}>FFT size +slider1:0<-100,100,0.1>dry mix +slider2:0<-100,100,0.1>C mix (Vocals) +slider3:0<-5,5,0.001>strength at Low Cut +slider4:0<-5,5,0.001>strength at High Cut +slider5:80<0,24000,10>Low Cut (Vocals) +slider6:24000<0,24000,10>High Cut (Vocals) +slider7:0<-90,90,0.1>Phase (Degrees) +slider8:0<-5,5,0.001>Phase width at Low Cut +slider9:0<-5,5,0.001>Phase width at High Cut +slider10:1<0,1,0.05>Attenuate if different volume +slider11:1<0,1,1{No,Yes}>undo input corrections +slider12:0<-180,180,0.05>Phase2 (Degrees) +//IFTEST */ //this section not used by testing code + +/*IFJSFX{ +in_pin:left input +in_pin:right input +out_pin:left output +out_pin:right output +}IFJSFX*/ + + + +//IFTEST /* +@init +//IFTEST */ +//IFJSFX /* +slider1=0; +slider2=0; +slider3=0; +slider4=0; +slider5=80; +slider6=24000; +slider7=0; +slider8=0; +slider9=0; +slider10=1; +slider11=1; +slider12=0; +//IFJSFX */ + +//IFTEST srate = 24000; + +_free_memory = 0; +function simple_alloc(amount) +local(_free_memory_old) +global(_free_memory) +( + _free_memory_old = _free_memory; + _free_memory += amount; + _free_memory_old; +); + +prev_fft_size = 0; + +function setup_stft_state(fft_size, first_time) ( + //////////////////////// Setup block + // This is called when playback starts, or when the FFT size is changed + memset(output_buffer, 0, buffer_length*2); + memset(input_buffer, 0, buffer_length*2); + // reset indexes and counters + buffer_index = 0; + output_index = 0; + fft_counter = 0; + // force silence for first overlaps where fft not yet run + silence = overlap_factor; + //////////////////////// +); + +function get_weight(strength, phase_width, l_real, l_imag, r_real, r_imag) local(norm_left, norm_right, w1, weight) ( + // cacluate energy for this bin + norm_left = sqrt(sqr(l_real) + sqr(l_imag)); + norm_right = sqrt(sqr(r_real) + sqr(r_imag)); + + // calculate phase difference between left & right, divide by phase_width + w1 = (l_real * r_real + l_imag * r_imag) / (norm_left * norm_right * phase_width * 2); + weight = ((max(0, min(1, w1+0.5))) ^ strength); + // calculate weight: truncate w1 to [0, 1] and apply strength, then take 1 - the result, and multiply + // by 1 - the square of the difference between the two norm values divided by the sum of the two, moderated by strength * slider10 + slider10 > 0 ? weight *= ( + 1 - sqr(norm_left - norm_right)/sqr(norm_left + norm_right) + ) ^ (slider10 / strength) * 0.5; + weight; +); + +setup_stft_state(fft_size, prev_fft_size == 0); + +function process_stft_segment(fft_buffer, fft_size) local(fft_bin, left_real, left_imag, right_real, right_imag) ( + fft_bin = 0; // FFT bin number + loop(fft_size/2+1, + fft_bin2 = fft_bin ? (fft_size - fft_bin) : 0; + + // Unfold complex spectrum into two real spectra + left_real = fft_buffer[2*fft_bin] + fft_buffer[2*fft_bin2]; + left_imag = fft_buffer[2*fft_bin + 1] - fft_buffer[2*fft_bin2 + 1]; + right_real = fft_buffer[2*fft_bin + 1] + fft_buffer[2*fft_bin2 + 1]; + right_imag = -fft_buffer[2*fft_bin] + fft_buffer[2*fft_bin2]; + + // input corrections + // first, change the phase of L based on phase2: + l_real_twisted = left_real*cosine2 + left_imag*sine2; + l_imag_twisted = left_imag*cosine2 - left_real*sine2; + + // now mix L&R together based on phase + r_real_rotated = right_real*cosine + l_real_twisted*sine; + r_imag_rotated = right_imag*cosine + l_imag_twisted*sine; + l_real_rotated = l_real_twisted*cosine - right_real*sine; + l_imag_rotated = l_imag_twisted*cosine - right_imag*sine; + + //////////////////////// Main STFT block + // The 'meat' of the algorithm, the code in this block will most resemble the code from vocalrediso.ny + + weight = 0; + // first, apply vocal reduction algorithm only in the right bands + fft_bin >= low_bin && fft_bin <= high_bin ? ( + weight = get_weight(strength_buffer[fft_bin], phase_width_buffer[fft_bin], l_real_rotated, l_imag_rotated, r_real_rotated, r_imag_rotated); + // find the c channel, the sum of the two complex numbers + c_real = l_real_rotated + r_real_rotated; + c_imag = l_imag_rotated + r_imag_rotated; + + // apply weight to c channel + c_real *= weight; + c_imag *= weight; + ) : + ( + // let wet signal have 0 for fft coefficients when out of bounds + c_real = 0; + c_imag = 0; + ); + //////////////////////// END MAIN STFT block + + // apply wet dry mix + out_l_real = l_real_rotated * dry_mix + c_real * wet_mix; + out_l_imag = l_imag_rotated * dry_mix + c_imag * wet_mix; + out_r_real = r_real_rotated * dry_mix + c_real * wet_mix; + out_r_imag = r_imag_rotated * dry_mix + c_imag * wet_mix; + + // output corrections + slider11 > 0.5 ? ( + // if requested, undo input corrections + + // unmix by phase + l_real_out_twisted = out_l_real*cosine - out_r_real*sine; + l_imag_out_twisted = out_l_imag*cosine - out_r_imag*sine; + right_real = out_r_real*cosine + out_l_real*sine; + right_imag = out_r_imag*cosine + out_l_imag*sine; + + left_real = l_real_out_twisted * cosine2 - l_imag_out_twisted * sine2; + left_imag = l_imag_out_twisted * cosine2 + l_real_out_twisted * sine2; + ) : + ( + // else, just copy the values + left_real = out_l_real; + left_imag = out_l_imag; + right_real = out_r_real; + right_imag = out_r_imag; + ); + + // Re-fold back into complex spectrum + fft_buffer[2*fft_bin] = (left_real - right_imag)*0.5; + fft_buffer[2*fft_bin + 1] = (left_imag + right_real)*0.5; + fft_buffer[2*fft_bin2] = (left_real + right_imag)*0.5; + fft_buffer[2*fft_bin2 + 1] = (-left_imag + right_real)*0.5; + + fft_bin += 1; + ); +); + +MAX_FFT_SIZE = 32768; +fft_size = 8192; + +freemem = 0; +fft_buffer = simple_alloc(MAX_FFT_SIZE*2); +window_buffer = simple_alloc(MAX_FFT_SIZE); + +strength_buffer = simple_alloc(MAX_FFT_SIZE); +phase_width_buffer = simple_alloc(MAX_FFT_SIZE); + +buffer_length = srate; +buffer_index = 0; +input_buffer = simple_alloc(buffer_length*2); +output_buffer = simple_alloc(buffer_length*2); + +function window(r) ( + // When squared, the Hann window adds up perfectly for overlap >= 4, so it's suitable for perfect reconstruction + //(0.5 - 0.5*cos(r*2*$pi))/sqrt(0.375); + // the MLT sine window also appears to add up correctly, with sigma = sqrt(2). + sin(r*$pi)*sqrt(2); +); + +overlap_factor = 4; +fft_interval = fft_size/overlap_factor; +fft_scaling_factor = 1/overlap_factor/fft_size; + +fft_size != prev_fft_size ? ( + setup_stft_state(fft_size, prev_fft_size == 0); + prev_fft_size = fft_size; + // Fill window buffer + i = 0; + loop(fft_size, + r = i/fft_size; + window_buffer[i] = window(r); + i += 1; + ); +); + +pdc_delay = fft_size; +pdc_bot_ch = 0; +pdc_top_ch = 2; + +/*IFJSFX{ +freembuf(freemem); + +@slider +}IFJSFX*/ +//IFTEST function slider_code() ( +// convert low cut and high cut to bins every time a slider is changed +low_bin = min(slider5, slider6) / srate * fft_size; +high_bin = max(slider6, slider5) / srate * fft_size; +// convert to radians +rotation = slider7*$pi/180; +// convert percentage to raw scale factor +dry_mix = slider1/100; +wet_mix = slider2/100; +low_strength = slider3; +high_strength = slider4; +phase_width_low = slider8; +phase_width_high = slider9; +cosine = cos(rotation); +sine = sin(rotation); +cosine2 = cos(slider12*$pi/180); +sine2 = sin(slider12*$pi/180); +// fill strength_buffer and phase_width_buffer +band_index = 0; +loop(fft_size/2+1, + band_index >= low_bin && band_index <= high_bin ? + ( + // only set values for the appropriate frequency range + frac = (band_index - low_bin)/(high_bin - low_bin + 1); + // fraction of progress through range [low_bin, high_bin) + strength = low_strength* (1 - frac) + high_strength * frac; + strength_buffer[band_index] = exp(strength); + // precaculate strength (actual value should be positive, so it makes + // sense to take the power of ten, but only after the + // linear mapping over the spectrum is done. + // precalculate phase width + phase_width = phase_width_low * (1 - frac) + phase_width_high * frac; + phase_width_buffer[band_index] = exp(phase_width); + ); + + band_index += 1; + // next index +); +//IFTEST ); // slider_code() ( + +//IFTEST /* +@sample +//IFTEST */ +//IFTEST function sample_code() ( +input_buffer[buffer_index*2] = spl0; +input_buffer[buffer_index*2 + 1] = spl1; + +output_index = buffer_index - fft_size; +output_index < 0 ? output_index += buffer_length; +silence > 0 ? ( + spl0 = spl1 = 0; + // silence for fft init +) : ( + spl0 = output_buffer[output_index*2]; + spl1 = output_buffer[output_index*2 + 1]; +); +output_buffer[output_index*2] = 0; // clear the sample we just read +output_buffer[output_index*2 + 1] = 0; + +fft_counter += 1; +fft_counter >= fft_interval ? ( + fft_counter = 0; + + // Copy input to buffer + bi = buffer_index - fft_size + 1; + i = 0; + loop(fft_size, + i2 = bi + i; + i2 < 0 ? i2 += buffer_length; + + fft_buffer[2*i] = input_buffer[2*i2]*window_buffer[i]; + fft_buffer[2*i + 1] = input_buffer[2*i2 + 1]*window_buffer[i]; + + i += 1; + ); + + // Process buffer + fft(fft_buffer, fft_size); + fft_permute(fft_buffer, fft_size); + + process_stft_segment(fft_buffer, fft_size); + + fft_ipermute(fft_buffer, fft_size); + ifft(fft_buffer, fft_size); + + // Add to output + bi = buffer_index - fft_size + 1; + i = 0; + loop(fft_size, + i2 = bi + i; + (i2 < 0) ? i2 += buffer_length; + + output_buffer[2*i2] += fft_buffer[2*i]*fft_scaling_factor*window_buffer[i]; + output_buffer[2*i2 + 1] += fft_buffer[2*i + 1]*fft_scaling_factor*window_buffer[i]; + + i += 1; + ); + silence > 0 ? silence -= 1; +); + +buffer_index = (buffer_index + 1)%buffer_length; + +//IFTEST ); // function sample_code() + + +/*IFJSFX{ +@serialize +}IFJSFX*/ +//IFEEL /* +//IFJSFX /* +function file_var(file, val) ( + printf("FILE_VAR, FILE: %d, VAL: %g\n", file, val) +); +//IFJSFX */ +serial_version = 1.00; +file_var(0, serial_version); +//IFEEL */ +// nothing serialized yet, but keep track of the serial_version +// for the preset state of the original plugin, serial_version would now be euqal to 0. + +/*IFTEST{ // main test block + +// helpers +function sum_first_pdc_samples(s0val, s1val) ( + setup_stft_state(fft_size, 1); + spl0sum = 0; + spl1sum = 0; + loop(pdc_delay, + spl0=s0val; + spl1=s1val; + sample_code(); + spl0sum += abs(spl0); + spl1sum += abs(spl1); + ); +); + +printf("SETUP: fft range is [0, srate/2] dry=100 wet=100 attenuate_diff_vol=1\n"); +slider1 = slider2 = 100; +slider5 = 0; +slider6 = srate/2; +slider10 = 1; +slider_code(); +printf("SETUP: got low_bin=%g, high_bin=%g\n", low_bin, high_bin); + +printf("get_weight_ALL_ZERO\n"); +assert_equal_exact(0, get_weight(1, 1, 0, 0, 0, 0)); + +printf("get_weight_PERPENDICULAR\n"); +assert_equal_exact(0.25, get_weight(1, 1, 0, 1, 1, 0)); + +printf("get_weight_EQUIVALENT\n"); +assert_equal_exact(0.5, get_weight(1, 1, 1, 0, 1, 0)); + +printf("get_weight_INVERSE\n"); +assert_equal_exact(0, get_weight(1, 1, 1, 0, -1, 0)); + +printf("get_weight_ONE_ONLY\n"); +assert_equal_exact(0, get_weight(1, 1, 1, 0, 0, 0)); + +printf("get_weight_HALF\n"); +assert_equal_exact(0.5, get_weight(1, 1, 0.5, 0, 0.5, 0)); + +printf("get_weight_STRANGE\n"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -33.28894848631978, -662.24573314730117, -33.28894848631978, -662.24573314730117), "STRANGE01"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, 452.28715419177155, 53.67874245364300, 452.28715419177155, 53.67874245364300), "STRANGE02"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, 267.79446776740104, 17.14616043377475, 267.79446776740104, 17.14616043377475), "STRANGE03"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -32.55549374859478, 261.41171292127513, -32.55549374859478, 261.41171292127513), "STRANGE04"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -2.75689966923304, -217.67592724424455, -2.75689966923304, -217.67592724424455), "STRANGE05"); + +printf("SAMPLE_0_0_TEST\n"); +// spl0=0 spl1=0 for 100 output samples +sum_first_pdc_samples(0, 0); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +loop(100, + spl0=1; + spl1=1; + sample_code(); + spl0sum += abs(spl0); + spl1sum += abs(spl1); +); +assert_near_equal(200, 0.001, spl0sum, "SAMPLE_1_1_TEST: spl0sum was not as expected"); +assert_near_equal(200, 0.001, spl1sum, "SAMPLE_1_1_TEST: spl1sum was not as expected"); +printf("spl0sum=%g, spl1sum=%g\n", spl0sum, spl1sum); + + +printf("SETUP: fft range is [srate/4, srate/2] dry=100 wet=100\n"); +// check no values before first output sample +slider1 = slider2 = 100; +slider5 = srate/4; +slider6 = srate/2; +printf("PDC_SILENCE_ON_PARTIAL_RANGE_TEST\n"); +slider_code(); +sum_first_pdc_samples(1, 1); +assert_equal_exact(0, spl0sum, "PDC_SILENCE_ON_PARTIAL_RANGE_TEST failed on spl0"); +assert_equal_exact(0, spl1sum, "PDC_SILENCE_ON_PARTIAL_RANGE_TEST failed on spl1"); + +test_summary(); +}IFTEST*/ // main test block diff --git a/vocalredisoBlurry.eel b/vocalredisoBlurry.eel new file mode 100644 index 0000000..8be84e1 --- /dev/null +++ b/vocalredisoBlurry.eel @@ -0,0 +1,517 @@ +// Based on Vocalrediso.ny, a nyquist filter for audacity, currently released under the GNU GPL V3 +// see https://www.gnu.org/licenses/gpl-3.0.en.html, https://www.audacityteam.org/ +// credits to Neil Bickford for his denoiser, https://github.com/Nbickford/REAPERDenoiser, +// which was used as a starting point for this script +// which now uses the STFT template code by geraintluff: https://forum.cockos.com/showthread.php?t=225955 + +//IFTEST /* //this section not used by testing code +desc: vocal removal/isolation +//tags: processing vocals stereo +//author: Michael Pannekoek + +//slider1:fft_size_index=2<0,7,1{256,512,1024,2048,4096,8192,16384,32768}>FFT size +slider1:0<-100,100,0.1>dry mix +slider2:0<-100,100,0.1>C mix (Vocals) +slider3:0<-5,5,0.001>strength at Low Cut +slider4:0<-5,5,0.001>strength at High Cut +slider5:80<0,24000,10>Low Cut (Vocals) +slider6:24000<0,24000,10>High Cut (Vocals) +slider7:0<-90,90,0.1>Phase (Degrees) +slider8:0<-5,5,0.001>Phase width at Low Cut +slider9:0<-5,5,0.001>Phase width at High Cut +slider10:1<0,1,0.05>Attenuate if different volume +slider11:1<0,1,1{No,Yes}>undo input corrections +slider12:0<-180,180,0.05>Phase2 (Degrees) +slider13:0<0,10,0.1>Blurring factor +//IFTEST */ //this section not used by testing code + +/*IFJSFX{ +in_pin:left input +in_pin:right input +out_pin:left output +out_pin:right output +}IFJSFX*/ + + + +//IFTEST /* +@init +//IFTEST */ +//IFJSFX /* +slider1=0; +slider2=0; +slider3=0; +slider4=0; +slider5=80; +slider6=24000; +slider7=0; +slider8=0; +slider9=0; +slider10=1; +slider11=1; +slider12=0; +slider13=0; +//IFJSFX */ + +//IFTEST srate = 24000; + +_free_memory = 0; +function simple_alloc(amount) +local(_free_memory_old) +global(_free_memory) +( + _free_memory_old = _free_memory; + _free_memory += amount; + _free_memory_old; +); + +prev_fft_size = 0; + +function setup_stft_state(fft_size, first_time) ( + //////////////////////// Setup block + // This is called when playback starts, or when the FFT size is changed + memset(output_buffer, 0, buffer_length*2); + memset(input_buffer, 0, buffer_length*2); + // reset indexes and counters + buffer_index = 0; + output_index = 0; + fft_counter = 0; + // force silence for first overlaps where fft not yet run + silence = overlap_factor; + //////////////////////// +); + +function get_weight(strength, phase_width, l_real, l_imag, r_real, r_imag) local(norm_left, norm_right, w1, weight) ( + // cacluate energy for this bin + norm_left = sqrt(sqr(l_real) + sqr(l_imag)); + norm_right = sqrt(sqr(r_real) + sqr(r_imag)); + + // calculate phase difference between left & right, divide by phase_width + w1 = (l_real * r_real + l_imag * r_imag) / (norm_left * norm_right * phase_width * 2); + weight = ((max(0, min(1, w1+0.5))) ^ strength); + // calculate weight: truncate w1 to [0, 1] and apply strength, then take 1 - the result, and multiply + // by 1 - the square of the difference between the two norm values divided by the sum of the two, moderated by strength * slider10 + slider10 > 0 ? weight *= ( + 1 - sqr(norm_left - norm_right)/sqr(norm_left + norm_right) + ) ^ (slider10 / strength) * 0.5; + weight; +); + +setup_stft_state(fft_size, prev_fft_size == 0); + +function process_stft_segment(fft_buffer, fft_size) local(fft_bin, left_real, left_imag, right_real, right_imag) ( + fft_bin = 0; // FFT bin number + loop(fft_size/2+1, + fft_bin2 = fft_bin ? (fft_size - fft_bin) : 0; + + // Unfold complex spectrum into two real spectra + left_real = fft_buffer[2*fft_bin] + fft_buffer[2*fft_bin2]; + left_imag = fft_buffer[2*fft_bin + 1] - fft_buffer[2*fft_bin2 + 1]; + right_real = fft_buffer[2*fft_bin + 1] + fft_buffer[2*fft_bin2 + 1]; + right_imag = -fft_buffer[2*fft_bin] + fft_buffer[2*fft_bin2]; + + // input corrections + // first, change the phase of L based on phase2: + l_real_twisted = left_real*cosine2 + left_imag*sine2; + l_imag_twisted = left_imag*cosine2 - left_real*sine2; + + // now mix L&R together based on phase + r_real_rotated = right_real*cosine + l_real_twisted*sine; + r_imag_rotated = right_imag*cosine + l_imag_twisted*sine; + l_real_rotated = l_real_twisted*cosine - right_real*sine; + l_imag_rotated = l_imag_twisted*cosine - right_imag*sine; + + //////////////////////// Main STFT block + // The 'meat' of the algorithm, the code in this block will most resemble the code from vocalrediso.ny + + weight = 0; + // first, apply vocal reduction algorithm only in the right bands + fft_bin >= low_bin && fft_bin <= high_bin ? ( + weight = get_weight(strength_buffer[fft_bin], phase_width_buffer[fft_bin], l_real_rotated, l_imag_rotated, r_real_rotated, r_imag_rotated); + ); + + // memorize this bin's weight and coeffs + blurring_buffer[fft_bin] = weight; + + left_fft_buffer[fft_bin*2] = l_real_rotated; + left_fft_buffer[fft_bin*2+1] = l_imag_rotated; + right_fft_buffer[fft_bin*2] = r_real_rotated; + right_fft_buffer[fft_bin*2+1] = r_imag_rotated; + + fft_bin += 1; + ); + fft_bin = 0; // FFT bin number + loop(fft_size/2+1, + fft_bin2 = fft_bin ? (fft_size - fft_bin) : 0; + l_real_rotated = left_fft_buffer[fft_bin*2]; + l_imag_rotated = left_fft_buffer[fft_bin*2+1]; + r_real_rotated = right_fft_buffer[fft_bin*2]; + r_imag_rotated = right_fft_buffer[fft_bin*2+1]; + fft_bin >= low_bin && fft_bin <= high_bin ? ( + // find the c channel, the sum of the two complex numbers + c_real = l_real_rotated + r_real_rotated; + c_imag = l_imag_rotated + r_imag_rotated; + + // Find the blurred attenuation factor for this coefficient + // find out how much to blur by in this + blur_amount = (slider13*log(1+fft_bin))~0; + attenuation_factor = 0; + + index = max(0, fft_bin-blur_amount); + max+index = min(fft_size/2-1, fft_bin+blur_amount); + elems = max+index - index + 1; + loop(elems, + attenuation_factor += blurring_buffer[index]; + ); + + weight = attenuation_factor / elems; + + // apply weight to c channel + c_real *= weight; + c_imag *= weight; + ) : + ( + // let wet signal have 0 for fft coefficients when out of bounds + c_real = 0; + c_imag = 0; + ); + //////////////////////// END MAIN STFT block + + // apply wet dry mix + out_l_real = l_real_rotated * dry_mix + c_real * wet_mix; + out_l_imag = l_imag_rotated * dry_mix + c_imag * wet_mix; + out_r_real = r_real_rotated * dry_mix + c_real * wet_mix; + out_r_imag = r_imag_rotated * dry_mix + c_imag * wet_mix; + + // output corrections + slider11 > 0.5 ? ( + // if requested, undo input corrections + + // unmix by phase + l_real_out_twisted = out_l_real*cosine - out_r_real*sine; + l_imag_out_twisted = out_l_imag*cosine - out_r_imag*sine; + right_real = out_r_real*cosine + out_l_real*sine; + right_imag = out_r_imag*cosine + out_l_imag*sine; + + left_real = l_real_out_twisted * cosine2 - l_imag_out_twisted * sine2; + left_imag = l_imag_out_twisted * cosine2 + l_real_out_twisted * sine2; + ) : + ( + // else, just copy the values + left_real = out_l_real; + left_imag = out_l_imag; + right_real = out_r_real; + right_imag = out_r_imag; + ); + + // Re-fold back into complex spectrum + fft_buffer[2*fft_bin] = (left_real - right_imag)*0.5; + fft_buffer[2*fft_bin + 1] = (left_imag + right_real)*0.5; + fft_buffer[2*fft_bin2] = (left_real + right_imag)*0.5; + fft_buffer[2*fft_bin2 + 1] = (-left_imag + right_real)*0.5; + + fft_bin += 1; + ); +); + +MAX_FFT_SIZE = 32768; +fft_size = 8192; + +freemem = 0; +fft_buffer = simple_alloc(MAX_FFT_SIZE*2); +window_buffer = simple_alloc(MAX_FFT_SIZE); + +strength_buffer = simple_alloc(MAX_FFT_SIZE); +phase_width_buffer = simple_alloc(MAX_FFT_SIZE); +blurring_buffer = simple_alloc(MAX_FFT_SIZE); +left_fft_buffer = simple_alloc(MAX_FFT_SIZE*2+2); +right_fft_buffer = simple_alloc(MAX_FFT_SIZE*2+2); + +buffer_length = srate; +buffer_index = 0; +input_buffer = simple_alloc(buffer_length*2); +output_buffer = simple_alloc(buffer_length*2); + +function window(r) ( + // When squared, the Hann window adds up perfectly for overlap >= 4, so it's suitable for perfect reconstruction + //(0.5 - 0.5*cos(r*2*$pi))/sqrt(0.375); + // the MLT sine window also appears to add up correctly, with sigma = sqrt(2). + sin(r*$pi)*sqrt(2); +); + +overlap_factor = 4; +fft_interval = fft_size/overlap_factor; +fft_scaling_factor = 1/overlap_factor/fft_size; + +fft_size != prev_fft_size ? ( + setup_stft_state(fft_size, prev_fft_size == 0); + prev_fft_size = fft_size; + // Fill window buffer + i = 0; + loop(fft_size, + r = i/fft_size; + window_buffer[i] = window(r); + i += 1; + ); +); + +pdc_delay = fft_size; +pdc_bot_ch = 0; +pdc_top_ch = 2; + +/*IFJSFX{ +freembuf(freemem); + +@slider +}IFJSFX*/ +//IFTEST function slider_code() ( +// convert low cut and high cut to bins every time a slider is changed +low_bin = min(slider5, slider6) / srate * fft_size; +high_bin = max(slider6, slider5) / srate * fft_size; +// convert to radians +rotation = slider7*$pi/180; +// convert percentage to raw scale factor +dry_mix = slider1/100; +wet_mix = slider2/100; +low_strength = slider3; +high_strength = slider4; +phase_width_low = slider8; +phase_width_high = slider9; +cosine = cos(rotation); +sine = sin(rotation); +cosine2 = cos(slider12*$pi/180); +sine2 = sin(slider12*$pi/180); +// fill strength_buffer and phase_width_buffer +band_index = 0; +loop(fft_size/2+1, + band_index >= low_bin && band_index <= high_bin ? + ( + // only set values for the appropriate frequency range + frac = (band_index - low_bin)/(high_bin - low_bin + 1); + // fraction of progress through range [low_bin, high_bin) + strength = low_strength* (1 - frac) + high_strength * frac; + strength_buffer[band_index] = exp(strength); + // precaculate strength (actual value should be positive, so it makes + // sense to take the power of ten, but only after the + // linear mapping over the spectrum is done. + // precalculate phase width + phase_width = phase_width_low * (1 - frac) + phase_width_high * frac; + phase_width_buffer[band_index] = exp(phase_width); + ); + + band_index += 1; + // next index +); +//IFTEST ); // slider_code() ( + +//IFTEST /* +@sample +//IFTEST */ +//IFTEST function sample_code() ( +input_buffer[buffer_index*2] = spl0; +input_buffer[buffer_index*2 + 1] = spl1; + +output_index = buffer_index - fft_size; +output_index < 0 ? output_index += buffer_length; +silence > 0 ? ( + spl0 = spl1 = 0; + // silence for fft init +) : ( + spl0 = output_buffer[output_index*2]; + spl1 = output_buffer[output_index*2 + 1]; +); +output_buffer[output_index*2] = 0; // clear the sample we just read +output_buffer[output_index*2 + 1] = 0; + +fft_counter += 1; +fft_counter >= fft_interval ? ( + fft_counter = 0; + + // Copy input to buffer + bi = buffer_index - fft_size + 1; + i = 0; + loop(fft_size, + i2 = bi + i; + i2 < 0 ? i2 += buffer_length; + + fft_buffer[2*i] = input_buffer[2*i2]*window_buffer[i]; + fft_buffer[2*i + 1] = input_buffer[2*i2 + 1]*window_buffer[i]; + + i += 1; + ); + + // Process buffer + fft(fft_buffer, fft_size); + fft_permute(fft_buffer, fft_size); + + process_stft_segment(fft_buffer, fft_size); + + fft_ipermute(fft_buffer, fft_size); + ifft(fft_buffer, fft_size); + + // Add to output + bi = buffer_index - fft_size + 1; + i = 0; + loop(fft_size, + i2 = bi + i; + (i2 < 0) ? i2 += buffer_length; + + output_buffer[2*i2] += fft_buffer[2*i]*fft_scaling_factor*window_buffer[i]; + output_buffer[2*i2 + 1] += fft_buffer[2*i + 1]*fft_scaling_factor*window_buffer[i]; + + i += 1; + ); + silence > 0 ? silence -= 1; +); + +buffer_index = (buffer_index + 1)%buffer_length; + +//IFTEST ); // function sample_code() + + +/*IFJSFX{ +@serialize +}IFJSFX*/ + /* +//IFJSFX /* +function file_var(file, val) ( + printf("FILE_VAR, FILE: %d, VAL: %g\n", file, val) +); +//IFJSFX */ +serial_version = 1.00; +file_var(0, serial_version); + */ +// nothing serialized yet, but keep track of the serial_version +// for the preset state of the original plugin, serial_version would now be euqal to 0. + +/*IFTEST{ // main test block + +// helpers +function sum_first_pdc_samples(s0val, s1val) ( + setup_stft_state(fft_size, 1); + spl0sum = 0; + spl1sum = 0; + loop(pdc_delay, + spl0=s0val; + spl1=s1val; + sample_code(); + spl0sum += abs(spl0); + spl1sum += abs(spl1); + ); +); + +printf("SETUP: fft range is [0, srate/2] dry=100 wet=100 attenuate_diff_vol=1\n"); +slider1 = slider2 = 100; +slider5 = 0; +slider6 = srate/2; +slider10 = 1; +slider_code(); +printf("SETUP: got low_bin=%g, high_bin=%g\n", low_bin, high_bin); + +printf("get_weight_ALL_ZERO\n"); +assert_equal_exact(0, get_weight(1, 1, 0, 0, 0, 0)); + +printf("get_weight_PERPENDICULAR\n"); +assert_equal_exact(0.25, get_weight(1, 1, 0, 1, 1, 0)); + +printf("get_weight_EQUIVALENT\n"); +assert_equal_exact(0.5, get_weight(1, 1, 1, 0, 1, 0)); + +printf("get_weight_INVERSE\n"); +assert_equal_exact(0, get_weight(1, 1, 1, 0, -1, 0)); + +printf("get_weight_ONE_ONLY\n"); +assert_equal_exact(0, get_weight(1, 1, 1, 0, 0, 0)); + +printf("get_weight_HALF\n"); +assert_equal_exact(0.5, get_weight(1, 1, 0.5, 0, 0.5, 0)); + +printf("get_weight_STRANGE\n"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -33.28894848631978, -662.24573314730117, -33.28894848631978, -662.24573314730117), "STRANGE01"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, 452.28715419177155, 53.67874245364300, 452.28715419177155, 53.67874245364300), "STRANGE02"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, 267.79446776740104, 17.14616043377475, 267.79446776740104, 17.14616043377475), "STRANGE03"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -32.55549374859478, 261.41171292127513, -32.55549374859478, 261.41171292127513), "STRANGE04"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -2.75689966923304, -217.67592724424455, -2.75689966923304, -217.67592724424455), "STRANGE05"); + +printf("SAMPLE_0_0_TEST\n"); +// spl0=0 spl1=0 for 100 output samples +sum_first_pdc_samples(0, 0); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +loop(100, + spl0=1; + spl1=1; + sample_code(); + spl0sum += abs(spl0); + spl1sum += abs(spl1); +); +assert_near_equal(200, 0.001, spl0sum, "SAMPLE_1_1_TEST: spl0sum was not as expected"); +assert_near_equal(200, 0.001, spl1sum, "SAMPLE_1_1_TEST: spl1sum was not as expected"); +printf("spl0sum=%g, spl1sum=%g\n", spl0sum, spl1sum); + + +printf("SETUP: fft range is [srate/4, srate/2] dry=100 wet=100\n"); +// check no values before first output sample +slider1 = slider2 = 100; +slider5 = srate/4; +slider6 = srate/2; +printf("PDC_SILENCE_ON_PARTIAL_RANGE_TEST\n"); +slider_code(); +sum_first_pdc_samples(1, 1); +assert_equal_exact(0, spl0sum, "PDC_SILENCE_ON_PARTIAL_RANGE_TEST failed on spl0"); +assert_equal_exact(0, spl1sum, "PDC_SILENCE_ON_PARTIAL_RANGE_TEST failed on spl1"); + +test_summary(); +}IFTEST*/ // main test block diff --git a/vocalredisoBlurry.jsfx b/vocalredisoBlurry.jsfx index 59c90cf..33a31b4 100644 --- a/vocalredisoBlurry.jsfx +++ b/vocalredisoBlurry.jsfx @@ -4,10 +4,12 @@ // which was used as a starting point for this script // which now uses the STFT template code by geraintluff: https://forum.cockos.com/showthread.php?t=225955 +//IFTEST /* //this section not used by testing code desc: vocal removal/isolation //tags: processing vocals stereo //author: Michael Pannekoek +//slider1:fft_size_index=2<0,7,1{256,512,1024,2048,4096,8192,16384,32768}>FFT size slider1:0<-100,100,0.1>dry mix slider2:0<-100,100,0.1>C mix (Vocals) slider3:0<-5,5,0.001>strength at Low Cut @@ -15,12 +17,14 @@ slider4:0<-5,5,0.001>strength at High Cut slider5:80<0,24000,10>Low Cut (Vocals) slider6:24000<0,24000,10>High Cut (Vocals) slider7:0<-90,90,0.1>Phase (Degrees) -slider8:90<1,180,0.1>Phase width at Low Cut (Degrees) -slider9:90<1,180,0.1>Phase width at High Cut (Degrees) +slider8:0<-5,5,0.001>Phase width at Low Cut +slider9:0<-5,5,0.001>Phase width at High Cut slider10:1<0,1,0.05>Attenuate if different volume slider11:1<0,1,1{No,Yes}>undo input corrections slider12:0<-180,180,0.05>Phase2 (Degrees) slider13:0<0,10,0.1>Blurring factor +//IFTEST */ //this section not used by testing code + in_pin:left input in_pin:right input @@ -28,7 +32,29 @@ out_pin:left output out_pin:right output + + +//IFTEST /* @init +//IFTEST */ + /* +slider1=0; +slider2=0; +slider3=0; +slider4=0; +slider5=80; +slider6=24000; +slider7=0; +slider8=0; +slider9=0; +slider10=1; +slider11=1; +slider12=0; +slider13=0; + */ + +//IFTEST srate = 24000; + _free_memory = 0; function simple_alloc(amount) local(_free_memory_old) @@ -44,11 +70,33 @@ prev_fft_size = 0; function setup_stft_state(fft_size, first_time) ( //////////////////////// Setup block // This is called when playback starts, or when the FFT size is changed - memset(fft_buffer, 0.0, MAX_FFT_SIZE*2); - memset(output_buffer, 0.0, buffer_length*2); + memset(output_buffer, 0, buffer_length*2); + memset(input_buffer, 0, buffer_length*2); + // reset indexes and counters + buffer_index = 0; + output_index = 0; + fft_counter = 0; + // force silence for first overlaps where fft not yet run + silence = overlap_factor; //////////////////////// ); +function get_weight(strength, phase_width, l_real, l_imag, r_real, r_imag) local(norm_left, norm_right, w1, weight) ( + // cacluate energy for this bin + norm_left = sqrt(sqr(l_real) + sqr(l_imag)); + norm_right = sqrt(sqr(r_real) + sqr(r_imag)); + + // calculate phase difference between left & right, divide by phase_width + w1 = (l_real * r_real + l_imag * r_imag) / (norm_left * norm_right * phase_width * 2); + weight = ((max(0, min(1, w1+0.5))) ^ strength); + // calculate weight: truncate w1 to [0, 1] and apply strength, then take 1 - the result, and multiply + // by 1 - the square of the difference between the two norm values divided by the sum of the two, moderated by strength * slider10 + slider10 > 0 ? weight *= ( + 1 - sqr(norm_left - norm_right)/sqr(norm_left + norm_right) + ) ^ (slider10 / strength) * 0.5; + weight; +); + setup_stft_state(fft_size, prev_fft_size == 0); function process_stft_segment(fft_buffer, fft_size) local(fft_bin, left_real, left_imag, right_real, right_imag) ( @@ -76,44 +124,29 @@ function process_stft_segment(fft_buffer, fft_size) local(fft_bin, left_real, le //////////////////////// Main STFT block // The 'meat' of the algorithm, the code in this block will most resemble the code from vocalrediso.ny - weight = 1; - - // first, apply vocal reduction algorithm only in the right bands + weight = 0; + // first, apply vocal reduction algorithm only in the right bands fft_bin >= low_bin && fft_bin <= high_bin ? ( - // get constants for this bin - strength = strength_buffer[fft_bin]; - phase_width = phase_width_buffer[fft_bin]; - - // cacluate energy for this bin - norm_left = sqrt(sqr(l_real_rotated) + sqr(l_imag_rotated)); - norm_right = sqrt(sqr(r_real_rotated) + sqr(r_imag_rotated)); - - // calculate phase difference between left & right, divide by phase_width - w1 = acos((l_real_rotated * r_real_rotated + l_imag_rotated * r_imag_rotated) / (norm_left * norm_right)) / phase_width; - // calculate weight: truncate w1 to [0, 1] and apply strength, then take 1 - the result, and multiply - // by 1 - the square of the difference between the two norm values divided by the sum of the two, moderated by strength * slider10 - weight = (1 - min(1, max(0, w1)) ^ strength) * ( - 1 - (sqr(norm_left - norm_right)/sqr(norm_left + norm_right)) - ) ^ (slider10 / strength) / 2; - - // memorize this bin's weight + weight = get_weight(strength_buffer[fft_bin], phase_width_buffer[fft_bin], l_real_rotated, l_imag_rotated, r_real_rotated, r_imag_rotated); ); - blurringBuffer[fft_bin] = weight; + // memorize this bin's weight and coeffs + blurring_buffer[fft_bin] = weight; - fft_buffer[2*fft_bin] = l_real_rotated; - fft_buffer[2*fft_bin + 1] = l_imag_rotated; - fft_buffer[2*fft_bin2] = r_real_rotated; - fft_buffer[2*fft_bin2 + 1] = r_imag_rotated; + left_fft_buffer[fft_bin*2] = l_real_rotated; + left_fft_buffer[fft_bin*2+1] = l_imag_rotated; + right_fft_buffer[fft_bin*2] = r_real_rotated; + right_fft_buffer[fft_bin*2+1] = r_imag_rotated; fft_bin += 1; ); fft_bin = 0; // FFT bin number loop(fft_size/2+1, - l_real_rotated = fft_buffer[2*fft_bin] - l_imag_rotated = fft_buffer[2*fft_bin + 1] - r_real_rotated = fft_buffer[2*fft_bin2] - r_imag_rotated = fft_buffer[2*fft_bin2 + 1] + fft_bin2 = fft_bin ? (fft_size - fft_bin) : 0; + l_real_rotated = left_fft_buffer[fft_bin*2]; + l_imag_rotated = left_fft_buffer[fft_bin*2+1]; + r_real_rotated = right_fft_buffer[fft_bin*2]; + r_imag_rotated = right_fft_buffer[fft_bin*2+1]; fft_bin >= low_bin && fft_bin <= high_bin ? ( // find the c channel, the sum of the two complex numbers c_real = l_real_rotated + r_real_rotated; @@ -121,24 +154,26 @@ function process_stft_segment(fft_buffer, fft_size) local(fft_bin, left_real, le // Find the blurred attenuation factor for this coefficient // find out how much to blur by in this - blurAmount = (slider13*log(1+bandIndex))~0; - attenuationFactor = 0; - divisor = 1; - - index = max(0, bandIndex-blurAmount); - maxIndex = min(SIZE/2-1, bandIndex+blurAmount); - elems = maxIndex - index + 1; + blur_amount = (slider13*log(1+fft_bin))~0; + attenuation_factor = 0; + + index = max(0, fft_bin-blur_amount); + max+index = min(fft_size/2-1, fft_bin+blur_amount); + elems = max+index - index + 1; loop(elems, - attenuationFactor = max(attenuationFactor, blurringBuffer[index]); + attenuation_factor += blurring_buffer[index]; ); - weight = attenuationFactor; - - centerReal = (LrealDry + RrealDry) * weight; - centerImag = (LimagDry + RimagDry) * weight; - ) : ( + + weight = attenuation_factor / elems; + + // apply weight to c channel + c_real *= weight; + c_imag *= weight; + ) : + ( // let wet signal have 0 for fft coefficients when out of bounds - centerReal = 0; - centerImag = 0; + c_real = 0; + c_imag = 0; ); //////////////////////// END MAIN STFT block @@ -188,13 +223,16 @@ window_buffer = simple_alloc(MAX_FFT_SIZE); strength_buffer = simple_alloc(MAX_FFT_SIZE); phase_width_buffer = simple_alloc(MAX_FFT_SIZE); +blurring_buffer = simple_alloc(MAX_FFT_SIZE); +left_fft_buffer = simple_alloc(MAX_FFT_SIZE*2+2); +right_fft_buffer = simple_alloc(MAX_FFT_SIZE*2+2); buffer_length = srate; buffer_index = 0; input_buffer = simple_alloc(buffer_length*2); output_buffer = simple_alloc(buffer_length*2); -function window(r) local(s, s2, gaussian_width, x) ( +function window(r) ( // When squared, the Hann window adds up perfectly for overlap >= 4, so it's suitable for perfect reconstruction //(0.5 - 0.5*cos(r*2*$pi))/sqrt(0.375); // the MLT sine window also appears to add up correctly, with sigma = sqrt(2). @@ -221,10 +259,12 @@ pdc_delay = fft_size; pdc_bot_ch = 0; pdc_top_ch = 2; -freembuf(freemem); +freembuf(freemem); @slider + +//IFTEST function slider_code() ( // convert low cut and high cut to bins every time a slider is changed low_bin = min(slider5, slider6) / srate * fft_size; high_bin = max(slider6, slider5) / srate * fft_size; @@ -235,40 +275,54 @@ dry_mix = slider1/100; wet_mix = slider2/100; low_strength = slider3; high_strength = slider4; -phase_width_low = slider8*$pi/180; -phase_width_high = slider9*$pi/180; +phase_width_low = slider8; +phase_width_high = slider9; cosine = cos(rotation); sine = sin(rotation); cosine2 = cos(slider12*$pi/180); sine2 = sin(slider12*$pi/180); // fill strength_buffer and phase_width_buffer band_index = 0; -loop(fft_size, +loop(fft_size/2+1, band_index >= low_bin && band_index <= high_bin ? ( // only set values for the appropriate frequency range - frac = (band_index - low_bin)/(high_bin - low_bin - 1); - frac = max(0, min(1, frac)); + frac = (band_index - low_bin)/(high_bin - low_bin + 1); // fraction of progress through range [low_bin, high_bin) strength = low_strength* (1 - frac) + high_strength * frac; - strength_buffer[band_index] = 10^strength; + strength_buffer[band_index] = exp(strength); // precaculate strength (actual value should be positive, so it makes // sense to take the power of ten, but only after the // linear mapping over the spectrum is done. // precalculate phase width - phase_width_buffer[band_index] = phase_width_low * (1 - frac) + phase_width_high * frac; + phase_width = phase_width_low * (1 - frac) + phase_width_high * frac; + phase_width_buffer[band_index] = exp(phase_width); ); band_index += 1; // next index ); +//IFTEST ); // slider_code() ( - +//IFTEST /* @sample - +//IFTEST */ +//IFTEST function sample_code() ( input_buffer[buffer_index*2] = spl0; input_buffer[buffer_index*2 + 1] = spl1; +output_index = buffer_index - fft_size; +output_index < 0 ? output_index += buffer_length; +silence > 0 ? ( + spl0 = spl1 = 0; + // silence for fft init +) : ( + spl0 = output_buffer[output_index*2]; + spl1 = output_buffer[output_index*2 + 1]; +); +output_buffer[output_index*2] = 0; // clear the sample we just read +output_buffer[output_index*2 + 1] = 0; + fft_counter += 1; fft_counter >= fft_interval ? ( fft_counter = 0; @@ -307,19 +361,157 @@ fft_counter >= fft_interval ? ( i += 1; ); + silence > 0 ? silence -= 1; ); -output_index = buffer_index - fft_size; -output_index < 0 ? output_index += buffer_length; -spl0 = output_buffer[output_index*2]; -spl1 = output_buffer[output_index*2 + 1]; -output_buffer[output_index*2] = 0; // clear the sample we just read -output_buffer[output_index*2 + 1] = 0; - buffer_index = (buffer_index + 1)%buffer_length; +//IFTEST ); // function sample_code() + + + @serialize + +//IFEEL /* + /* +function file_var(file, val) ( + printf("FILE_VAR, FILE: %d, VAL: %g\n", file, val) +); + */ serial_version = 1.00; file_var(0, serial_version); +//IFEEL */ // nothing serialized yet, but keep track of the serial_version // for the preset state of the original plugin, serial_version would now be euqal to 0. + +/*IFTEST{ // main test block + +// helpers +function sum_first_pdc_samples(s0val, s1val) ( + setup_stft_state(fft_size, 1); + spl0sum = 0; + spl1sum = 0; + loop(pdc_delay, + spl0=s0val; + spl1=s1val; + sample_code(); + spl0sum += abs(spl0); + spl1sum += abs(spl1); + ); +); + +printf("SETUP: fft range is [0, srate/2] dry=100 wet=100 attenuate_diff_vol=1\n"); +slider1 = slider2 = 100; +slider5 = 0; +slider6 = srate/2; +slider10 = 1; +slider_code(); +printf("SETUP: got low_bin=%g, high_bin=%g\n", low_bin, high_bin); + +printf("get_weight_ALL_ZERO\n"); +assert_equal_exact(0, get_weight(1, 1, 0, 0, 0, 0)); + +printf("get_weight_PERPENDICULAR\n"); +assert_equal_exact(0.25, get_weight(1, 1, 0, 1, 1, 0)); + +printf("get_weight_EQUIVALENT\n"); +assert_equal_exact(0.5, get_weight(1, 1, 1, 0, 1, 0)); + +printf("get_weight_INVERSE\n"); +assert_equal_exact(0, get_weight(1, 1, 1, 0, -1, 0)); + +printf("get_weight_ONE_ONLY\n"); +assert_equal_exact(0, get_weight(1, 1, 1, 0, 0, 0)); + +printf("get_weight_HALF\n"); +assert_equal_exact(0.5, get_weight(1, 1, 0.5, 0, 0.5, 0)); + +printf("get_weight_STRANGE\n"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -33.28894848631978, -662.24573314730117, -33.28894848631978, -662.24573314730117), "STRANGE01"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, 452.28715419177155, 53.67874245364300, 452.28715419177155, 53.67874245364300), "STRANGE02"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, 267.79446776740104, 17.14616043377475, 267.79446776740104, 17.14616043377475), "STRANGE03"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -32.55549374859478, 261.41171292127513, -32.55549374859478, 261.41171292127513), "STRANGE04"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -2.75689966923304, -217.67592724424455, -2.75689966923304, -217.67592724424455), "STRANGE05"); + +printf("SAMPLE_0_0_TEST\n"); +// spl0=0 spl1=0 for 100 output samples +sum_first_pdc_samples(0, 0); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +loop(100, + spl0=1; + spl1=1; + sample_code(); + spl0sum += abs(spl0); + spl1sum += abs(spl1); +); +assert_near_equal(200, 0.001, spl0sum, "SAMPLE_1_1_TEST: spl0sum was not as expected"); +assert_near_equal(200, 0.001, spl1sum, "SAMPLE_1_1_TEST: spl1sum was not as expected"); +printf("spl0sum=%g, spl1sum=%g\n", spl0sum, spl1sum); + + +printf("SETUP: fft range is [srate/4, srate/2] dry=100 wet=100\n"); +// check no values before first output sample +slider1 = slider2 = 100; +slider5 = srate/4; +slider6 = srate/2; +printf("PDC_SILENCE_ON_PARTIAL_RANGE_TEST\n"); +slider_code(); +sum_first_pdc_samples(1, 1); +assert_equal_exact(0, spl0sum, "PDC_SILENCE_ON_PARTIAL_RANGE_TEST failed on spl0"); +assert_equal_exact(0, spl1sum, "PDC_SILENCE_ON_PARTIAL_RANGE_TEST failed on spl1"); + +test_summary(); +}IFTEST*/ // main test block diff --git a/vocalredisoBlurry.jsfx-template b/vocalredisoBlurry.jsfx-template new file mode 100644 index 0000000..35aded3 --- /dev/null +++ b/vocalredisoBlurry.jsfx-template @@ -0,0 +1,517 @@ +// Based on Vocalrediso.ny, a nyquist filter for audacity, currently released under the GNU GPL V3 +// see https://www.gnu.org/licenses/gpl-3.0.en.html, https://www.audacityteam.org/ +// credits to Neil Bickford for his denoiser, https://github.com/Nbickford/REAPERDenoiser, +// which was used as a starting point for this script +// which now uses the STFT template code by geraintluff: https://forum.cockos.com/showthread.php?t=225955 + +//IFTEST /* //this section not used by testing code +desc: vocal removal/isolation +//tags: processing vocals stereo +//author: Michael Pannekoek + +//slider1:fft_size_index=2<0,7,1{256,512,1024,2048,4096,8192,16384,32768}>FFT size +slider1:0<-100,100,0.1>dry mix +slider2:0<-100,100,0.1>C mix (Vocals) +slider3:0<-5,5,0.001>strength at Low Cut +slider4:0<-5,5,0.001>strength at High Cut +slider5:80<0,24000,10>Low Cut (Vocals) +slider6:24000<0,24000,10>High Cut (Vocals) +slider7:0<-90,90,0.1>Phase (Degrees) +slider8:0<-5,5,0.001>Phase width at Low Cut +slider9:0<-5,5,0.001>Phase width at High Cut +slider10:1<0,1,0.05>Attenuate if different volume +slider11:1<0,1,1{No,Yes}>undo input corrections +slider12:0<-180,180,0.05>Phase2 (Degrees) +slider13:0<0,10,0.1>Blurring factor +//IFTEST */ //this section not used by testing code + +/*IFJSFX{ +in_pin:left input +in_pin:right input +out_pin:left output +out_pin:right output +}IFJSFX*/ + + + +//IFTEST /* +@init +//IFTEST */ +//IFJSFX /* +slider1=0; +slider2=0; +slider3=0; +slider4=0; +slider5=80; +slider6=24000; +slider7=0; +slider8=0; +slider9=0; +slider10=1; +slider11=1; +slider12=0; +slider13=0; +//IFJSFX */ + +//IFTEST srate = 24000; + +_free_memory = 0; +function simple_alloc(amount) +local(_free_memory_old) +global(_free_memory) +( + _free_memory_old = _free_memory; + _free_memory += amount; + _free_memory_old; +); + +prev_fft_size = 0; + +function setup_stft_state(fft_size, first_time) ( + //////////////////////// Setup block + // This is called when playback starts, or when the FFT size is changed + memset(output_buffer, 0, buffer_length*2); + memset(input_buffer, 0, buffer_length*2); + // reset indexes and counters + buffer_index = 0; + output_index = 0; + fft_counter = 0; + // force silence for first overlaps where fft not yet run + silence = overlap_factor; + //////////////////////// +); + +function get_weight(strength, phase_width, l_real, l_imag, r_real, r_imag) local(norm_left, norm_right, w1, weight) ( + // cacluate energy for this bin + norm_left = sqrt(sqr(l_real) + sqr(l_imag)); + norm_right = sqrt(sqr(r_real) + sqr(r_imag)); + + // calculate phase difference between left & right, divide by phase_width + w1 = (l_real * r_real + l_imag * r_imag) / (norm_left * norm_right * phase_width * 2); + weight = ((max(0, min(1, w1+0.5))) ^ strength); + // calculate weight: truncate w1 to [0, 1] and apply strength, then take 1 - the result, and multiply + // by 1 - the square of the difference between the two norm values divided by the sum of the two, moderated by strength * slider10 + slider10 > 0 ? weight *= ( + 1 - sqr(norm_left - norm_right)/sqr(norm_left + norm_right) + ) ^ (slider10 / strength) * 0.5; + weight; +); + +setup_stft_state(fft_size, prev_fft_size == 0); + +function process_stft_segment(fft_buffer, fft_size) local(fft_bin, left_real, left_imag, right_real, right_imag) ( + fft_bin = 0; // FFT bin number + loop(fft_size/2+1, + fft_bin2 = fft_bin ? (fft_size - fft_bin) : 0; + + // Unfold complex spectrum into two real spectra + left_real = fft_buffer[2*fft_bin] + fft_buffer[2*fft_bin2]; + left_imag = fft_buffer[2*fft_bin + 1] - fft_buffer[2*fft_bin2 + 1]; + right_real = fft_buffer[2*fft_bin + 1] + fft_buffer[2*fft_bin2 + 1]; + right_imag = -fft_buffer[2*fft_bin] + fft_buffer[2*fft_bin2]; + + // input corrections + // first, change the phase of L based on phase2: + l_real_twisted = left_real*cosine2 + left_imag*sine2; + l_imag_twisted = left_imag*cosine2 - left_real*sine2; + + // now mix L&R together based on phase + r_real_rotated = right_real*cosine + l_real_twisted*sine; + r_imag_rotated = right_imag*cosine + l_imag_twisted*sine; + l_real_rotated = l_real_twisted*cosine - right_real*sine; + l_imag_rotated = l_imag_twisted*cosine - right_imag*sine; + + //////////////////////// Main STFT block + // The 'meat' of the algorithm, the code in this block will most resemble the code from vocalrediso.ny + + weight = 0; + // first, apply vocal reduction algorithm only in the right bands + fft_bin >= low_bin && fft_bin <= high_bin ? ( + weight = get_weight(strength_buffer[fft_bin], phase_width_buffer[fft_bin], l_real_rotated, l_imag_rotated, r_real_rotated, r_imag_rotated); + ); + + // memorize this bin's weight and coeffs + blurring_buffer[fft_bin] = weight; + + left_fft_buffer[fft_bin*2] = l_real_rotated; + left_fft_buffer[fft_bin*2+1] = l_imag_rotated; + right_fft_buffer[fft_bin*2] = r_real_rotated; + right_fft_buffer[fft_bin*2+1] = r_imag_rotated; + + fft_bin += 1; + ); + fft_bin = 0; // FFT bin number + loop(fft_size/2+1, + fft_bin2 = fft_bin ? (fft_size - fft_bin) : 0; + l_real_rotated = left_fft_buffer[fft_bin*2]; + l_imag_rotated = left_fft_buffer[fft_bin*2+1]; + r_real_rotated = right_fft_buffer[fft_bin*2]; + r_imag_rotated = right_fft_buffer[fft_bin*2+1]; + fft_bin >= low_bin && fft_bin <= high_bin ? ( + // find the c channel, the sum of the two complex numbers + c_real = l_real_rotated + r_real_rotated; + c_imag = l_imag_rotated + r_imag_rotated; + + // Find the blurred attenuation factor for this coefficient + // find out how much to blur by in this + blur_amount = (slider13*log(1+fft_bin))~0; + attenuation_factor = 0; + + index = max(0, fft_bin-blur_amount); + max+index = min(fft_size/2-1, fft_bin+blur_amount); + elems = max+index - index + 1; + loop(elems, + attenuation_factor += blurring_buffer[index]; + ); + + weight = attenuation_factor / elems; + + // apply weight to c channel + c_real *= weight; + c_imag *= weight; + ) : + ( + // let wet signal have 0 for fft coefficients when out of bounds + c_real = 0; + c_imag = 0; + ); + //////////////////////// END MAIN STFT block + + // apply wet dry mix + out_l_real = l_real_rotated * dry_mix + c_real * wet_mix; + out_l_imag = l_imag_rotated * dry_mix + c_imag * wet_mix; + out_r_real = r_real_rotated * dry_mix + c_real * wet_mix; + out_r_imag = r_imag_rotated * dry_mix + c_imag * wet_mix; + + // output corrections + slider11 > 0.5 ? ( + // if requested, undo input corrections + + // unmix by phase + l_real_out_twisted = out_l_real*cosine - out_r_real*sine; + l_imag_out_twisted = out_l_imag*cosine - out_r_imag*sine; + right_real = out_r_real*cosine + out_l_real*sine; + right_imag = out_r_imag*cosine + out_l_imag*sine; + + left_real = l_real_out_twisted * cosine2 - l_imag_out_twisted * sine2; + left_imag = l_imag_out_twisted * cosine2 + l_real_out_twisted * sine2; + ) : + ( + // else, just copy the values + left_real = out_l_real; + left_imag = out_l_imag; + right_real = out_r_real; + right_imag = out_r_imag; + ); + + // Re-fold back into complex spectrum + fft_buffer[2*fft_bin] = (left_real - right_imag)*0.5; + fft_buffer[2*fft_bin + 1] = (left_imag + right_real)*0.5; + fft_buffer[2*fft_bin2] = (left_real + right_imag)*0.5; + fft_buffer[2*fft_bin2 + 1] = (-left_imag + right_real)*0.5; + + fft_bin += 1; + ); +); + +MAX_FFT_SIZE = 32768; +fft_size = 8192; + +freemem = 0; +fft_buffer = simple_alloc(MAX_FFT_SIZE*2); +window_buffer = simple_alloc(MAX_FFT_SIZE); + +strength_buffer = simple_alloc(MAX_FFT_SIZE); +phase_width_buffer = simple_alloc(MAX_FFT_SIZE); +blurring_buffer = simple_alloc(MAX_FFT_SIZE); +left_fft_buffer = simple_alloc(MAX_FFT_SIZE*2+2); +right_fft_buffer = simple_alloc(MAX_FFT_SIZE*2+2); + +buffer_length = srate; +buffer_index = 0; +input_buffer = simple_alloc(buffer_length*2); +output_buffer = simple_alloc(buffer_length*2); + +function window(r) ( + // When squared, the Hann window adds up perfectly for overlap >= 4, so it's suitable for perfect reconstruction + //(0.5 - 0.5*cos(r*2*$pi))/sqrt(0.375); + // the MLT sine window also appears to add up correctly, with sigma = sqrt(2). + sin(r*$pi)*sqrt(2); +); + +overlap_factor = 4; +fft_interval = fft_size/overlap_factor; +fft_scaling_factor = 1/overlap_factor/fft_size; + +fft_size != prev_fft_size ? ( + setup_stft_state(fft_size, prev_fft_size == 0); + prev_fft_size = fft_size; + // Fill window buffer + i = 0; + loop(fft_size, + r = i/fft_size; + window_buffer[i] = window(r); + i += 1; + ); +); + +pdc_delay = fft_size; +pdc_bot_ch = 0; +pdc_top_ch = 2; + +/*IFJSFX{ +freembuf(freemem); + +@slider +}IFJSFX*/ +//IFTEST function slider_code() ( +// convert low cut and high cut to bins every time a slider is changed +low_bin = min(slider5, slider6) / srate * fft_size; +high_bin = max(slider6, slider5) / srate * fft_size; +// convert to radians +rotation = slider7*$pi/180; +// convert percentage to raw scale factor +dry_mix = slider1/100; +wet_mix = slider2/100; +low_strength = slider3; +high_strength = slider4; +phase_width_low = slider8; +phase_width_high = slider9; +cosine = cos(rotation); +sine = sin(rotation); +cosine2 = cos(slider12*$pi/180); +sine2 = sin(slider12*$pi/180); +// fill strength_buffer and phase_width_buffer +band_index = 0; +loop(fft_size/2+1, + band_index >= low_bin && band_index <= high_bin ? + ( + // only set values for the appropriate frequency range + frac = (band_index - low_bin)/(high_bin - low_bin + 1); + // fraction of progress through range [low_bin, high_bin) + strength = low_strength* (1 - frac) + high_strength * frac; + strength_buffer[band_index] = exp(strength); + // precaculate strength (actual value should be positive, so it makes + // sense to take the power of ten, but only after the + // linear mapping over the spectrum is done. + // precalculate phase width + phase_width = phase_width_low * (1 - frac) + phase_width_high * frac; + phase_width_buffer[band_index] = exp(phase_width); + ); + + band_index += 1; + // next index +); +//IFTEST ); // slider_code() ( + +//IFTEST /* +@sample +//IFTEST */ +//IFTEST function sample_code() ( +input_buffer[buffer_index*2] = spl0; +input_buffer[buffer_index*2 + 1] = spl1; + +output_index = buffer_index - fft_size; +output_index < 0 ? output_index += buffer_length; +silence > 0 ? ( + spl0 = spl1 = 0; + // silence for fft init +) : ( + spl0 = output_buffer[output_index*2]; + spl1 = output_buffer[output_index*2 + 1]; +); +output_buffer[output_index*2] = 0; // clear the sample we just read +output_buffer[output_index*2 + 1] = 0; + +fft_counter += 1; +fft_counter >= fft_interval ? ( + fft_counter = 0; + + // Copy input to buffer + bi = buffer_index - fft_size + 1; + i = 0; + loop(fft_size, + i2 = bi + i; + i2 < 0 ? i2 += buffer_length; + + fft_buffer[2*i] = input_buffer[2*i2]*window_buffer[i]; + fft_buffer[2*i + 1] = input_buffer[2*i2 + 1]*window_buffer[i]; + + i += 1; + ); + + // Process buffer + fft(fft_buffer, fft_size); + fft_permute(fft_buffer, fft_size); + + process_stft_segment(fft_buffer, fft_size); + + fft_ipermute(fft_buffer, fft_size); + ifft(fft_buffer, fft_size); + + // Add to output + bi = buffer_index - fft_size + 1; + i = 0; + loop(fft_size, + i2 = bi + i; + (i2 < 0) ? i2 += buffer_length; + + output_buffer[2*i2] += fft_buffer[2*i]*fft_scaling_factor*window_buffer[i]; + output_buffer[2*i2 + 1] += fft_buffer[2*i + 1]*fft_scaling_factor*window_buffer[i]; + + i += 1; + ); + silence > 0 ? silence -= 1; +); + +buffer_index = (buffer_index + 1)%buffer_length; + +//IFTEST ); // function sample_code() + + +/*IFJSFX{ +@serialize +}IFJSFX*/ +//IFEEL /* +//IFJSFX /* +function file_var(file, val) ( + printf("FILE_VAR, FILE: %d, VAL: %g\n", file, val) +); +//IFJSFX */ +serial_version = 1.00; +file_var(0, serial_version); +//IFEEL */ +// nothing serialized yet, but keep track of the serial_version +// for the preset state of the original plugin, serial_version would now be euqal to 0. + +/*IFTEST{ // main test block + +// helpers +function sum_first_pdc_samples(s0val, s1val) ( + setup_stft_state(fft_size, 1); + spl0sum = 0; + spl1sum = 0; + loop(pdc_delay, + spl0=s0val; + spl1=s1val; + sample_code(); + spl0sum += abs(spl0); + spl1sum += abs(spl1); + ); +); + +printf("SETUP: fft range is [0, srate/2] dry=100 wet=100 attenuate_diff_vol=1\n"); +slider1 = slider2 = 100; +slider5 = 0; +slider6 = srate/2; +slider10 = 1; +slider_code(); +printf("SETUP: got low_bin=%g, high_bin=%g\n", low_bin, high_bin); + +printf("get_weight_ALL_ZERO\n"); +assert_equal_exact(0, get_weight(1, 1, 0, 0, 0, 0)); + +printf("get_weight_PERPENDICULAR\n"); +assert_equal_exact(0.25, get_weight(1, 1, 0, 1, 1, 0)); + +printf("get_weight_EQUIVALENT\n"); +assert_equal_exact(0.5, get_weight(1, 1, 1, 0, 1, 0)); + +printf("get_weight_INVERSE\n"); +assert_equal_exact(0, get_weight(1, 1, 1, 0, -1, 0)); + +printf("get_weight_ONE_ONLY\n"); +assert_equal_exact(0, get_weight(1, 1, 1, 0, 0, 0)); + +printf("get_weight_HALF\n"); +assert_equal_exact(0.5, get_weight(1, 1, 0.5, 0, 0.5, 0)); + +printf("get_weight_STRANGE\n"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -33.28894848631978, -662.24573314730117, -33.28894848631978, -662.24573314730117), "STRANGE01"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, 452.28715419177155, 53.67874245364300, 452.28715419177155, 53.67874245364300), "STRANGE02"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, 267.79446776740104, 17.14616043377475, 267.79446776740104, 17.14616043377475), "STRANGE03"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -32.55549374859478, 261.41171292127513, -32.55549374859478, 261.41171292127513), "STRANGE04"); +assert_near_equal(0.5, 0.0001, get_weight(1, 1, -2.75689966923304, -217.67592724424455, -2.75689966923304, -217.67592724424455), "STRANGE05"); + +printf("SAMPLE_0_0_TEST\n"); +// spl0=0 spl1=0 for 100 output samples +sum_first_pdc_samples(0, 0); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +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); +assert_equal_exact(0, spl0sum, "spl0sum for init"); +assert_equal_exact(0, spl1sum, "spl1sum for init"); + +loop(100, + spl0=1; + spl1=1; + sample_code(); + spl0sum += abs(spl0); + spl1sum += abs(spl1); +); +assert_near_equal(200, 0.001, spl0sum, "SAMPLE_1_1_TEST: spl0sum was not as expected"); +assert_near_equal(200, 0.001, spl1sum, "SAMPLE_1_1_TEST: spl1sum was not as expected"); +printf("spl0sum=%g, spl1sum=%g\n", spl0sum, spl1sum); + + +printf("SETUP: fft range is [srate/4, srate/2] dry=100 wet=100\n"); +// check no values before first output sample +slider1 = slider2 = 100; +slider5 = srate/4; +slider6 = srate/2; +printf("PDC_SILENCE_ON_PARTIAL_RANGE_TEST\n"); +slider_code(); +sum_first_pdc_samples(1, 1); +assert_equal_exact(0, spl0sum, "PDC_SILENCE_ON_PARTIAL_RANGE_TEST failed on spl0"); +assert_equal_exact(0, spl1sum, "PDC_SILENCE_ON_PARTIAL_RANGE_TEST failed on spl1"); + +test_summary(); +}IFTEST*/ // main test block