数字信号实验3：FFT算法

数字信号实验3：FFT算法

一、DIT-FFT(时域抽取法基2FFT)
1、原理：
2、代码

function [Result] = my_fft(xn,M)
%DIT-FFT(时域抽取法基2FFT)
%正序输入，正序输出！
%调整补零
%(序列长度不够，后面补0,扩展序列长度为2^M,可以更快的实现FFT，提高计算效率)
b=length(xn);               
if(b<2^M)
    xn=[xn,zeros(1,(2^M)-b)];
end
b=length(xn);                %获取补零后xn长度
A=zeros(1,b);                %设置xn转换数组
N=2^M;                       %计算点数
nxd=bin2dec(fliplr(dec2bin([1:N]-1,M)))+1;%倒序排列序号
%dec2bin()：把非负的数字向下取整转成二进制的字符串输出
%fliplr():翻转矩阵
%bin2dec():把字符串形式的二进制数转换成十进制数
xn=xn(nxd);                               %倒序xn         
for i=1:N                                 %N个1点DFT，xn本身，赋值到数组A
   A(i)=xn(i);
end
for L = 1:M                               %DIT-FFT变换，M级蝶形变换
    B = 2^(L-1);                          %两个输入数据距离
    for J = 0:B-1;%J从0-B-1，歩长1         %旋转因子处理
        %B相同执行的是同一级的蝶形运算
        P=2^(M-L)*J;
        %P变化，执行不同的wn^P
        for k=(J+1):2^L:N;%k从J+1-N，歩长2^L        %本次蝶形运算跨越时间
            %P相同，环执行同一个wn^P
            %每个蝶形计算要运算一个复数乘法和两次复数加减法
            W=exp(-j*2*pi*P/N);           %对应旋转因子
            T=A(k)+A(k+B)*W;              %进行蝶形运算
            A(k+B)=A(k)-A(k+B)*W;
            A(k)=T;
            % A(k+B)为X2(k)
            %A(k)为X1(k)
            %T为x1(k)+x2(k)*wn^P
        end
    end
end
Result=A;                                     %输出Result为A
end

3、验证

%绘制出x(n)图形
M=128;     %x(n)的长度
n=0:M-1;
x1=cos((pi*n)/6);
xn=4*cos((pi*n)/6)+3*cos((pi*n)/3)+4*cos((pi*n)/7); 
subplot(2,1,1);stem(n,xn);
xlabel('n');ylabel('x(n)');title('x(n)');grid on;

%对x(n)作128点DFT运算
N1=128; %N1=128个抽样点
k1=0:N1-1;
Xk1=xn*(exp((-j)*((2*pi)/N1))).^(n'*k1); %DFT[x(n)]得到X(k1)
magXk1=abs(Xk1);
angXk1=angle(Xk1);
figure
subplot(2,1,1);stem(k1,magXk1);
xlabel('k1');ylabel('magXk1');title('magXk1');grid on;
subplot(2,1,2);stem(k1,angXk1);
xlabel('k1');ylabel('angXk1');title('angXk1');grid on;

%利用MATLAB内部的fft(x,N)实现序列x(n)作128点DFT运算
Xk2=fft(xn,128);
[a2,b2]=size(Xk2);
o2=1:b2;
magXk2=abs(Xk2);
angXk2=angle(Xk2);
figure
subplot(2,1,1);stem(o2,magXk2);
xlabel('b2');ylabel('magXk2');title('magXk2');grid on;
subplot(2,1,2);stem(o2,angXk2/pi);
xlabel('b2');ylabel('angXk2');title('angXk2');grid on;

%利用my_fft(x,N)实现序列x(n)作2^7=128点DFT运算
Xk3=my_fft(xn,7);
[a3,b3]=size(Xk3);
o3=1:b3;
magXk3=abs(Xk3);
angXk3=angle(Xk3);
figure
subplot(2,1,1);stem(o3,magXk3);
xlabel('b3');ylabel('magXk3');title('magXk3');grid on;
subplot(2,1,2);stem(o3,angXk3);
xlabel('b3');ylabel('angXk3');title('angXk3');grid on;

4、图形展示

所得的图形结果一样！