% INICIALIZA��O
% N�mero de est�gios. Este par�metro define todos os outros na implementa��o da FFT.
est=4;
N=2^est

% vector de teste para o caso de uma rampa de magnitude normalizada.
re=1/N*([0:N-1]);

im=zeros(1,N);



% vector de teste para uma sinus�ide pura � entrada da FFT.
% (comentar se se preferir o sinal de rampa)
bin_teste=2  % 0 <= bin_teste <= N/2
re=sin(bin_teste*2*pi/N*[0:N-1]); im=zeros(1,N);

%
% prepara��o do vector que define as trocas do "bit reversal"
% (n�o � necess�rio analisar em detalhe, nem preciso responder �s 3 perguntas)
%
nchg=0; N2=N/2; flag=zeros(1,N);
for i=1:N-2,                         % excluem-se os casos de i==0 e i==N-1, porqu� ?
   if(flag(i)==0)
      dv=N2; ic=1; final=0; pos=i;
      for k=1:est,                   % que faz este ciclo "for" ?
         if (pos/dv >= 1)
            pos=pos-dv;
            final=final+ic;
         end
         dv=dv/2; ic=ic*2;
      end
      if (i~=final)                  % que se pretende evitar aqui ?
         rev(1+nchg)=i;
         rev(N2+nchg)=final;
         flag(final)=1;
         nchg=nchg+1;
      end
   end   
end

% este � o n�mero de trocas a fazer
nchg

% este � o vector que detalha as trocas a fazer
rev(1:nchg)        % pontos de um lado
rev(N2:N2+nchg-1)  % pontos do outro lado ("bit reversed")


%re   % antes da FFT
%im   % antes da FFT

grup=1; butf=N; phi=2*pi/N;
for i=0:est-1,                            % avan�a nos est�gios
   butf=butf/2;
   for j=0:grup-1,                        % avan�a nos grupos
      for k=0:butf-1,                     % avan�a nas borboletas
         ptr1=2*j*butf+k;
         ptr2=ptr1+butf;
         arg=phi*grup*k;
         
         retmp = re(1+ptr1) - re(1+ptr2);
         imtmp = im(1+ptr1) - im(1+ptr2);
         
         re(1+ptr1)=re(1+ptr1)+re(1+ptr2);
         im(1+ptr1)=im(1+ptr1)+im(1+ptr2);
         
         C=cos(arg);
         S=sin(arg);
         
         re(1+ptr2) = retmp*C + imtmp*S;
         im(1+ptr2) = imtmp*C - retmp*S;
       
      end
  end
    grup=grup*2;
end

%
% "bit reversal"
%
for i=1:nchg,
   tmp=re(1+rev(i)); re(1+rev(i))=re(1+rev(N2+i-1)); re(1+rev(N2+i-1))=tmp;
   tmp=im(1+rev(i)); im(1+rev(i))=im(1+rev(N2+i-1)); im(1+rev(N2+i-1))=tmp;
end



%re  % depois da FFT
%im  % depois da FFT

y=(re.^2+im.^2)/N^2;
plot([0:N-1],y);
title('Densidade Espectral com FFT DIF');
xlabel('Linha Espectral');
ylabel('Quadrado da Magnitude');