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;
x_new=x(new_senum,new_senum);
%开始蝶形变换
%变换用到的常数
L=log2(N); %变换层数
Wn=exp(-j*2*pi/N); %复指数因子的常数部分
%首先对列进行变换
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);
end
