function x = fft_dit(x)
% FFT_DIT decimation in time fast fourier transform
% x = fft_dit(x)
%
% Compute fourier transform of x by the decimation in time FFT.
%
% fft(x) = fft_dit(bitrevorder(x))
%
% See also fft_dif.
%

p = nextpow2(length(x));
Np = 2^p;
if length(x) ~= Np
        error('Input vector length is not a power of two.');
end
T = zeros(Np/2, 1);
for n = 0:Np/2-1
        T(n+1) = exp(-1i*2*pi*n/Np);
end

B = Np/2;
step_T = Np/2;
N = 2;

for P = 1:p
        base_t = 0;
        for b = 1:B
                base_b = base_t + N/2;
                nT = 1;
                for k = 1:N/2
                        t = x(base_t + k);
                        b = T(nT) * x(base_b + k);
                        x(base_t + k) = t + b;
                        x(base_b + k) = t - b;
                        nT = nT + step_T;
                end
                base_t = base_t + N;
        end
        B = B/2;
        N = 2*N;
        step_T = step_T/2;
end