function A = CooleyTukey( a )
% A VERY simple implementation of the Cooley & Tukey algorithm
% for computing the FFT. Only works if length(v) is a power of 2.
%
% Alex Townsend, August 2014. Updated September 2014.

N = size( a, 1 );

if ( abs( round( log2( N ) ) - log2( N ) ) )
    % Error if n is not a power of 2.
    error('I need a vector of length a power of 2. Sorry.')
end

% If n < 16, then do not recurse:
if ( log2( N ) < 4 )
    A = fft( a );
    return
end

%%%  REMEMBER THE 3 FACTS FROM CLASS  %%%
%
%%%   FACT 1:  %%%
%
%         You can think of the DFT as polynomial evaluation of a degree N-1 at
%         the N roots of unity (in reverse order):
%
%
%         X_k = A(z_k) = sum_n  a_n z_k^n,    z_k = e^(-2*pi*1i*j*n/N)
%
%        (This fact is purely to make the DFT easy to remember.)
%
%
%
%%%   FACT 2:  %%%
%
%        If N is even, then you can decompose A(z) into two polynomials
%
%       A(z) = a_0z^0         a_2z^2           a_4z^4         
%                      a_1z^1          a_3z^3          a_5z^5
%   
%   =>  A(z) =   A_1(z^2)  +   zA_2(z^2),   (Check the degrees are consistent.)
%
%       (You'll remember Fact 2 if N is a power of 3.)
%
%
%
%%%   FACT 3: %%%
%
%       If N is even, the N-roots of unity enjoy the remarkable SQUARING
%       PROPERTY, i.e., the N-roots of unity squared are the N/2-roots of
%       unity:
%
%       N = 8;
%       theta = linspace(0, 2*pi, N+1)'; theta(end)=[];
%       zN = exp( -1i*theta );       % N-roots of unity in reverse order.
%       theta = linspace(0, 2*pi, N/2+1)'; theta(end)=[];
%       zNhalf = exp( -1i*theta );   % N/2-roots of unity in reverse order.
%          
%       zN.^2 - [zNhalf ; zNhalf]    % SQUARING PROPERTY (why are they
%                                                               repeated twice?)
%
%
%%%  BONUS FACT:  %%%
%
%       If N is even, the first N/2 of the N-roots of unity equal minus the last
%       N/2. That is,
%
%       zN(1:N/2) - (-zN(N/2+1:N))    % = 0 (up to machine precision)
%
%       (Do you remember the bonus fact for N = power of 3?)
%
%
%
%  ALGORITHM:  JUST USE THE FACTS ABOVE... YOU'LL NEED THE BONUS FACT!
%
%     A(zN) =    [ I   I ] [     A_1( zNhalf )     ]
%                [ I  -I ] [ zN(1:N/2)*A_2(zNhalf) ]
%

theta = linspace(0,2*pi,N+1)'; theta(end)=[];
zN = exp(-1i*theta);                    % N roots of unity (reverse order)
A1 = CooleyTukey( a( 1:2:N ) );         % evaluate A_1(zNhalf)
A2 = CooleyTukey( a( 2:2:N ) );         % evaluate A_2(zNhalf)

I = eye(N/2);
A = [I I ; I -I] * [A1 ; zN(1:N/2).*A2];   % evaluate A(zN)

end