% Copyright Mobile Satellite Services Corporation 2001
% file: FFT2KSIM.M
% Author: Alain Zarembowitch / MSS
% Last edited: 12-5-05
% Revision: 2
% FFT simulation
%
% Rev 2. 12/5/05 AZ
% Corrected spectral inversion.
%--------------------------------------------------
% place reference signal here
f = 1/200;
% X = complex(zeros(2048,1),zeros(2048,1));
% X(2) = complex(1.0,0.0);

% reference
Yref = fft(X, 2048);

% 1st stage = radix 2
for k = 1:1024
   Wn = cos(-2*pi*(k-1)/2048) + j*sin(-2*pi*(k-1)/2048);
   X1(k) = X(k) + X(k+1024);
   X1(k+1024) = (X(k) - X(k+1024))*Wn;
end   


% 2nd stage = radix 4, 2 parts, 1024 ppt FFT
for l = 0:1
for k = 1+l*1024:256+l*1024
   Cb = cos(2*pi*(k-1)/1024);
   Sb = sin(2*pi*(k-1)/1024);
   Cc = cos(2*pi*2*(k-1)/1024);
   Sc = sin(2*pi*2*(k-1)/1024);
   Cd = cos(2*pi*3*(k-1)/1024);
   Sd = sin(2*pi*3*(k-1)/1024);
   xa = real(X1(k));
   ya = imag(X1(k));
   xb = real(X1(k+256));
   yb = imag(X1(k+256));
   xc = real(X1(k+512));
   yc = imag(X1(k+512));
   xd = real(X1(k+768));
   yd = imag(X1(k+768));
   
   x2a = xa + xb + xc + xd;
   y2a = ya + yb + yc + yd;
	x2b = (xa+yb-xc -yd)*Cb - (ya-xb-yc+xd)*(-Sb);   
	y2b = (ya-xb-yc+xd)*Cb + (xa+yb-xc-yd)*(-Sb);   
	x2c = (xa-xb+xc-xd)*Cc - (ya-yb+yc-yd)*(-Sc);   
	y2c = (ya-yb+yc-yd)*Cc + (xa-xb+xc-xd)*(-Sc);   
	x2d = (xa-yb-xc+yd)*Cd - (ya+xb-yc-xd)*(-Sd);   
	y2d = (ya+xb-yc-xd)*Cd + (xa-yb-xc+yd)*(-Sd);   
   
   X2(k) = x2a + j*y2a;
   X2(k+256) = x2b + j*y2b;
   X2(k+512) = x2c + j*y2c;
   X2(k+768) = x2d + j*y2d;
end   
end

% 3rd stage = radix 4, 4 parts, 256 ppt FFT
for l = 0:7
for k = 1+l*256:64+l*256
   Cb = cos(2*pi*(k-1)/256);
   Sb = sin(2*pi*(k-1)/256);
   Cc = cos(2*pi*2*(k-1)/256);
   Sc = sin(2*pi*2*(k-1)/256);
   Cd = cos(2*pi*3*(k-1)/256);
   Sd = sin(2*pi*3*(k-1)/256);
   xa = real(X2(k));
   ya = imag(X2(k));
   xb = real(X2(k+64));
   yb = imag(X2(k+64));
   xc = real(X2(k+128));
   yc = imag(X2(k+128));
   xd = real(X2(k+192));
   yd = imag(X2(k+192));
   
   x2a = xa + xb + xc + xd;
   y2a = ya + yb + yc + yd;
	x2b = (xa+yb-xc -yd)*Cb - (ya-xb-yc+xd)*(-Sb);   
	y2b = (ya-xb-yc+xd)*Cb + (xa+yb-xc-yd)*(-Sb);   
	x2c = (xa-xb+xc-xd)*Cc - (ya-yb+yc-yd)*(-Sc);   
	y2c = (ya-yb+yc-yd)*Cc + (xa-xb+xc-xd)*(-Sc);   
	x2d = (xa-yb-xc+yd)*Cd - (ya+xb-yc-xd)*(-Sd);   
	y2d = (ya+xb-yc-xd)*Cd + (xa-yb-xc+yd)*(-Sd);   
   
   X3(k) = x2a + j*y2a;
   X3(k+64) = x2b + j*y2b;
   X3(k+128) = x2c + j*y2c;
   X3(k+192) = x2d + j*y2d;
end   
end

% 4rd stage = radix 4, 16 parts, 64 ppt FFT
for l = 0:31
for k = 1+l*64:16+l*64
   Cb = cos(2*pi*(k-1)/64);
   Sb = sin(2*pi*(k-1)/64);
   Cc = cos(2*pi*2*(k-1)/64);
   Sc = sin(2*pi*2*(k-1)/64);
   Cd = cos(2*pi*3*(k-1)/64);
   Sd = sin(2*pi*3*(k-1)/64);
   xa = real(X3(k));
   ya = imag(X3(k));
   xb = real(X3(k+16));
   yb = imag(X3(k+16));
   xc = real(X3(k+32));
   yc = imag(X3(k+32));
   xd = real(X3(k+48));
   yd = imag(X3(k+48));
   
   x2a = xa + xb + xc + xd;
   y2a = ya + yb + yc + yd;
	x2b = (xa+yb-xc -yd)*Cb - (ya-xb-yc+xd)*(-Sb);   
	y2b = (ya-xb-yc+xd)*Cb + (xa+yb-xc-yd)*(-Sb);   
	x2c = (xa-xb+xc-xd)*Cc - (ya-yb+yc-yd)*(-Sc);   
	y2c = (ya-yb+yc-yd)*Cc + (xa-xb+xc-xd)*(-Sc);   
	x2d = (xa-yb-xc+yd)*Cd - (ya+xb-yc-xd)*(-Sd);   
	y2d = (ya+xb-yc-xd)*Cd + (xa-yb-xc+yd)*(-Sd);   
   
   X4(k) = x2a + j*y2a;
   X4(k+16) = x2b + j*y2b;
   X4(k+32) = x2c + j*y2c;
   X4(k+48) = x2d + j*y2d;
end   
end



% 5th stage = radix 4, 64 parts, 16 ppt FFT
for l = 0:127
for k = 1+l*16:4+l*16
   Cb = cos(2*pi*(k-1)/16);
   Sb = sin(2*pi*(k-1)/16);
   Cc = cos(2*pi*2*(k-1)/16);
   Sc = sin(2*pi*2*(k-1)/16);
   Cd = cos(2*pi*3*(k-1)/16);
   Sd = sin(2*pi*3*(k-1)/16);
   xa = real(X4(k));
   ya = imag(X4(k));
   xb = real(X4(k+4));
   yb = imag(X4(k+4));
   xc = real(X4(k+8));
   yc = imag(X4(k+8));
   xd = real(X4(k+12));
   yd = imag(X4(k+12));
   
   x2a = xa + xb + xc + xd;
   y2a = ya + yb + yc + yd;
	x2b = (xa+yb-xc -yd)*Cb - (ya-xb-yc+xd)*(-Sb);   
	y2b = (ya-xb-yc+xd)*Cb + (xa+yb-xc-yd)*(-Sb);   
	x2c = (xa-xb+xc-xd)*Cc - (ya-yb+yc-yd)*(-Sc);   
	y2c = (ya-yb+yc-yd)*Cc + (xa-xb+xc-xd)*(-Sc);   
	x2d = (xa-yb-xc+yd)*Cd - (ya+xb-yc-xd)*(-Sd);   
	y2d = (ya+xb-yc-xd)*Cd + (xa-yb-xc+yd)*(-Sd);   
   
   X5(k) = x2a + j*y2a;
   X5(k+4) = x2b + j*y2b;
   X5(k+8) = x2c + j*y2c;
   X5(k+12) = x2d + j*y2d;
end   
end

% 6th stage = radix 4, 256 parts, 4 ppt FFT
for l = 0:511
for k = 1+l*4:1+l*4
   Cb = cos(2*pi*(k-1)/4);
   Sb = sin(2*pi*(k-1)/4);
   Cc = cos(2*pi*2*(k-1)/4);
   Sc = sin(2*pi*2*(k-1)/4);
   Cd = cos(2*pi*3*(k-1)/4);
   Sd = sin(2*pi*3*(k-1)/4);
   xa = real(X5(k));
   ya = imag(X5(k));
   xb = real(X5(k+1));
   yb = imag(X5(k+1));
   xc = real(X5(k+2));
   yc = imag(X5(k+2));
   xd = real(X5(k+3));
   yd = imag(X5(k+3));
   
   x2a = xa + xb + xc + xd;
   y2a = ya + yb + yc + yd;
	x2b = (xa+yb-xc -yd)*Cb - (ya-xb-yc+xd)*(-Sb);   
	y2b = (ya-xb-yc+xd)*Cb + (xa+yb-xc-yd)*(-Sb);   
	x2c = (xa-xb+xc-xd)*Cc - (ya-yb+yc-yd)*(-Sc);   
	y2c = (ya-yb+yc-yd)*Cc + (xa-xb+xc-xd)*(-Sc);   
	x2d = (xa-yb-xc+yd)*Cd - (ya+xb-yc-xd)*(-Sd);   
	y2d = (ya+xb-yc-xd)*Cd + (xa-yb-xc+yd)*(-Sd);   
   
   X6(k) = x2a + j*y2a;
   X6(k+1) = x2b + j*y2b;
   X6(k+2) = x2c + j*y2c;
   X6(k+3) = x2d + j*y2d;
end   
end

% reshuffle indices
for k = 1:2048
   input_index = dec2bin(k-1,11);
   output_index(1) = input_index(10);
   output_index(2) = input_index(11);
   output_index(3) = input_index(8);
   output_index(4) = input_index(9);
   output_index(5) = input_index(6);
   output_index(6) = input_index(7);
   output_index(7) = input_index(4);
   output_index(8) = input_index(5);
   output_index(9) = input_index(2);
   output_index(10) = input_index(3);
   output_index(11) = input_index(1);
   k2 = bin2dec(output_index) + 1;
   X7(k2) = complex(real(X6(k)), -imag(X6(k))); % spectral inversion
end   

% compare Matlab FFT with this radix2/radix4 FFT
plot(abs(X7));
hold;
plot(abs(Yref),'r');

