fft - Decimation in time
4 views (last 30 days)
Show older comments
Hello
I have this code of a fast fourier transform decimation in time(fft_DIT). I need to change into a fft-decimation in frequency. the difference between DIT and DIF is: -the order of the samples: -On DIT the input is bit-reversed order and the output is natural order; -On DIF the input is natural order and the output is bit-reversed order;
-the butterfy commutation: -On DIT Multiplication is done before additions; -On DIF Multiplication is done after additions;
I think these are the differences. My problem is understand this code and adapt it into a fft decimation in frequency..
I appreciate any help, thank you
if true
%
est=4;
N=2^est;
re=1/N*([0:N-1]);
im=zeros(1,N);
xnovo=re+j*im;
nchg=0;
N2=N/2; flag=zeros(1,N);
for i=1:N-2,
if(flag(i)==0)
dv=N2; ic=1; final=0; pos=i;
for k=1:est,
if (pos/dv >= 1)
pos=pos-dv;
final=final+ic;
end
dv=dv/2; ic=ic*2;
end
if (i~=final)
rev(1+nchg)=i;
rev(N2+nchg)=final;
flag(final)=1;
nchg=nchg+1;
end
end
end
rev(1:nchg);
rev(N2:N2+nchg-1);
for i=1:nchg,
tmp=re(1+rev(i));
re(1+rev(i))=re(1+rev(N2+i-1));
re(1+rev(N2+i-1))=tmp;
tmp=im(1+rev(i));
im(1+rev(i))=im(1+rev(N2+i-1));
im(1+rev(N2+i-1))=tmp;
end
grup=N;
butf=1;
phi=2*pi/N;
for i=0:est-1,
grup=grup/2;
for j=0:grup-1,
for k=0:butf-1,
ptr1=2*j*butf+k;
ptr2=ptr1+butf;
arg=phi*grup*k;
C=cos(arg); S=sin(arg);
retmp=C*re(1+ptr2)+S*im(1+ptr2);
imtmp=C*im(1+ptr2)-S*re(1+ptr2);
re(1+ptr2)=re(1+ptr1)-retmp;
im(1+ptr2)=im(1+ptr1)-imtmp;
re(1+ptr1)=re(1+ptr1)+retmp;
im(1+ptr1)=im(1+ptr1)+imtmp;
end
end
butf=butf*2;
end
for i=1:nchg,
tmp=re(1+rev(i));
re(1+rev(i))=re(1+rev(N2+i-1));
re(1+rev(N2+i-1))=tmp;
tmp=im(1+rev(i));
im(1+rev(i))=im(1+rev(N2+i-1));
im(1+rev(N2+i-1))=tmp;
end
y=(re.^2+im.^2)/N^2;
plot([0:N-1],y);
title('Densidade Espectral');
xlabel('Linha Espectral');
ylabel('Quadrado da Magnitude');
end
0 Comments
Answers (1)
Sunil
on 9 Jan 2024
if true
%
est=4;
N=2^est;
re=1/N*([0:N-1]);
im=zeros(1,N);
xnovo=re+j*im;
nchg=0;
N2=N/2; flag=zeros(1,N);
for i=1:N-2,
if(flag(i)==0)
dv=N2; ic=1; final=0; pos=i;
for k=1:est,
if (pos/dv >= 1)
pos=pos-dv;
final=final+ic;
end
dv=dv/2; ic=ic*2;
end
if (i~=final)
rev(1+nchg)=i;
rev(N2+nchg)=final;
flag(final)=1;
nchg=nchg+1;
end
end
end
rev(1:nchg);
rev(N2:N2+nchg-1);
for i=1:nchg,
tmp=re(1+rev(i));
re(1+rev(i))=re(1+rev(N2+i-1));
re(1+rev(N2+i-1))=tmp;
tmp=im(1+rev(i));
im(1+rev(i))=im(1+rev(N2+i-1));
im(1+rev(N2+i-1))=tmp;
end
grup=N;
butf=1;
phi=2*pi/N;
for i=0:est-1,
grup=grup/2;
for j=0:grup-1,
for k=0:butf-1,
ptr1=2*j*butf+k;
ptr2=ptr1+butf;
arg=phi*grup*k;
C=cos(arg); S=sin(arg);
retmp=C*re(1+ptr2)+S*im(1+ptr2);
imtmp=C*im(1+ptr2)-S*re(1+ptr2);
re(1+ptr2)=re(1+ptr1)-retmp;
im(1+ptr2)=im(1+ptr1)-imtmp;
re(1+ptr1)=re(1+ptr1)+retmp;
im(1+ptr1)=im(1+ptr1)+imtmp;
end
end
butf=butf*2;
end
for i=1:nchg,
tmp=re(1+rev(i));
re(1+rev(i))=re(1+rev(N2+i-1));
re(1+rev(N2+i-1))=tmp;
tmp=im(1+rev(i));
im(1+rev(i))=im(1+rev(N2+i-1));
im(1+rev(N2+i-1))=tmp;
end
y=(re.^2+im.^2)/N^2;
plot([0:N-1],y);
title('Densidade Espectral');
xlabel('Linha Espectral');
ylabel('Quadrado da Magnitude');
end
0 Comments
See Also
Categories
Find more on Multirate Signal Processing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!