signal.m

FFT算法

F=input('输入信号的频率:');
N=input('输入采样长度:');
T=input('输入采样间隔时间:');
state=input('是否补零? "1"代表是, "0"代表否:');
if state
    num=input('输入补零的个数:');
end
x=0:T:T*(N-1);
A=sin(2*pi*F*x);
if state
    A=[A,zeros(1,num)];
end %补零
%以上是信号的产生，下面是按时间抽取的FFT算法
N=length(A);
M=log2(N);
AA=ones(1,N); %验证位，供码位倒置时使用
for a=0:N-1
    b=dec2bin(a,M);
    b=fliplr(b);
    b=bin2dec(b);
    if AA(a+1)&&AA(b+1) %已倒置的位标0，未倒置的标0
        temp=A(a+1);A(a+1)=A(b+1);A(b+1)=temp;
        AA(a+1)=0;AA(b+1)=0;
    end
end %码位倒置
od=1:N;mul=1;
for a=1:M
    for b=1:mul
        k(b)=(b-1)*N/(mul*2);
    end %碟形系数中W因子的指数
    mul=mul*2;
    m=1;
    for b=1:2:N-1
        if m>length(k)
            m=1;
        end
        temp=A(od(b));
        A(od(b))=temp+A(od(b+1))*exp(-j*2*pi/N*k(m));
        A(od(b+1))=temp-A(od(b+1))*exp(-j*2*pi/N*k(m));
        m=m+1;
    end %一级FFT运算
    od1=1:N;
    for b=1:2:N-1
        c=od1(1);d=c+mul;
        od(b)=od1(1);od(b+1)=d;
        od1(1)=[];od1(find(od1==d))=[];
    end %FFT运算序列中点的运算顺序
    k=[];
end
%以下是频谱的模特性图
t=0:N-1;
A=abs(A); %取模
A=A./max(A); %归一化
stem(t,A)