desc: FFT Multi Tool. Read docs. in ALPHA
//tags: requested by James HE , modded by geraintluff, modded by junh1024,
//via https://forum.cockos.com/showpost.php?p=1962348&postcount=18

// TODO: performance improvements like moving FFT stuff to @block, frequency limiting

slider1:8<0,9,1{Max,Min (Denoise),Align (FPA) A,Subtract A,Phase limit (unused),Phase reflect,Gap filler 2 ARF,Sustainizer A,Split-Combine Frequency A,Min2 A}>Mode
slider2:100<0,200,1>Amount
slider3:0.5<0,1,0.01>Response
slider4:2<0,4,1{1024,2048,4096,8192}>FFT Size


// import surroundlib3.txt

@init 

fftsize = -1;
pdc_bot_ch = 0;
pdc_top_ch = 4;

// array pointers
window=0;

sixteen_k=16384;


buf0 = sixteen_k*1;
buf1 = sixteen_k*2;

buf2 = sixteen_k*3;
buf3 = sixteen_k*4;

pos=0;

//my arrays for analysis
Magnitudes1=sixteen_k*5;
Magnitudes2=sixteen_k*6;
Magnitudes3=sixteen_k*7;
Magnitudes4=sixteen_k*8;
Magnitudes3_old=sixteen_k*9;
Magnitudes4_old=sixteen_k*10;


Angles1=sixteen_k*11;
Angles2=sixteen_k*12;
Angles3=sixteen_k*13;
Angles4=sixteen_k*14;

Angle_Difference=sixteen_k*31;

old_magnitude=0;

gfx_r = gfx_g = gfx_b = 1;

mode=Amount=0;

left_r = left_i = right_r= right_i= 0;
left_r_two = left_i_two = right_r_two= right_i_two= 0;
	mid_r  = mid_i  = 0;
	side_r = side_i = 0;

tmp=a=b=a2=b2=0;

L_scaling=R_scaling=0;
pi=3.141592653589793238;

ad1=ad2=p1=p2=0;

phaseadj_l = cos_adj_l =sin_adj_l = 0;

inspector1=0;
inspector2=0;

@slider

sliderfft = (2^(slider4+10))|0;
fftsize != sliderfft ?
(
	pos = 0;
	memset(buf0, 0, 2*fftsize);
	memset(buf1, 0, 2*fftsize);
	memset(buf2, 0, 2*fftsize);
	memset(buf3, 0, 2*fftsize);
	memset(window, 0, fftsize);
	fftsize = sliderfft;
	w = 2.0*$pi/fftsize;
	i = 0;
	
	loop(fftsize/2,
		window[i] = 0.42-0.50*cos(i*w)+0.08*cos(2.0*i*w); //can make this a triangle?
		i += 1;
	); 
	pdc_delay = fftsize;
);




Amount=slider2/100;

mode=slider1+1;//cuz I don't like 0-based mode/slider

//V2 Re-Im modes sound different to original Ma-Ang modes

Response=slider3;




@sample

pos >= fftsize ?
(
	tmp = buf0;
	buf0 = buf1;
	buf1 = tmp;

	tmp = buf2;
	buf2 = buf3;
	buf3 = tmp;

	fft(buf0, fftsize);
	fft_permute(buf0, fftsize);

	fft(buf2, fftsize);
	fft_permute(buf2, fftsize);

	i = 2;
	loop(fftsize/2 - 1, 
		
		//loop indexes
		a = i;
		b = a+1;
		a2 = 2*fftsize-i;
		b2 = a2+1;


		//fft() works on complex numbers, but the input is real, so if you put only L in fft, it's wasteful. So they put L & R into Re & Imag to be more efficient. But you need to untangle it.

		x = buf0[a];
		y = buf0[b];
		x2 = buf0[a2];
		y2 = buf0[b2];

		 left_r = (x + x2)*0.5;
		 left_i = (y - y2)*0.5;
		right_r = (y + y2)*0.5;
		right_i = (x2 - x)*0.5;

		x_two = buf2[a];
		y_two = buf2[b];
		x2_two = buf2[a2];
		y2_two = buf2[b2];



		left_r_two = ( x_two + x2_two)*0.5;
		left_i_two = ( y_two - y2_two)*0.5;
		right_r_two= ( y_two + y2_two)*0.5;
		right_i_two= (x2_two -  x_two )*0.5;

		Magnitudes1[i]=sqrt(left_r^2+left_i^2);
		Magnitudes2[i]=sqrt(right_r^2+right_i^2);

		// Angles1[i]=atan2(left_r,left_i);
		// Angles2[i]=atan2(right_r,right_i);

		Magnitudes3[i]=sqrt(left_r_two^2+left_i_two^2);
		Magnitudes4[i]=sqrt(right_r_two^2+right_i_two^2);

		// Angles3[i]=atan2(left_r_two,left_i_two);
		// Angles4[i]=atan2(right_r_two,right_i_two);

		p1=abs(Magnitudes3[i]-Magnitudes1[i])/(Magnitudes3[i]+Magnitudes1[i]); //for FPA & min2

		(mode==1)? //Max mode, sounds ok with no overlap, but a bit of artefact
		(
			//independently test mags of L&R
			(Magnitudes1[i]<Magnitudes3[i])?
			(
			left_r=left_r_two;
			left_i=left_i_two;
			);
			
			(Magnitudes2[i]<Magnitudes4[i])?
			(
			right_r=right_r_two;
			right_i=right_i_two;
			);

		);
		
		(mode==2)? //Min mode
		(
			(Magnitudes1[i]>Magnitudes3[i])?
			(
			left_r=left_r_two;
			left_i=left_i_two;
			);
			
			(Magnitudes2[i]>Magnitudes4[i])?
			(
			right_r=right_r_two;
			right_i=right_i_two;
			);
		);

		(	//FPA mode, align only close magnitudes
			mode==3 &&
			p1 <Amount ;
		)
		? 
		(
			Angles3[i]=atan2(left_r_two,left_i_two);
			Angles4[i]=atan2(right_r_two,right_i_two);

			// Angle_Difference1[i]=Angles3[i]-Angles1[i];
			// Angle_Difference2[i]=Angles4[i]-Angles2[i];

			// Angles1[i]+=Angle_Difference1[i]*slider2/100;
			// Angles2[i]+=Angle_Difference2[i]*slider2/100;

			//just synthesize Re & Im with mag set1 & angles set2
			left_r = Magnitudes1[i]*sin(Angles3[i]) ;
			left_i = Magnitudes1[i]*cos(Angles3[i]) ;
			right_r = Magnitudes2[i]*sin(Angles4[i]) ;
			right_i = Magnitudes2[i]*cos(Angles4[i]) ;

		);

		mode==4? //subtract
		(
			//clamp mag to >0, otherwise will make more sound
			L_scaling=max(Magnitudes1[i]-Magnitudes3[i]*Amount,0)/Magnitudes1[i];
			R_scaling=max(Magnitudes2[i]-Magnitudes4[i]*Amount,0)/Magnitudes2[i];
			
			left_r *= L_scaling;
			left_i *= L_scaling;
			right_r *= R_scaling;
			right_i *= R_scaling;
		);
		
		mode==6? //phase reflect, wider bins sound worse w/ no overlap
		(

			Angles1[i]=atan2(left_r,left_i);
			Angles2[i]=atan2(right_r,right_i);
			
			ad1=abs(Angles1[i]-Angles2[i]);//different to below
			
			ad1>pi?ad1=2*pi-ad1;//fix angle difference cuz should be under 180*
			
			//reverse polarity of R if >90*, but sounds bad
			// ad1>(pi/2)?
			// (
			// right_r *= -1;
			// right_i *= -1;
			// );
			
			ad1>(pi/2)? //if >90*, but sounds good
			(
				//NB: M/S reversed!!! Swap M&S.
				side_r = (left_r+right_r)*0.5;
				side_i = (left_i+right_i)*0.5;
				mid_r  = (left_r-right_r)*0.5;
				mid_i  = (left_i-right_i)*0.5;
				
				left_r  = mid_r + side_r ;
				left_i  = mid_i + side_i ;
				right_r = mid_r - side_r ;
				right_i = mid_i - side_i ;	

			);
				
			
		);
		
		
		(mode==7&&Magnitudes1[i]<old_magnitude)?// Gap filler
		(

			L_scaling=(old_magnitude*Amount)/Magnitudes1[i];
			R_scaling=(old_magnitude*Amount)/Magnitudes2[i];
			
			left_r *= L_scaling;
			left_i *= L_scaling;
			right_r *= R_scaling;
			right_i *= R_scaling;
			
		);
		
		old_magnitude=Magnitudes1[i]*(Response/10)+old_magnitude*(1-Response/10);
		
		(mode==8)? // Sustainizer, does not sound good w/ no overlap
		(

			// left_r*=(-1)^(i/2);
			// right_r*=(-1)^(i/2);
			// left_i*=(-1)^i;
			// right_i*=(-1)^i;
			
			// rand1=ceil(rand(10));
			// rand2=ceil(rand(10));
			
			// left_r*=(-1)^rand1;
			// right_r*=(-1)^rand1;
			// left_i*=(-1)^rand2;
			// right_i*=(-1)^rand2;
			
			//randomise phase, inspired by DTBlkFx
			phaseadj_l = rand(pi)*Amount;
			
			cos_adj_l = cos(phaseadj_l);
			sin_adj_l = sin(phaseadj_l);
			
			//need to put into temporary vars 1st cuz l_i etc get overwritten b4 they're used again= incorrect results
			 left_r_two = left_r* cos_adj_l - left_i* sin_adj_l;
			 left_i_two = left_r* sin_adj_l + left_i* cos_adj_l;
			right_r_two = right_r*cos_adj_l - right_i*sin_adj_l;
			right_i_two = right_r*sin_adj_l + right_i*cos_adj_l;
			
			 left_r =  left_r_two ;
			 left_i =  left_i_two ;
			right_r = right_r_two ;
			right_i = right_i_two ;

		);
		
// inspector1=(i)/fftsize;

		// (mode==9&&((i)/fftsize) > Amount)? //Split-Combine Frequency, freq split (linear scale up to nyquist) depends on Amount
		//Amount controls a log frequency. sqrt of 100 is 10. log can make -ve number so numbers are enlarged 4 comparing
		//does not sound good w/ no overlap. artefacts @ border frequency
		(mode==9&& sqrt((i*100)/fftsize) > Amount*10)?
		(
		
		
			left_r=left_r_two;
			left_i=left_i_two;
			right_r=right_r_two;
			right_i=right_i_two;
		
		);
		
		(mode==10)? //Min2
		(
			Angles1[i]=atan2(left_r,left_i);
			Angles2[i]=atan2(right_r,right_i);
			Angles3[i]=atan2(left_r_two,left_i_two);
			Angles4[i]=atan2(right_r_two,right_i_two);
			
			//compute pans as per my V3 upmix
			// p1=(Magnitudes3[i]-Magnitudes1[i])/(Magnitudes3[i]+Magnitudes1[i]);
			p2=(Magnitudes4[i]-Magnitudes2[i])/(Magnitudes4[i]+Magnitudes2[i]);
		//pan isn't that useful for dual stereo-different noise recordings, but useful for High Touch! (different singers on different versions)
		

			
			ad1=abs(Angles1[i]-Angles3[i]);
			ad2=abs(Angles2[i]-Angles4[i]);

			
			ad1>pi?ad1=2*pi-ad1;//fix angle difference cuz should be under 180*
			ad2>pi?ad2=2*pi-ad2; 

			L_scaling=(1-(min(ad1*2,pi)/pi))*(1-   (p1));//double ad cuz 90* ad means it's already purely stereo noise
			R_scaling=(1-(min(ad2*2,pi)/pi))*(1-abs(p2));//Reduce mags the more it's panned & phased

			left_r *= L_scaling^Amount; //raised to power because doing linear multiple 1- is tricky
			left_i *= L_scaling^Amount;
			right_r *= R_scaling^Amount;
			right_i *= R_scaling^Amount;

		);  
		
		//TODO: center pan, generic de distort via harmonics, excite/pitchshift, denoise, transient designer, Tonal/noise separator, select notes via MIDI input, combine low frequency from set1, high frequency from set2.
		
		
		// convert to cartesian
		//x=R cos (theta)



		buf0[a]  = left_r - right_i;
		buf0[b]  = left_i + right_r;
		buf0[a2] = left_r + right_i;
		buf0[b2] = right_r - left_i;

		i += 2
	);

	fft_ipermute(buf0, fftsize);
	ifft(buf0, fftsize);

	pos=0;
);

//windowing?
w1 = window[pos/2];
w2 = window[(fftsize-pos)/2-1];
sw = (w1+w2)*fftsize;

//i think this is FFT-size inverse scaling
out0 = (buf0[pos]+buf1[fftsize+pos])/sw;
out1 = (buf0[pos+1]+buf1[fftsize+pos+1])/sw;


// out2 = (buf2[pos]+buf1[fftsize+pos])/sw;
// out3 = (buf2[pos+1]+buf1[fftsize+pos+1])/sw;

buf0[pos] = w1*spl0;
buf0[pos+1] = w1*spl1;
buf1[fftsize+pos] = w2*spl0;
buf1[fftsize+pos+1] = w2*spl1;

//more stuff for 2nd input
buf2[pos] = w1*spl2;
buf2[pos+1] = w1*spl3;
buf3[fftsize+pos] = w2*spl2;
buf3[fftsize+pos+1] = w2*spl3;

spl0 = out0;
spl1 = out1;

spl2=spl3=0;//silence set2

pos += 2;

// @gfx 500 16

// box_width=8;

// gfx_a = 0.8; 
// drawbox(floor((pos/fftsize)*gfx_w),8);