%
% algoritmo FFT 1D
%

clear i;  % cosi' siamo sicuri che viene interpretata come unita' immaginaria; potevano anche fare "i = sqrt(-1)"

% sequenza di dati complessi
xi = [1.0+i*1.0 1.3+i*1.8 1.7+i*4.0 2.0+i*1.3 1.0+i*4.0 4.4+i*1.0 1.1+i*3.0 7.0+i*2.0 ...
      2.2+i*2.0 3.3+i*2.2 2.0+i*4.5 1.5+i*6.0 6.0+i*1.4 1.4+i*4.2 1.7+i*1.7 1.3+i*1.2 ];
N = length(xi);

% calcolo diretto della TFD in un secondo vettore "X"
x = xi;
X = zeros(1,N) + i*zeros(1,N);
tic
for k = 0:N-1
  for n = 0:N-1
    Wkn = exp(-i*2*pi/N*k*n);  % ATTENZIONE ALL'USO DI "i" COME VARIABILE !!!
    X(k+1) = X(k+1) + x(n+1)*Wkn;
    % equivalente, mostrando per esteso la moltiplicazione tra valori complessi:
    %X(k+1) = X(k+1) + real(x(n+1))*real(Wkn) - imag(x(n+1))*imag(Wkn) + i*( real(x(n+1))*imag(Wkn) + imag(x(n+1))*real(Wkn) );
  end;
end;
texec_diretto = toc;
%disp(['tempo di esecuzione calcolo TFD con algoritmo "diretto": ' num2str(texec_diretto) ' [sec]']);
disp(['verifica: ' num2str(sum(abs(real(X - fft(xi)))) + sum(abs(imag(X - fft(xi)))))]);

%
% calcolo della FFT
%
x = xi;
%
% passo 1: permutazione dei campioni della sequenza secondo l'ordinamento a disposizione invertita dei bit
%
j=0;
tic
for k=0:N-1
  disp(['k=' dec2bin(k,log2(N)) '   j=' dec2bin(j,log2(N))]);
  if j > k   % faccio la permutazione delle componenti
    temp = x(j+1); % "n+1" e "j+1" perchè in Matlab i vettori cominciano dall'indice "1" e non "0"
    x(j+1) = x(k+1);
    x(k+1) = temp;
  end;
  %
  % calcolo il prossimo "j" secondo l'ordinamento a disposizione invertita dei bit
  n = N/2;
  while n <= j
    j = j - n;
    n = n / 2;
  end;
  j = j + n;
end;
%
% passo 2: calcolo "sul posto" della sequenza trasformata  NB: qui "i" e "j" sono usate come variabili !
%
%tic  % NB: scommentare quest'istruzione per vedere quanto incide il riordinamento degli indici sul tempo di esecuzione !
for k=1:log2(N)  % "k" è l'indice dello stadio di calcolo
  u = complex(1,0);  %  W_N^0
  w = complex(cos(pi/2^(k-1)),-sin(pi/2^(k-1))); 
  for j=1:2^(k-1) 
    for ip=j:2^k:N-1
      % calcolo della "farfalla":
      iq = ip + 2^(k-1);
      %disp(['i=' num2str(i) '  ip=' num2str(ip) '  u=' num2str(u)]);
      temp = x(iq) * u;
      x(iq) = x(ip) - temp;
      x(ip) = x(ip) + temp;
    end;
    u = u * w;
  end;
end;
texec_FFT = toc;
%disp(['tempo di esecuzione calcolo TFD con algoritmo "FFT": ' num2str(texec_FFT) ' [sec]']);
disp(['verifica: ' num2str(sum(abs(real(x - fft(xi)))) + sum(abs(imag(x - fft(xi)))))]);

