% fft_anim.m  --- 
% animate the Danielson Lanczos FFT algorithm
% HP Gavin, Duke University, 2018-09-21


RealSignal = 1;  % 1: freqs match f_q, 2: freq's do not match, 3: unit impluse

nn  = 128;
dt  = 1;
t1 = [0:nn-1]*dt;
t(1:2:2*nn-1) = t1;
t(2:2:2*nn)   = t1;
f_max = 0.5/dt;
df = 1/max(t1);
f = [ 0:nn/2 , -nn/2+1:-1] *df;

switch RealSignal
  case 1
    datar = sin(2*pi*f_max/8*t1) + cos(2*pi*f_max/4*t1);    % aligned freq's
    title_str = 'signal frequencies aligned with FFT frequencies';
  case 2
    datar = sin(2*pi*f_max/7*t1) + cos(2*pi*f_max/3*t1);    % not aligned
    title_str = 'signal frequencies not aligned with FFT frequencies';
  case 3
    datar = exp(-0.1*t1).*sin(2*pi*f_max/4*t1);             % unit impulse
    title_str = 'impulse response';
end

datai = zeros(1,nn);     % imaginary data

data = zeros(1,2*nn);
data(1:2:2*nn-1) = datar;
data(2:2:2*nn)   = datai;

fft_matlab = fft(datar + sqrt(-1)*datai);

formatPlot(18,2,3);
figure(1)
  clf
  hold on
  hr = plot(t1,data(1:2:2*nn-1),'-or');
  hi = plot(t1,data(2:2:2*nn)  ,'-ov');
  hx = xlabel('time, s');
  legend('real','imaginary')
  legend('boxoff')
  legend('left')
  title(title_str);

pause(5)

% ======================================================================
% Danielson-Lanczos FFT algorithm 
% translated from routine four1.c 
% Press et al., Numerical Recipes in C, Cambridge University Press, 1988, Ch. 12
% https://www.cec.uchile.cl/cinetica/pcordero/MC_libros/NumericalRecipesinC.pdf
% ======================================================================

isign = -1;          % forward FFT :  time to frequency

n = length(data);    % must be a power of two for this algorithm
j=1;

set(hx,'string','index')
for i=1:2:n-1           % bit reversal
  if (j > i)            % exchange two complex numbers
    temp = data(i);   data(i)   = data(j);   data(j)   = temp; % swap real parts
    temp = data(i+1); data(i+1) = data(j+1); data(j+1) = temp; % swap imag parts
    set (hr,'Ydata',data(1:2:2*nn-1)); 
    set (hi,'Ydata',data(2:2:2*nn));  
    title('bit reversal in progress : exchanging complex values')
    drawnow;
    pause(0.10)
  end
  m = n/2;              % bit operator    % m = n >> 1
  while (m >= 2 && j > m) 
    j = j-m;
    m = m / 2;          % bit operator    % m = m >>= 1;
  end
  j = j+m;
end

% Danielson-Lanczos algorithm
mmax = 2; 
while (n > mmax)         % outer loop evaluated log2(nn) times
  istep = mmax * 2;      % bit operator    % mmax << 1;
  theta = isign * (2*pi / mmax);
  wtemp = sin(0.5*theta); 
  wpr   = -2.0*wtemp*wtemp;
  wpi   =  sin(theta);
  wr   = 1.0;
  wi   = 0.0;
  title('Danielson-Lanczos in progress')
  for m=1:2:mmax-1
    for i=m:istep:n         % Danielson-Lanczos formula
      j         = i + mmax;  
      tempr     = wr*data(j)   - wi*data(j+1);
      tempi     = wr*data(j+1) + wi*data(j);
      data(j)   = data(i)   - tempr;
      data(j+1) = data(i+1) - tempi;
      data(i)   = data(i)   + tempr;
      data(i+1) = data(i+1) + tempi;
      set (hr,'Ydata',data(1:2:2*nn-1)); 
      set (hi,'Ydata',data(2:2:2*nn));  
      drawnow;
      pause(0.10)
    end                   
    wtemp = wr;    % trigonometric recurrence 
    wr = wtemp*wpr - wi*wpi + wr;
    wi = wi*wpr + wtemp*wpi + wi;
  end 
  mmax = istep;
end

pause(1)

clf
hold on
plot(f(1:nn/2+1),data(1:2:nn+2),'-ob');
plot(f(1:nn/2+1),data(2:2:nn+2),'-or');
plot(f(nn/2+2:nn),data(nn+3:2:2*nn),'-ob');
plot(f(nn/2+2:nn),data(nn+4:2:2*nn),'-or');
plot(f,real(fft_matlab),'xp', f,imag(fft_matlab),'xp') % compare to Matlab fft
title(title_str);
if ( RealSignal == 3 ) , title('frequency response'); end
legend('real    (even)','imaginary (odd) ')
legend('boxoff')
legend('left')
xlabel('frequency, Hz')
axis('tight')
drawnow;

