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

145 次查看（过去 30 天）
Nur Kaynar
Nur Kaynar 2015-5-13
评论： Ocean Van Rooi 2021-11-3
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 个评论
Geoff Hayes
Geoff Hayes 2015-5-18
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（6 个）

Jan Afridi
Jan Afridi 2017-9-6
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
Ryan Black
Ryan Black 2019-3-11
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.
Deepthi Ravula
Deepthi Ravula 2020-2-19
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
mohammed alenazy
mohammed alenazy 2021-1-24
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
ROBIN SINGH SOM
ROBIN SINGH SOM 2021-9-1
Ran in:
% 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
ROBIN SINGH SOM
ROBIN SINGH SOM 2021-9-7
Ran in:
% 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

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by