c=rand(1,512);% 定制512点抽样%
%%%%求反序,实现反序输入,正序输出%%%%
c(1)=0;c(2)=256;
for i=1:8 % log(512)-1=8%
c((2^i+1):2^(i+1))=c(1:2^i)+2^(8-i);
end
d=rand(1,512);
t=linspace(1,100,512)/100;b=[t;sin(2*pi*t*5)]; %输入5HZ正弦波,
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 取均分512点的正弦值作为fft的初始值
for i=1:512
d(i)=b(2,(c(i)+1));
end%完成初始值的反序%
%%%%求出所有的蝶算指数因子%%%%
x=rand(512,10);a=rand(1,256);q=1;w=rand(1,256);
for o=1:256
w(o)=exp(-j*2*pi*(o-1)/512);
end
%%%%fft变换%%%%
x(:,1)=d';
for k=2:10
i=1;p=1;
for m=1:2^(k-2)%求出每级所需的蝶算指数因子%
q=(m-1)*512/(2^(k-1));
a(p)=w(q+1);p=p+1;
end
for n=1:(512/(2^(k-1)))%对偶节点个数的确定及变换计算%
p=1;
for h=i:(i+(2^(k-2)-1))
x(h,k)=x(h,k-1)+a(p)*x(h+2^(k-2),k-1);
x(h+2^(k-2),k)=x(h,k-1)-a(p)*x(h+2^(k-2),k-1);
p=p+1;
end
i=(2^(k-1))*n+1;
end
end
%%%%取出fft变换后的点数值,取1到50点的值进行作图%%%%
S=abs(x(:,10));
plot((0:49)/50*50,abs(S(1:50)),'linewidth',2)
我是根据一个DITFFT(时域抽选)程序进行修改的DIFFFT(频域抽选),但是画图的时候发现信号和原先他的程序差异很大,我怀疑是我的程序编写的有问题,请大家帮忙改改!
% DIF-FFT
N=256;
x1=[1 1 1 1];
x=[x1 zeros(1,N-length(x1))];
L=1;
M=log2(N);
while L<=M
LE=2^L;
LE1=LE/2;
U=1;
W=exp(-j*2*pi/LE);
J=0;
while J<=LE1-1;
I=J;
while I<=N-1
IP=I+LE1;
T=x(I+1)+x(IP+1);
x(IP+1)=(x(I+1)-x(IP+1))*U;
x(I+1)=T;
I=I+LE;
end
U=U*W;
J=J+1;
end
L=L+1;
end
x;
NV2=N/2;
NM1=N-1;
I=0;
J=0;
while I<NM1
if I<J
T=x(J+1);
x(J+1)=x(I+1);
x(I+1)=T;
end
K=NV2;
while K<=J
J=J-K;
K=K/2;
end
J=J+K;
I=I+1;
end
x;
n=0:N-1;
q=n*2*pi/N;
stem(q,angle(x),'.')
---------------------------------------------
这是DIT-FFT的程序,我就是根据这个来修改的,测试通过,而DIF-FFT则波形不对。
N=256;
x1=[1 1 1 1];
x=[x1 zeros(1,N-length(x1))];
NV2=N/2;
NM1=N-1;
I=0;
J=0;
while I<NM1
if I<J
T=x(J+1);
x(J+1)=x(I+1);
x(I+1)=T;
end
K=NV2;
while K<=J
J=J-K;
K=K/2;
end
J=J+K;
I=I+1;
end
x;
L=1;
M=log2(N);
while L<=M
LE=2^L;
LE1=LE/2;
U=1;
W=exp(-j*2*pi/LE);
J=0;
while J<=LE1-1;
I=J;
while I<=N-1
IP=I+LE1;
T=x(IP+1)*U;
x(IP+1)=x(I+1)-T;
x(I+1)=x(I+1)+T;
I=I+LE;
end
U=U*W;
J=J+1;
end
L=L+1;
end
x;
n=0:N-1;
q=n*2*pi/N;
stem(q,abs(x),'.')
作者:fang5566 发表时间:2007-4-29 15:09:00
第1楼
哪位兄弟帮我改改啊,急需啊!!!!!!!!
作者:zhoujb 发表时间:2007-5-2 13:45:00
第2楼
你好,不知道你还在不在?
作者:fang5566 发表时间:2007-5-4 10:53:00
第3楼
多谢楼上关心 我已经自行解决了
就是LE=2^L;改成LE=2^(M-L+1)
基于Matlab的时间抽取基2FFT算法
function y=myditfft(x)
%本程序对输入序列实现DIT-FFT基2算法,点数取大于等于长度的2的幂次
%------------------------------------
% Leo's fft program(改编网上的一个程序)
%------------------------------------
m=log2(2^nextpow2(length(x))); %求的x长度对应的2的最低幂次m
N=2^m;
if length(x)<N
x=[x,zeros(1,N-length(x))]; %若长度不是2的幂,补0到2的整数幂
end
x;
%--------------------------------------------------------------------------
%对输入序列进行倒序
%如果输入序列的自然顺序号I用二进制数(例如n2n1n0)表示
%则其倒位序J对应的二进制数就是(n0n1n2),这样,在原来自然顺序时应该放x(I)的
%单元,现在倒位序后应放x(J)。
%--------------------------------------------------------------------------
%以下程序相当于以下程序:
%nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; %求1:2^m数列的倒序
%y=x(nxd); %将倒序排列作为初始值
%--------------------------------------------------------------------------
NV2=N/2;
NM1=N-1;
I=0;
J=0;
while I<NM1
if I<J
T=x(J+1);
x(J+1)=x(I+1);
x(I+1)=T;
end
K=NV2;
while K<=J
J=J-K;
K=K/2;
end
J=J+K;
I=I+1;
end
x;
%--------------------------------------------------------------------------
%以下程序解释:
%第一级从x(0)开始,跨接一阶蝶形,再取每条对称
%第二级从x(0)开始,跨接两阶蝶形,再取每条对称
%第m级从x(0)开始,跨接2^(m-1)阶蝶形,再取每条对称....
%--------------------------------------------------------------------------
for mm=1:m %将DFT做m次基2分解,从左到右,对每次分解作DFT运算
Nmr=2^mm;
u=1; %旋转因子u初始化
WN=exp(-j*2*pi/Nmr); %本次分解的基本DFT因子WN=exp(-i*2*pi/Nmr)
for n=1:Nmr/2 %本次跨越间隔内的各次碟形运算
for k=n:Nmr:N %本次碟形运算的跨越间隔为Nmr=2^mm
kp=k+Nmr/2; %确定碟形运算的对应单元下标(对称性)
t=x(kp)*u; %碟形运算的乘积项
x(kp)=x(k)-t; %碟形运算的加法项
x(k)=x(k)+t;
end
u=u*WN; %修改旋转因子,多乘一个基本DFT因子WN
end
end
y=x; %输出