myfft2.m

该程序实现二维的FFT变换

function X=myfft2(x,N)
%二维fft变换
%N必须是2的整数次幂,表示把图像扩展成N*N的矩阵
[row,col]=size(x);
%如果设定的N比图像原始尺寸小，那么设定的尺寸不合理
if(N<row)|(N<col)
    disp('setted range isnot reasonable');
end
%补零
x(row+1:N,:)=0;
x(1:row,col+1:N)=0;
%序号变换，二进制反转
serial_number=1:N;      %序号
serial_number=serial_number-1;
binary_senum=dec2bin(serial_number);     %序号变为二进制
new_binsenum=fliplr(binary_senum);       %序号的二进制反转
new_senum=bin2dec(new_binsenum);
new_senum=new_senum';
new_senum=new_senum+1;
%开始蝶形变换
%变换用到的常数
L=log2(N);       %变换层数
Wn=exp(-j*2*pi/N);        %复指数因子的常数部分
x_new=x(new_senum,new_senum);    %排序后的矩阵，行和列都已经进行了排序

%首先对列进行变换
X_yonly=zeros(N,N,L+1);   %L+1存放各级变换得到的中间结果
X_yonly(:,:,1)=x_new;     %第一级变换的输入，即原始矩阵
n=1:N;     %对各列进行变换
for m=1:L      %各层变换
    for i=1:N/2^m     %每层内的变换，每层包含N/2^m组
        for k1=1:2^(m-1);    %每组包含2^(m-1)个元素    
            k=k1+(i-1)*2^m;    %结果中的元素编号
            r=(k1-1)*2^(L-m);    %复指数因子的指数
            %以下为两部分的迭代公式，分求和和求差两种
            X_yonly(k,n,m+1)=X_yonly(k,n,m)+X_yonly(k+2^(m-1),n,m)*Wn^r;
            X_yonly(k+2^(m-1),n,m+1)=X_yonly(k,n,m)-X_yonly(k+2^(m-1),n,m)*Wn^r;
        end
    end
end

%再对行进行变换
X_xonly=zeros(N,N,L+1);
X_xonly(:,:,1)=X_yonly(:,:,L+1);    %行变换的第一级输入为列变换最后得到的矩阵
n=1:N;     %对各行进行变换
%以下过程同列变换
for m=1:L      
    for i=1:N/2^m     
        for k1=1:2^(m-1);     
            k=k1+(i-1)*2^m;
            r=(k1-1)*2^(L-m);
            %公式同列变换，注意行号，列号顺序
            X_xonly(n,k,m+1)=X_xonly(n,k,m)+X_xonly(n,k+2^(m-1),m)*Wn^r;
            X_xonly(n,k+2^(m-1),m+1)=X_xonly(n,k,m)-X_xonly(n,k+2^(m-1),m)*Wn^r;
        end
    end
end

%最后FFT变换的结果即行变换最后一级得到的矩阵
X=X_xonly(:,:,L+1);