How to write fast fourier transform function without using fft function ?
89 views (last 30 days)
Show older comments
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
1 Comment
Geoff Hayes
on 18 May 2015
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.
Answers (6)
Jan Afridi
on 6 Sep 2017
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
0 Comments
Ryan Black
on 11 Mar 2019
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.
0 Comments
Deepthi Ravula
on 19 Feb 2020
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
0 Comments
mohammed alenazy
on 24 Jan 2021
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
2 Comments
ABHISHEK GHOSH
on 19 Apr 2021
For a discrete values how I implement the fft code. I use above code and data is retrieved from a text file.
ROBIN SINGH SOM
on 1 Sep 2021
% 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_fft % fft calc. using inbuild function
out % fft calc. using myfunction (DIT)
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
0 Comments
ROBIN SINGH SOM
on 7 Sep 2021
% 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
1 Comment
Ocean Van Rooi
on 3 Nov 2021
Hi what FFT algorithm method is used in this solution is it decimation in time?
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!