fenji.m

我的FFT分级调试程序

clear;           
f=100;            %待测信号频率
fs=12800;         %采样频率
N=128;           %采样点数
m=log2(N);        %总级数
t=linspace(0,N/fs,N);
count=0:1023;
y=cos(2*pi*count*1000/1024);
yy=fft(y,1024);
plot(abs(yy));
stem(y);
%y=sin(2*pi*f*t);  %输入序列
%q=ones(1,64);
%h=-ones(1,64);
%y=ones(1,128);
%y(1:64)=q;
%y(65:128)=h;
%y=ones(1,128);

p=7;             %第p级运算结果
%***********************方波
% q=ones(1,64);
% h=-ones(1,64);
% y=ones(1,128);
% y(1:64)=q;
% y(65:128)=h;
% y=ones(1,128);
%***********************正弦波
% for count=1:128
%     y(count)=1000*cos(2*pi*(count-1)/128);
% end
l=1;
d=2^(l-1);   %第l级 输入间隔为d  j取值0~d-1 k=j*2^(11-l) exp(-j*2*pi*k/2048)
%**************************锯齿波
% for count=1:128
%     y(count)=count-1;
% end

%for count=1:128
    %y(count)=1000*cos(2*pi*(count-1)/128);
%end
% q=ones(1,64);
% h=-ones(1,64);
% y=ones(1,128);
% y(1:64)=q;
% y(65:128)=h;
%****************************按二进制倒序输入序列
for h=0:(N-1)    
   o=dec2bin(h,m);
   e=fliplr(o);
    g=bin2dec(e);
  x(h+1)=y(g+1);
end
%****************************进行前p级蝶形运算
for l=1:p        
    d=2^(l-1);
    %***********************计算旋转因子
    for j=0:(d-1) 
        k=j*2^(m-l);
        w(j+1)= exp(-i*2*pi*k/N);
    end

    for n=1:(2^l):((N+1)-2^l);
        %*****************蝶形运算基本单元
        for j=0:(d-1)
            a(n)=x(n)+x(n+d)*w(1+j);    
            a(n+d)=x(n)-x(n+d)*w(1+j);  %a存中间结果
            x(n)=a(n);
            x(n+d)=a(n+d);
            n=n+1;                      %进入下一个交叉了的蝶形运算
        end
            %n=n-d;                     %进入下一个不交叉的蝶形运算
    end
end
aaa=abs(x);
rp=real(x);                             %取最后结果实部
ip=imag(x);                             %取最后结果虚部
lp(1:N)=rp; 
lp((N+1):2*N)=ip;

[fu,pos]=max(abs(x));
freq=fs/N*(pos-1);                     %计算到最后一级求出频率
%stem(lp);
figure;
stem(rp);
stem(aaa);
%title('分级运算结果');