How to write fast fourier transform function without using fft function ?

function [A]=DFT(x)
N=length(x);
for k=1:N
    A(k)=0;
    for n=1:N
        A(k)=A(k)+x(n).*exp((-1j).*2.*pi.*(n-1).*(k-1)./N);
    end
end

[Geoff Hayes]

Nur - your above code for the discrete Fourier transform seems correct though I would pre-size A as

 A = zeros(N,1);

prior to entering the outer for loop. As for writing a function equivalent to the MATLAB fft then you could try implementing the Radix-2 FFT which is relatively straightforward though is used for block sizes N that are powers of two.

Answer 1

[Jan Afridi]

I tried to write some code but as I am not a good programmer so the code little bit lengthy... I wish that it will help you. And if you found any error or bugs then also help me to fix it ... Radix-2 FFT Matlab code

Answer 2

[Ryan Black]

Function optimizes the DFT or iDFT for any composite-length signal.

The forward transform is triggered by -1 and takes the time-signal as a row vector:

[Y] = branchandnodefft(y,-1)

The inverse transform is triggered by 1 and takes the frequency-signal as a row vector:

[y] = branchandnodefft(Y,1)

%% branch and node FFT 
%% optimizes any composite-length signal
%% define signal
function [Y] = branchandnodefft(y,typef)
N = length(y);                               %signal length
%% define tree structure
bl = factor(N);                              %branches/node
if length(bl) == 1
    disp('No optimized tree structure exists for prime N... ')
    disp('Sorry, for large N, this may take minutes!');
    Y = DFT(y,N,typef);
else
bl = bl(1:end-1);                            %decimate until Ml is prime
nl = zeros(1,length(bl));                    %node per level preallocation
Ml = zeros(1,length(bl));                    %elements ber branch
nl(1) = 1; bl = cat(2,1,bl); Ml(1) = N;      %declare level 0
for l = 1:length(bl)-1                       %for each tree level
    nl(l+1) = bl(l) * nl(l);                 %find nodes per level
    Ml(l+1) = Ml(l) / bl(l+1);               %find elements per branch @l
end
El = bl.*Ml;                                 %elements per node 
L = l;                                       %define active levels
%% reindex to prime structure
INDEX = zeros(length(Ml)-1,N);               %preallocate reindexing levels
INDEX = cat(1,y,INDEX);                      %load in lvl 0
for l = 1:L                                  %for active tree levels
    n=1; k=1;                                %reset n,k for new level
    for ni = 1:nl(l+1)                       %for each node per level
        for bi = 1:bl(l+1)                   %for each branch per node
            INDEX(l+1,n:n+Ml(l+1)-1) = ...      %decimate branch
                INDEX(l,k:bl(l+1):El(l+1)+k-1);
            n = n + Ml(l+1);                 %hop to next solution location
            k = k+1;                         %iterate to adjacent branch
        end
            k=ni*El(l+1)+1;                  %hop to next node in level
    end
end
%% compute FFT
bl = flip(bl); nl = flip(nl); Ml = flip(Ml); %flip tree instructions 
El = flip(El);
BUTTER = zeros(length(Ml),N);                %preallocate FFT upsample tree
for l = 1:L+1                                %for all active levels
    n = 1; k =1;                             %reset n,k for new level
    if l == 1                                %First level is the DFT level
    for ni = 1:nl(l)                         %for each node per level
        for bi = 1:bl(l)                     %for each branch per node
            Fy = DFT(INDEX(end,n:n+Ml(1)-1),Ml(1),typef); %take DFT
            BUTTER(1,n:n+Ml(1)-1) = Fy;      %load to main array
            n = n + Ml(1);                   %hop to next solution location
        end
    end
    else                                     %Upsample through lvls > 1
    for ni = 1:nl(l-1)                       %for each node per level
        if bl(l-1) == 2                      %if bl == 2 use conj twiddle
            CONJ = ...                       %twiddle odd level
                BUTTER(l-1,n+Ml(l-1):n+2*Ml(l-1)-1) .* ...
                1i.^(typef*2*(0:1/Ml(l-1):1-1/Ml(l-1)));
            BUTTER(l,n:n+El(l-1)-1) = ...    %conj with even lvl
                [BUTTER(l-1,n:n+Ml(l-1)-1) + CONJ ,...
                BUTTER(l-1,n:n+Ml(l-1)-1) - CONJ];
        elseif bl(l-1) > 2                   %else use non-conj twiddle
            for bi = 1:bl(l-1)               %for each branch in node
            tempcomp = ...                   %send to non-conj twiddle fct
            NONCONJfct(BUTTER(l-1,k:k+Ml(l-1)-1),bi,Ml(l-1),bl(l-1),typef);
            BUTTER(l,n:n+El(l-1)-1) = tempcomp + ...
                BUTTER(l,n:n+El(l-1)-1);     %superimpose to main array
            k = k + Ml(l-1);                 %iterate computation location
            end
         end
            n = ni*El(l-1)+1;                %hop to next node in level
            k = n;                           %computations as well
    end
    end
end
Y = BUTTER(end,:);                           %extract FD spectrum
end
end
%% function calls
function [Fy] = DFT(y,N,typef)               %DFT main function                                 
Fy=zeros(1,N);                               %preallocate Fy                   
for Fit=1:1:N                                %test ALL combinations
    Fy(Fit)=sum(y.*1i.^(typef*4*(Fit-1)*(0:1/N:1-1/N))); 
end  
end
function [Fyp] = NONCONJfct(Y,bi,Ml,b,typef) %Non-conjugating twiddle fct
N = Ml*b;                                    %upsample size
FY = Y;                                      %hold local FD
for p = 1:b-1                                %periodically extend FD
FY = cat(2,Y,FY);
end
if bi == 1                                   %branch 1 has no twiddle
Fyp = FY;                    
else                                         %twiddle remaining branches
Fyp = FY.*1i.^(4*typef*(bi-1)*(0:1/N:1-1/N));
end
end
%author: Ryan Black

The algorithm decimates a signal to its prime factorization following the branches and nodes of a factor tree. It takes DFTs of the factored time-signals then up-samples with twiddle factors through the branches and nodes to the tree origin.

Answer 3

[Deepthi Ravula]

Clc Clear all N=length(x); for k=1:N A(k)=0; for n=1:N A(k)=A(k)+x(n).*exp(-1j) end end

Answer 4

[mohammed alenazy]

This may help.

%%%%%%%%%%%%%%%%%%%%% Fast Fourier Transform %%%%%%%%%%%%%%%%%%%%%%
% This Script shows how the Fast Fourier Transform algorithm works for
% better understanding.
clear all; close all; 
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
mm3COS=[];      mm3SIN=[];               % First we allocate two vectors that each
                                         % will enroll elements successively as explained below.
sf=1000;   Ts=1/sf;    t=(1:2*sf)*Ts;    % sf; sampling frequency.  t; signal epoch (duration).
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
y=10*sin(2*pi*20*t)+5*sin(2*pi*40*t)+12*sin(2*pi*30*t);    % As an example
                                                           % we have a signal that has 20,30, and 40 frequencies 
                                                           % with different amplitudes (10,12,5 respectively).
N=length(y);                  % N; the length of the signal.        
for k=1:(N/2)                 % We are using the Loop (iteration) to create pairs of sinusoid signals that have 
                              % the same length as the epoch above.Each paired signals have the same frequency and
                              % magnitude,but differ in the phase. The
                              % first iteration creates a pair that have
                              % a whole wave that occupies the epoch
                              % length. Each following iteration adds
                              % another whole wave.                       
    K=(sf/N)*k;               % K is to adjust k if the length is not one second.
    
    COS3=1*cos(2*pi*K*t);     % This is the first pair sinusoid signal created for each iteration.
    y3COS = y.*COS3;          % The first sinusoid signal is multiplied by the signal epoch.
    mm3COS=[mm3COS sum(y3COS)];   % Then, the summation of this multiplication is included in a growing vector. 
      
    SIN3=1*sin(2*pi*K*t);     % This is the second pair sinusoid signal created for each iteration.
    y3SIN = y.*SIN3;          % The second sinusoid signal is multiplied by the signal epoch.         
    mm3SIN=[mm3SIN sum(y3SIN)];   % Then, the summation of this multiplication is included in a growing vector.
end
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
% The iteration will yield two vectors. Each paired elements in these
% vectors will be squared, summed, and then square rooted (pythagorean equation).

mm3COS=(mm3COS).^2;   % The first vector is squared.
mm3SIN=(mm3SIN).^2;   % The first vector is squared.    

mm3COSIN=sqrt(mm3COS+mm3SIN);   % Then, the square root is taken for the summation.

mm3COSINF=2*abs(mm3COSIN)/N;    % This step is done for two reasons: To normalize the summation and to get the 
                                % magnitude for each frequency.

f=sf.*(1:N/2)/N;                % This is done to normalize the signal epoch to the sampling frequency.
figure(3);          
plot(f,mm3COSINF);
title('Frequency spectrum for a given signal');
xlabel('Frequency(Hz)');      ylabel('Amplitude(volt)');     axis([-60 60 0 20]); 
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

[ABHISHEK GHOSH]

For a discrete values how I implement the fft code. I use above code and data is retrieved from a text file.

[mohammed alenazy]

The (y) vactor has the data that you need to process.

Answer 5

[ROBIN SINGH SOM]

% sep2021
% fft operation without using inbuilt function
% only capable to compute fft upto N = 128 but it can be easily expendable
clc
clear
% generating random signal of length 30
x = randi([-5,5],1,30);
% length of the signal
N = length(x);
% number of stage required
M = log2(N);
% adding zeros 
if (rem(M,1) ~= 0)
    re =rem(M,1);
    M=M-re+1;
    Ne = 2^M;
    x = [x,zeros(1,Ne-N)];
else
    Ne = N;
end
% calc. using inbuild function
x_fft =fft(x);
% bitreversing
x = bitrevorder(x);
% instialization of variables used for different stages
temp = zeros(1,Ne);
temp2 = zeros(1,Ne);
temp3 = zeros(1,Ne);
temp4 = zeros(1,Ne);
temp5 = zeros(1,Ne);
temp6 = zeros(1,Ne);
temp7 = zeros(1,Ne);
% code 
for l = 1:M
    if l==1
        for t=0:2:Ne-1
            temp(t+1:t+2) = temp(t+1:t+2) + myfun(x(t+1:t+2),2^l);
        end
        out = temp;
    end
    if l==2
        for k=0:4:Ne-1
            temp2(k+1:k+4) = temp2(k+1:k+4) + myfun(temp(k+1:k+4),2^l);
        end  
        out= temp2;
    end
    if l==3
        for k=0:8:Ne-1
            temp3(k+1:k+8) = temp3(k+1:k+8) + myfun(temp2(k+1:k+8),2^l);
        end  
        out= temp3;
    end
    if l==4
        for k=0:16:Ne-1
            temp4(k+1:k+16) = temp4(k+1:k+16) + myfun(temp3(k+1:k+16),2^l);
        end  
        out= temp4;
    end
    if l==5
        for k=0:32:Ne-1
            temp5(k+1:k+32) = temp5(k+1:k+32) + myfun(temp4(k+1:k+32),2^l);
        end  
        out= temp5;
    end
    if l==6
        for k=0:64:Ne-1
            temp6(k+1:k+64) = temp6(k+1:k+64) + myfun(temp5(k+1:k+64),2^l);
        end  
        out= temp6;
    end
    if l==7
        for k=0:128:Ne-1
            temp7(k+1:k+128) = temp7(k+1:k+128) + myfun(temp6(k+1:k+128),2^l);
        end  
        out= temp7;
    end
end
x % values of x 

x = 1×32

1 -3 -1 -3 3 -2 -3 3 1 -3 -1 3 -4 1 1 0 0 -1 -1 5 -5 -2 5 2 0 1 5 -2 1 -4

x_fft % fft calc. using inbuild function

x_fft = 
   2.0000 + 0.0000i   3.5084 + 0.4818i   9.3638 +24.4138i  -0.5552 -12.2633i  -6.2929 - 5.5355i   3.0460 -11.6709i  -7.9385 +12.5444i  -8.2722 + 1.0552i  -3.0000 + 3.0000i  27.0506 -14.3335i  11.9385 - 0.7694i -14.7462 + 9.8704i  -7.7071 - 1.5355i  -2.8573 - 1.5504i  -5.3638 +15.1001i  24.8259 - 9.7354i -16.0000 + 0.0000i  24.8259 + 9.7354i  -5.3638 -15.1001i  -2.8573 + 1.5504i  -7.7071 + 1.5355i -14.7462 - 9.8704i  11.9385 + 0.7694i  27.0506 +14.3335i  -3.0000 - 3.0000i  -8.2722 - 1.0552i  -7.9385 -12.5444i   3.0460 +11.6709i  -6.2929 + 5.5355i  -0.5552 +12.2633i   9.3638 -24.4138i   3.5084 - 0.4818i

out   % fft calc. using myfunction (DIT)

out = 
   2.0000 + 0.0000i   3.5084 + 0.4818i   9.3638 +24.4138i  -0.5552 -12.2633i  -6.2929 - 5.5355i   3.0460 -11.6709i  -7.9385 +12.5444i  -8.2722 + 1.0552i  -3.0000 + 3.0000i  27.0506 -14.3335i  11.9385 - 0.7694i -14.7462 + 9.8704i  -7.7071 - 1.5355i  -2.8573 - 1.5504i  -5.3638 +15.1001i  24.8259 - 9.7354i -16.0000 - 0.0000i  24.8259 + 9.7354i  -5.3638 -15.1001i  -2.8573 + 1.5504i  -7.7071 + 1.5355i -14.7462 - 9.8704i  11.9385 + 0.7694i  27.0506 +14.3335i  -3.0000 - 3.0000i  -8.2722 - 1.0552i  -7.9385 -12.5444i   3.0460 +11.6709i  -6.2929 + 5.5355i  -0.5552 +12.2633i   9.3638 -24.4138i   3.5084 - 0.4818i

figure()
hold on
stem(abs(out),'filled','LineStyle',"--","Marker","diamond","Color",'b',"LineWidth",1)
title("Without Inbuilt Function","FontSize",15)
hold off

figure()
hold on
stem(abs(x_fft),'filled','LineStyle',"-.",'Marker',"*","color",'r',"LineWidth",1)
title("With Inbuilt Function","FontSize",15)
hold off

figure()
hold on
stem(abs(x_fft),'filled','LineStyle',"-.",'Marker',"*","color",'r',"LineWidth",1)
stem(abs(out),'filled','LineStyle',"none","Marker","diamond","Color",'b',"LineWidth",1)
title("Overlapping","FontSize",15)
hold off

   

function out = myfun(inp,N)
    twi = twiddle(N);
    inp1 = [inp(1:N/2),inp(1:N/2)];
    inp2 = [inp(N/2+1:N),inp(N/2+1:N)];
    out = zeros(1,N);
    for i=1:N
        out(i) = inp1(i) + twi(i)*inp2(i);
    end
end

function out = twiddle(N)
out = zeros(1,N);
for k=1:N
  out(k) = exp(-1j*2*pi*(k-1)/N);
end
end



% author: Robin Singh Som

Answer 6

[ROBIN SINGH SOM]

% 7sep2021
% fft operation without using inbuilt function 
clc
clear
% generating a random signale with 10000 samples
x = randi([-10,10],1,500);
% length of the signal
N = length(x);
M = log2(N);
% adding zeros 
if (rem(M,1) ~=0)
    re = rem(M,1);
    M = M-re+1;
    Ne = 2^M;
    x = [x,zeros(1,Ne-N)];
else
    Ne = N;
end
N = Ne;
% number of stages
M = log2(N);
% calc. using inbuild function
x_fft =fft(x);
% bitreversing
x = bitrevorder(x);
% code 

for l =1:M
    k = 2^l;
    w = exp((-1j*2*pi*(0:k-1))./(k));
         for t=0:k:N-1 
             z(t+1:t+k) = [x(t+1:t+k/2)+x(t+k/2+1:t+k) .* w(1:k/2) , x(t+1:t+k/2)+x(t+k/2+1:t+k).*w(k/2+1:k)] ;            
         end
         x = z;
         z = 0;         
 end
out = x;
figure()
hold on
stem((0:N-1/N),(abs(out)),'filled','LineStyle',"--","Marker","diamond","Color",'b',"LineWidth",1)
title("Without Inbuilt Function","FontSize",15)
hold off

figure()
hold on
stem((0:N-1/N),(abs(x_fft)),'filled','LineStyle',"-.",'Marker',"*","color",'r',"LineWidth",1)
title("With Inbuilt Function","FontSize",15)
hold off

figure()
hold on
stem(abs(x_fft),'filled','LineStyle',"-.",'Marker',"*","color",'r',"LineWidth",1)
stem(abs(out),'filled','LineStyle',"none","Marker","diamond","Color",'b',"LineWidth",1)
title("Overlapping","FontSize",15)
hold off

%robin singh


   

[Ocean Van Rooi]

Hi what FFT algorithm method is used in this solution is it decimation in time?