%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% by erinch
% of course this is not the best fft function written, but this function
% one of the best open source fft algorithm.
% the aim of this code is to explain how Discrite Fourier Transform (DFT)
% can be found by using Fast Fourier Transform (FFT) for 2^m number of bits
% using radix2
%
% DFT is used to find the frequencies used in a signal, if there are 
% more than one signal then DFT finds every frequencies in the input signal
% in the result its first value is the DC value of the signal.
% other outputs are frequency harmonics
% 
% -------------------------------------------------------------------------
% Ex: sampling frequency(fs)=50Hz
% output's 1st value is DC value of the signal
% output's 2nd value is 50Hz
% output's 3rd value is 100Hz
% output's 4th value is 150Hz
% .           .           .   
% .           .           .
% .           .           .
% .           .           .
% goes on...
% -------------------------------------------------------------------------
% FFT is used to find DFT faster, and is faster if you try to find most
% frequencies. 
% radix2 is an algorithm used in FFT, which takes values as a power of 2
% ex: at first it add 2 values, then add 4 values, then 8 values etc
% -------------------------------------------------------------------------
% 
% function can take two argument as input
% x=my_fft(y,m)
% m is the number of bits, which is used to take 2^m points from the signal
% if you give more than maximum possible bit or if you do not enter it then
% algorithm takes the maximum possible value for m. if you do not sure what
% to enter for m leave it. 
% -nargin is used to understand how many arguments you enter for 
% the function, it gives us an usage without entering some or/and all input
% arguments for the functions-
% y is our input signal, which we want to find DFT. it can be any signal
% if you just want to test the function you can leave it, in this condition
% input signal is taken as y = sin(2*pi*50*t)+sin(2*pi*120*t) as sample
% 
% reversed order must be used either at the beginning for the input 
% X(n) values or at the end for result Y(n) values. it is the first rule in
% FFT. i use at the beginning.
% a built-in function digitrevorder(numbers_in_decimal,base) is does this 
% work for us.
% then input X(n) values are copied in reversed order
% 
% 
% then we start multiplying and summation procedure
% we can say that state holds our addition number
% in each state it sums 2,4,8,16,32,...etc number of results as 2^m 
% summation in each state,which is come from multiplying
% 
% in each multiplying step it multiplies X(n) values with the coefficients
% called W(kn), which are calculated in each state. because number of W(kn)
% changes in each state and has different values. W(kN)=exp(-j*2*pi*k/N)
% 
% and also the algorithm has different number of summation or subtraction
% in different order in each state which are calculated in for loops
% called plus and minus. they are used in while loop because of their order
% and number
% 
% and finally at end of each state it equalize the results as new inputs
% 
% after the states finished the last output gives us DFT result of 
% the signal
% of course this includes complex values and maybe this values are nothing 
% for you. do not worry. all you need to take their absolute values and
% then to plot them.
% Ex:
% x=my_fft(y,m);
% x=abs(x);
% plot(x);
% bar(x);
% plotting with bar will give better result to understand FFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function x=my_fft(y,m)
% y is the signal which you want to take FFT
% m is the number of bits

%% initializing
% you can call the function without entering a signal or number of points
% because there is default values for y and m which are;
% y = sin(2*pi*50*t)+sin(2*pi*120*t);
% m=3;

if nargin==0
    % initializing
    % test for 8 points
    m=3;
    % test signal y = sin(2*pi*50*t)+sin(2*pi*120*t);
    t=1:.1:50;
    y = sin(2*pi*50*t)+sin(2*pi*120*t);
end

% controlling m bit length
% if m given more than the bits in signal than m equalizing to maximum
% possible value
if nargin ==1 | length(y)>2^m
    m=log2(length(y));
    m=floor(m);
end

% creating 2^m points from signal
y=y(1:2^m);
%     -----------------------------------------------------------------


%% reversed order
% making a decimal order from bit reversed order
% m = # of bits
decimal_order=0:(2^m)-1;
reversed_order=digitrevorder(decimal_order, 2);

%% copying the signal points
% copying original signal y to x in reversed bit order
for index=1:2^m
    x(index)=y(reversed_order(index)+1);
end
% reversed_order(index)+1 is because the matrix must be start from index of
% 1 not 0


%% algorithm

for state=1:m       % state is the each multiplying step
    %     in each state we combine 2,4,8,16 ...etc pt. of DFT
    
    %     -----------------------------------------------------------------
    %     creating W for each step
    %     there are (2^(state-1)) different W for each state
    for kN=1:2^(state-1)
        N=2^state;
        k=kN-1;                 
        % in real equation actually k=kn 
        % but in matlab kn can not be 'zero' because of the matrix index 
        W(kN)=exp(-j*2*pi*k/N);
    end     % end of W
    %     -----------------------------------------------------------------
    
    %     -----------------------------------------------------------------
    %     initializing variables
    index=1;        % index counter of the equation number
    diff=2^(state-1);
    kN=1;           % index counter of the W
    %     multiplying coefficients 
    
    while (index <= (2^m))                         % there are 2^m equation
        % i use "while" instead of "for" loop because i decide the number
        % of addition or subtraction by "for loop" inside and if i did not
        % use "while" it cause index exceed or index cancellation in matrix
        
        for plus=1:2^(state-1)                      
            % plus is the number of addition in a sequance
            % in every state the order changes 
            % (like one addition,2 addition,4 addition ...etc.)
            U(index)=x(index)+W(kN)*x(index+diff);
            index=index+1;
            % i have to increase index value becasue i use "while"
            if kN < 2^(state-1)
                kN=kN+1;
            else
                kN=1;
            end     % end of "if statement" of 'kN' index
            % i have to check the value of 'kN' because in every step there
            % are 2^(state-1) number of W and in every equation it changes
            % whether additon or subtraction, so i have to also increase
            % the index of array of 'W'
        end     % end of "for loop" of equations for additon
       
        for minus=1:2^(state-1)
            % minus is the number of subtraction in a sequance
            % in every state the order changes 
            % (like 1 subtraction,2 subtraction,4 subtraction ...etc.)
            U(index)=x(index-diff)-W(kN)*x(index);
            index=index+1;
            % i have to increase indes value becasue i use "while"
            if kN < 2^(state-1)
                kN=kN+1;
            else
                kN=1;
            end     % end of "if statement" of 'kN' index
        end     % end of "for loop" of equations for subtraction
        
    end         % end of "while" --- equation index
    %     -----------------------------------------------------------------
    
    %     passing the multiplier U() to x(), because i use matrix of x
    %     inside the "while" to represent former values
    for index=1:2^m
        x(index)=U(index);
    end         % end of "for loop" -- copying array to pass next state

end             % end of "for loop" of state value
%     ---------------------------------------------------------------------

x;
% x is the DFT of the signal y for the number of 2^m points