# Need FFT Code for Matlab (not built in)

조회 수: 407(최근 30일)
John 15 Mar 2013
Commented: Apostolescu Stefan 17 Jan 2020 11:37
Does anyone have FFT code for Matlab? I don't want to use the built-in Matlab function.

#### 댓글 수: 4

표시 이전 댓글 수: 1
John 15 Mar 2013
I want the radix 2 decimation in time code in Matlab to compare.
Walter Roberson 15 Mar 2013
radix 2 decimination is a DFT technique.
John 15 Mar 2013
Okay I am looking for the O(n log n) algorithm for DFT then =)

로그인 to comment.

### 답변 수 (4)

Youssef Khmou 15 Mar 2013
hi John
Yes there are many versions of Discrete Fourier Transform :
% function z=Fast_Fourier_Transform(x,nfft)
%
% N=length(x);
% z=zeros(1,nfft);
% Sum=0;
% for k=1:nfft
% for jj=1:N
% Sum=Sum+x(jj)*exp(-2*pi*j*(jj-1)*(k-1)/nfft);
% end
% z(k)=Sum;
% Sum=0;% Reset
% end
% return
this is a part of the code , try my submission it contains more details :

#### 댓글 수: 4

표시 이전 댓글 수: 1
Youssef Khmou 15 Mar 2013
FFT is finite Fourier transform, its fast when the length of vector on which is evaluated is ~ to 2^N where N is an integer .
% function z=Fast_Fourier_Transform(x)
%
% N=length(x);
% nfft=2^ceil(log2(N));
% z=zeros(1,nfft);
% Sum=0;
% for k=1:nfft
% for jj=1:N
% Sum=Sum+x(jj)*exp(-2*pi*j*(jj-1)*(k-1)/nfft);
% end
% z(k)=Sum;
% Sum=0;% Reset
% end
% return
but keep in mind that its not fast at all for vectors> 1000 points, built-in fft is very fast.
Abdullah Khan 28 Nov 2013
This is also DFT
The complexity of the code is still N x N
because
% nfft=2^ceil(log2(N)) = nfft (on upper DTFT code )
or
% 2^ceil(log2(N))= N
Walter Roberson 6 Sep 2017
"FFT" stands for "Fast Fourier Transform", which is a particular algorithm to make Discrete Fourier Transform (DFT) more efficient.

로그인 to comment.

Tobias 11 Feb 2014
Walter Roberson 님이 편집함. 6 Sep 2017
For anyone searching an educating matlab implementation, here is what I scribbled together just now:
function X = myFFT(x)
%only works if N = 2^k
N = numel(x);
xp = x(1:2:end);
xpp = x(2:2:end);
if N>=8
Xp = myFFT(xp);
Xpp = myFFT(xpp);
X = zeros(N,1);
Wn = exp(-1i*2*pi*((0:N/2-1)')/N);
tmp = Wn .* Xpp;
X = [(Xp + tmp);(Xp -tmp)];
else
switch N
case 2
X = [1 1;1 -1]*x;
case 4
X = [1 0 1 0; 0 1 0 -1i; 1 0 -1 0;0 1 0 1i]*[1 0 1 0;1 0 -1 0;0 1 0 1;0 1 0 -1]*x;
otherwise
error('N not correct.');
end
end
end

#### 댓글 수: 2

Alireza Salehi 6 Aug 2018
hello, Thanks for this code could you please explain how can i use this function? and explain this more for understanding? thanx a lot
Apostolescu Stefan 17 Jan 2020 11:37
Hello, I'm trying to understand what you have done there. Why is X equal to (Xp-tmp)? For Xp+tmp I understood, it's basically DFT, but why Xp-tmp?

로그인 to comment.

Jan Afridi 6 Sep 2017
Jan Afridi 님이 편집함. 6 Nov 2017

#### 댓글 수: 1

Jan Suchánek 6 Dec 2018
Hello i would like to ask you what had to be different if you would like to use it for harmonical function.
I tried to implement your solution but it was working only for random vectors but not for any type of harmonical.
i was using this for the harmonical function.
______________________
t2=0:1:255;
h1=sin(t2);
x=h1;
X11 = myFFT(x);
______________________
it gave me this error:
Error using *
Inner matrix dimensions must agree.
Error in myFFT (line 18)
X = [1 0 1 0; 0 1 0 -1i; 1 0 -1 0;0 1 0 1i]*[1 0 1 0;1
0 -1 0;0 1 0 1;0 1 0 -1]*x;
Error in myFFT (line 7)
Xp = myFFT(xp);
Error in myFFT (line 7)
Xp = myFFT(xp);
Error in myFFT (line 7)
Xp = myFFT(xp);
Error in myFFT (line 7)
Xp = myFFT(xp);
Error in myFFT (line 7)
Xp = myFFT(xp);
Error in myFFT (line 7)
Xp = myFFT(xp);
Error in ZBS_projekt (line 246)
X11 = myFFT(x);

로그인 to comment.

Ryan Black 11 Mar 2019
Ryan Black 님이 편집함. 14 Mar 2019
Function optimizes the DFT or iDFT for prime and composite signals.
The forward transform is triggered by -1 and takes the time-signal as a row vector:
[Y] = fftmodule(y,-1)
The inverse transform auto-normalizes by N, is triggered by 1 and takes the frequency-signal as a row vector:
[y] = fftmodule(Y,1)
%% FFT
%% optimizes any prime or composite signal
%% define signal
function [Y] = fftmodule(y,typef)
N = length(y); %signal length
%% define tree structure
bl = factor(N); %branches/node
if N <= 3 %just do a naive DFT
Y = DFT(y,N,typef);
elseif length(bl) == 1 %if prime, use Rader's
[r,igm,rm,Nz] = ... %find indices
else
Y = treeandnodefft(y,N,typef,bl); %tree FFT struct
end
if typef == 1 %for inverse
Y = Y / N; %normalize by N
end
end
%% function calls
function [Y] = treeandnodefft(y,N,typef,bl)
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
if Ml(1) > 4
[r,igm,rm,Nz] = ... %find indices
end
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
if Ml(1) > 4
Fy = ...
else
Fy = DFT(INDEX(end,n:n+Ml(1)-1),Ml(1),typef); %take DFT
end
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
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
yg = zeros(1,length(y)-1); %preallocate yg
for n = 1:N-1
yg(n) = y(round(r(n))+1); %index to yg
end
bl = factor(Nz); %new factorization
Y1 = treeandnodefft(yg,Nz,typef,bl); %FFT of yg
Y2 = treeandnodefft(igm,Nz,typef,bl); %FFT of igm
C = Y1.*Y2; %FD pt-wise mult
Ys = treeandnodefft(C,Nz,-typef,bl)/Nz + y(1); %iFFT of C + y(1)
Yn = Yreindex(Ys(1:N-1),rm); %reindex Ys(1:N-1)
Y = cat(2,sum(y),Yn); %cat 0-Hz
end
%find indices
for g = 1:N-1 %guess a generator
rtry = zeros(1,N-1); t=1; %preallocate index vector
for q = 1:N-1 %for all indices
t = t*g; %get new argument
rtry(q) = mod(t,N); %find argument mod
t = mod(t,N); %reduce t
end
if length(rtry) - length(unique(rtry)) == 0
r = rtry; %bijective solution reached
break
end
end
rm = [flip(r(1:end-1)),r(end)]; %m inv index perm
Nz = 2^(ceil(log2(2*N-2))); %find zero-pad length
ig = 1i.^(typef*4*rm/N); %reindex sinusoid
igm = ig; %hold for comp
while length(igm) < Nz
igm = cat(2,ig,igm); %extend igm to >Nz
end
igm = igm(1:Nz); %extract first Nz
end
function [Yn] = Yreindex(Ys,rm)
%find reindex scheme
Yn = zeros(1,length(Ys)); %preallocate Rader Solution
r = [rm(2:end),rm(1)]; %create index vector
for n=1:length(Ys)
Yn(round(r(n))) = Ys(n); %index to Yn
end
end
%author: Ryan Black
The algorithm decimates to N's prime factorization following the branches and nodes of a factor tree. At the prime tree level, algorithm either performs a naive DFT or if needed performs a single Rader's Algorithm Decomposition to (M-1), zero-pads to power-of-2, then proceeds to Rader's Convolution routine. Finally it upsamples through the nodes and branches with twiddle factors for the solution.

#### 댓글 수: 0

로그인 to comment.

이 질문에 답변하려면 로그인을(를) 수행하십시오.

Translated by