% Radix-2 DIT implementation in GNU Octave, version 3.4.3
% John Bryan, December 2016
echo off;
clear all;


function x=bitreversal(x,N)
    M1=2;
    r(1)=0;
    r(2)=1;
    % initialize r() in while loop
    while  (M1 < N)
         k = 0;
         while (k < M1)
               T1 = 2*r(k+1);
               r(k+1) = T1;
               r(k+M1+1) = T1+1;
               k = k+1;
         endwhile
         M1 = 2*M1;
    endwhile
    % end of initializing r()
    % r(n1) is the bitreversal of n1.
    % For n1=1,2,...,N-2 do: if r(n1) > n1, swap x(n1),x(r(n1)).
    n1=1;
    while (n1 < (N-1))
        if ((r(n1+1))>n1)
           temp=x(n1+1);
           x(n1+1)=x(r(n1+1)+1);
           x(r(n1+1)+1)=temp;
        endif
        n1=n1+1;
    endwhile
endfunction


function x=fft1(x,N)
    tss=N/2;        % initial twiddle step size
    p=log2(N);
    Bp=N/2;
    Np=2;
    for n=[1:1:(N/2)]
      T(n)=exp(2*pi*i*(n-1)/N);   % size N/2 twiddle table
    endfor
    for P=[0:1:(p-1)]        % pass loop
      Npp=Np/2;              % # butterflies per sub-block
      BaseT=0;
      for b=[0:1:(Bp-1)]     % sub=block loop
         BaseB=BaseT+Npp;
         for n=[1:1:Npp]     % butterfly loop
            top=x(BaseT+n);
            bot=x(BaseB+n)*T(((n-1)*tss)+1);
            x(BaseT+n)=top+bot;
            x(BaseB+n)=top-bot;
         endfor
         BaseT=BaseT+Np;
      endfor
      Bp=Bp/2;       % # subblocks halved
      Np=Np*2;       % # butterflies/sub-block doubles
      tss=tss/2;     % twiddle step size halved
    endfor
endfunction


function plotx(t,x,y,N)
    subplot(2,1,1);
    plot(t,x);
    hold
    stem(t,x,"r");
    title("Test of DIT FFT with 6-hz sine input");
    xlabel("t (seconds)");
    ylabel("amplitude");
    subplot(2,1,2);
    stem(t(1:32)*N,y,"r");
    text(15,0.45,"|X(f)| output");
    ylabel ("|X(f)|");
    xlabel ("frequency (hz)");
    axis([0,31,0.0,0.55]);
    print -dpng fft_dec6.png;
endfunction


N=64;
t=[0:1/N:(N-1)/N];
f=6;
fs=18;
x1=sin(2*pi*f*t);
x=bitreversal(x1,N);
y=fft1(x,N);
y1=abs(y)/N;
y2=y1(1:(N/2));
plotx(t,x1,y2,N);