function [Xk] = slow_fft(xn,N)
% "Slow" FFT
% Computes DIF Fast Fourier Transform but without
% using MATLAB built-in efficiencies
%  
% [Xk] = slow_fft(xn,N)
% Xk = FFT coeff. array over 0 <= k <= N-1
% xn = N-point finite-duration sequence
% N = Length of DFT
% Doesn't run on MATLAB V4.2 because of multiple function calls in same file.

% Coded by Terrence Sim, 10/97 [thanks, Terrence!]

% Comments and intro by Terrence:
%	I thought it might be instructive to look at actual MATLAB code that
%	implements FFT.  The ones given in OS2 are in FORTRAN, which is perhaps
%	alien to many of us.
%	
%	The code below is a simple straightforward rendition of the
%	decimation-in-freq algorithm presented in class.  It assumes that N is a
%	power of 2.  The recursive structure is apparent.  Of course, one can
%	optimize it further by precomputing the twiddle factors (the WN^n).

%	Feel free to copy and use or comment or improve.

% FFT
% assumes N is power of 2 
% Xk: DFT coefficients of the N-point DFT of xn
% xn: input signal

if rem(log2(N),1) ~=0
	disp('slow_fft: N must be an exact power of 2')
	return
	end

if N == 2
% base case, do 1 butterfly
 [even, odd] = butterfly(xn(1,1), xn(1,2));
  Xk = [even, odd];
else
  n2 = N/2;
  [even,odd] = butterfly(xn(1,1:n2), xn(1,(n2+1):N));
  % Make vector of wn^0 to wn^(n2-1)
  Wn = make_twiddle_vector(N);
  odd = odd .* Wn;
  Xeven = slow_fft(even, n2); % recursive call to compute N/2 DFT
  Xodd = slow_fft(odd, n2); % recursive call to compute N/2 DFT
  % now merge the 2 vectors to get sorted output
  Xk = merge(Xeven, Xodd);
end

%% Defintions of helper routines
function [sum, diff] = butterfly(x1, x2)
sum = x1 + x2;
diff = x1 - x2;

function v = make_twiddle_vector(n)
wn = exp(-j*2*pi/n);
v = 0:1:(n/2 - 1);
v = wn .^ v;

function v3 = merge(v1, v2)
v3 = [];
for i=1:length(v1)
  v3 = [v3, v1(i), v2(i)];
end