function [X] = dsp_fft(x,N,c,rind)
% X = dsp_fft(x,N)
% Compute X fft coeficients based on discrete input x
% If N is not power of two, increase N to power of 2
% Then pad x with zeros or truncate as N determines
% inputs:
% x    -> input data vector
% N    -> length of fft
% c    -> pre-calculated coef. for fft
% rind -> bit-reverse index for input vector
% outputs:
% X -> fft of x, length(X) = N
% dsp_fft.m
% author:    Nicholas Kottenstette
% date:      11/9/2004
% copyright: 2004 by Nicholas Kottenstette
  
% performance:
% N = 2^16;
% x = sin(0.25*pi*[0:N-1]);
% baseline to compare:
% tic;fft(x,N);toc            -> 0.006 seconds
% tic;dsp_fft(x,N,[],[]);toc  -> 0.086 seconds
% tic;dsp_fft(x,N,c,rind);toc -> 0.048 seconds
% Note when N = 2^15
% tic;dsp_fft(x,N,c,rind);toc -> 0.017 seconds 
% Note when N = 2^14
% tic;dsp_fft(x,N,c,rind);toc -> 0.006 seconds
% Note when N = 2^13
% tic;dsp_fft(x,N,c,rind);toc -> 0.0025 seconds
% Note when N = 2^12
% tic;dsp_fft(x,N,c,rind);toc -> 0.0012 seconds
% tic;fft(x,N);toc            -> 0.0002 seconds
  [m,n] = size(x);
  % Check x is a vecor and will keep as a single column vector
  if (m > 1) & (n > 1)
    error('dsp_fft:invalid_input','x not a vector m = %d, n = %d\n"',m,n);
  elseif m < n
    x = x';
    [m,n] = size(x);
  end
  if N < 2
    N = m;
  end
  N = bitshift(1,ceil(log2(N)));
  Nx = m;
  % zero pad or truncate x
  if Nx < N
    x = [x;zeros(N-Nx,1)];
  elseif Nx > N
    x = x(1:N,1);
  end
  % Shuffle inputs
  [xs,Smax] = shuffle(x,N,rind);
  % Calculate fft using sparse matrices
  X = mfft(xs,Smax,c);
  
function [xs,Smax] = shuffle(x,N,rind)
  Smax = log2(N);
  if length(rind) ~= N
    fp = fopen('fft_coef/rind.mat');
    if fp == -1
      for smax = 1:16
        ind = [0:2^smax-1]';
        rit{smax} = bitrevorder(ind);
      end
      save fft_coef/rind.mat rit
      rind = rit{Smax};
    else
      fclose(fp);
      fp = load('fft_coef/rind.mat');
      rind = fp.rit{Smax};
    end
  end
  xs = x(rind+1);
  
function X = mfft(x,Smax,c)
  if length(c) ~= 2^16 - 1
    fp = fopen('fft_coef/coef.mat');
    if fp == -1
      c = gen_coef();
    else
      fclose(fp);
      fp = load('fft_coef/coef.mat');
      c = fp.c;
    end
  end
  X = fft_mult_im(c,x,Smax);
  
function c = gen_coef()
  SMAX = 16;
  NcMax = 2^SMAX;
  ind = 1;
  for s = 1:SMAX
    L = 2^s;
    for i = 0:2^(s-1)-1
      c(ind,1) = exp(-j*2*pi*i/L);
      ind = ind + 1;
    end
  end
  save fft_coef/coef.mat c
  