%----------------------------------------------------------------
%----------------------------------------------------------------
%exa_fft16.m --- demonstrate how FFT is done for a size 16 vector
%                Example of direct FFT :  F = A_16 * f
%----------------------------------------------------------------
if ~exist('mygrid'), disp('Need M-file mygrid.m for plotting ...');return;end
rand('seed',0), format compact, help exafft16
    n=16
if ~exist('f'), f = round( rand(n,1)*n ); end
                 if length(f)~=n, disp('f has the wrong size ...');return;end
                 f=f(:); clear i
disp('Direct DFT ------------------------------')
A_16=[]; omega=exp(-2*pi*i/n)
for j=1:n
 for k=1:n
 A_16(j,k)=omega^( (j-1)*(k-1) );
 end
end
 F_16 = A_16 * f;
 disp('F_16 ='); mat_prt4(F_16)

disp('Step-by-Step FFT ------------------------')
r = [ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 ];
rr = [0:2:15 1:2:15], r0=1+[0:2:7 1:2:7]; % Round 1 permutation
      r2=rr(9:16); r1=rr(1:8);
r8 = [r1(r0) r2(r0)], r0=1+[0:2:3 1:2:3]; % Round 2 permutation
      r4=r8(13:16); r3=r8(9:12);
      r1=r8(1:4);   r2=r8(5:8);
r16 = [r1(r0) r2(r0) r3(r0) r4(r0)]        % Round 3 permutation
     I = eye(n); PT = I(:,r16+1)'; P_16=PT';
                   close all, figure; spy(PT,'x',12), 
                   title('exafft16.m --- P^T matrix')
                   set(gca,'Xtick',0:n), set(gca,'Ytick',0:n)
           
disp(['-----------Here is fr=P^T*f with the permutation matrix '...
      'P^T based on r16---------'])
fr = PT*f;  
                    f_orig=f', f_perm=fr'

%............... Step 1 >>> B_{2^1}
 d = exp( -2 *pi*i * (0:0) / 2^1 );     b2=[ eye(1) diag(d); eye(1) -diag(d)];
B2 = kron( eye(8), b2); f2 = B2*fr;
                     figure; spy(B2), nx=2; mygrid
                     title('exafft16 --- B_2 using \Omega_1')
                     B2_matrix = b2, f2_vector=f2' % Display only

%............... Step 2 >>> B_{2^2}
 d = exp( -2 *pi*i * (0:2^1-1) / 2^2 ); b4=[ eye(2) diag(d); eye(2) -diag(d)];
B4 = kron( eye(4), b4); f4 = B4*f2;
                     figure; spy(B4), nx=4; mygrid
                     title('exafft16 --- B_4 using \Omega_2')
                     B4_matrix = b4, f4_vector=f4' % Display only

%............... Step 3 >>> B_{2^3}
 d = exp( -2 *pi*i * (0:2^2-1) / 2^3 ); b8=[ eye(4) diag(d); eye(4) -diag(d)];
B8 = kron( eye(2), b8); f8 = B8*f4;
                     figure; spy(B8), nx=8; mygrid
                     title('exafft16 --- B_8 using \Omega_4')
                     B8_matrix = b8, f8_vector=f8' % Display only

%............... Step 4 >>> B_{2^4}
 d = exp( -2 *pi*i * (0:2^3-1) / 2^4 ); b1=[ eye(8) diag(d); eye(8) -diag(d)];
B16 = kron( eye(1), b1); f16 = B16*f8;
                     figure; spy(B16), nx=16; mygrid
                     title('exafft16 --- B_{16} using \Omega_8')
                     B16_matrix = b1, f16_vector=f16' % Display only
disp('____________________|__________________|__________________')
  f_fft_compare = [F_16 f16 fft(f)];
  mat_prt4( f_fft_compare )

disp('Direct A_16*f_______|Step FFT__________|Matlab FFT________')
disp('_____FFT Theorem checking____Notation:  Bn denotes I_{N/n} x B_n')
disp('                                        B16, B8, B4, B2, r16')
fprintf(' norm( A_16*P_16 - B16*B8*B4*B2 )\n') 
