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

N = size( a, 1 );

if ( abs( round( log( N )/log(3) ) - log( N )/log(3) ) > 10*eps )
    % Error if n is not a power of 3.
    error('I need a vector of length a power of 3. Sorry.')
end

% If n < 27, then do not recurse:
if ( log( N )/log(3) < 3 )
    A = fft( a );
    return
end

%%%  REMEMBER THE 3 FACTS FROM CLASS  %%%
% 
%  We need FACTS 1-3.  FACT 1 is independent of N. 
% 
%%%   FACT 2:  %%%
%
%        If N is a multiple of 3, then you can decompose A(z) into 3 polynomials
%
%       A(z) = a_0z^0                 a_3z^3                  a_6z^6
%                      a_1z^1                 a_4z^4   
%                             a_2z^2                   a_5z^5
%   
%   =>  A(z) =   A_1(z^3) + zA_2(z^3) + z^2A_3(z^3), (Check the degrees are consistent.)
%
%
%%%   FACT 3: %%%
%
%       If N is a multiple of 3, the N-roots of unity enjoy the remarkable
%       CUBING PROPERTY, i.e., the N-roots of unity cubed are the N/3-roots of
%       unity:
%
%       N = 27; 
%       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/3+1)'; theta(end)=[];
%       zNthird = exp( -1i*theta );   % N/2-roots of unity in reverse order.
%        
%       zN.^3 - [ zNthird ; zNthird ; zNthird ]  % = 0 (up to machine precision)
%
%
%%%  BONUS FACT:  %%%
%
%       If N is a multiple of 3, then we have    
%
%         w = exp(-1i*2*pi/3);            % [1 w w^2] are the 3-roots of unity
%         zN(1:N/3) - w*zN(N/3+1:2*N/3)                      (in reverse order)
%         zN(1:N/3) - w^2*zN(2*N/3+1:N)   
%  
%       (Draw a diagram if you cannot remember why these relations hold.)
%
%
%  ALGORITHM:  JUST USE THE FACTS ABOVE... YOU'LL NEED THE BONUS FACT!
%
%     A(zN) =    [ I    I       I  ]  [     A_1( zNthird )     ]
%                [ I    wI     w^2I]  [ zN(1:N/3)*A_2(zNthird) ]
%                [ I   w^2I    w^4I]  [ zN(1:N/3)^2*A_3(zNthird) ]
%

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

I = eye(N/3);
w = exp(-1i*2*pi/3);
% Evaluate A(zN):
A = [I I I ; I w*I w^2*I; I w^2*I w^4*I] * [A1 ; zN(1:N/3).*A2 ; zN(1:N/3).^2.*A3]; 

% Equivalently... This is what Cooley and Tukey did not realize!: 
% A = fft([A1 zN(1:N/3).*A2 zN(1:N/3).^2.*A3].').'; 
% A = A(:);
end