基于二维图像的FFT算法实现

基于二维图像的FFT算法实现1摘要FFT算法基本分为两大类：时域抽取法FFT(Decimation-In-TimeFFT，简称DIT-FFT)和频域抽取法FFT(Decimation-In-FrequencyFFT，简称DIF-FFT)。本文选取时域抽取法，即DIT-FFT算法，利用matlab编程实现基于二维图像的FFT算法，并选取二维图片，对该图片进行加噪和滤波处理，最后使用逆傅里叶变换恢复原始图片，从而检验该算法的有效性。2算法描述设序列()xn的长度为N，且满足2MN，M为自然数按n的奇偶把()xn分解为两个/2N点的子序列1()(2)xrxr，0,1,...12Nr2()(21)xrxr，0,1,...12Nr则()xn的DFT为()()()knknNNnnXkxnWxnW偶数奇数/21/2122001()2()NNkrkkrNNNrrxrWWxrW由于222222jkrNjkrkrkrNNNWeeW所以k()1()2()NXkXkWXk，0,1,...,1kN其中1()Xk和2()Xk分别是1()2xrx和（r)的/2N点DFT，即1()[1()]2()[2()]XkDFTxrXkDFTxr由于1()Xk和2()Xk均以/2N为周期，且2NkkNNWW，所以()Xk又可以表示为()1()2(),0,1,...,12kNNXkXkWXkk(4.1)()1()2(),0,1,...,122kNNNXkXkWXkk(4.2)这样，就将N点DFT分解为两个/2N点的DFT和4.1式以及4.2式的运算。4.1式和4.2式的运算可用图1所示的流图符号表示，根据其形状称其为蝶形运算。图1由于2MN，/2N仍然是偶数，故可以对/2N点DFT再做进一步分解。由DIT-FFT算法的分解过程可知，2MN时，其运算流图应有M级蝶形，每一级都由/2N个蝶形运算构成。通过一次分解，可以使运算量减少近一半，从而，DIT-FFT算法比直接计算DFT的运算次数大大减少，可以使运算效率极大提高。而逆傅里叶只要将DITA-DFT运算式中的系数knNW改变为knNW，最后乘以1/N，就是DIT-IDFT的运算公式。ABCA+BCA-BC3程序流程送入灰度图矩阵x开始112MN，222MN构造全零矩阵m，将x送入m1,LM12LB0,1JB2MLPJ1,,2LkJN(:,)WXmkbWN输出结束倒序，temp(:,)(:,)(:,)(:,)mkbmkWXmkmkWX0temp(,:)WXmkbWN(,:)(,:)(,:)(,:)mkbmkWXmkmkWXYN图2DIT-FFT运算和程序框图/2112MNjMNN2,11iNij()()()()tempxixixjxjtempKM1jK(/2)jjKKroundKjjKYYN图3倒序程序框图4实验结果实验中，我们分三组对图片进行不同情况的处理，从而检验该算法的有效性。4.1未加噪且未滤波如图4为原始的二维图片，图5为经过快速傅里叶变换后的频谱图，通过与matlab自带的基二快速傅里叶函数fft2变换后的频谱图比较（图6），可以看出fft2_m（自己编写的快速傅里叶变换函数）函数已经可以很好的实现基二的快速傅里叶变换。图4原始灰度图图5原灰度图的频谱图原始灰度图原灰度图的频谱图图6fft2_m实现的原图的频谱图与fft2实现的原图的频谱图如图7是由ifft2_m（自己编写的快速逆傅里叶变换函数）函数恢复的原始灰度图，由于我们计算时以二为基，对图片矩阵进行了适当的扩大处理，所以在恢复图的边缘加上一些黑色边缘，通过与原始灰度图对比可知，该算法很好的恢复了原始二维图片，证明该算法是有效的。图7ifft2_m的恢复图原灰度图的频谱图(由fft2-m实现)原灰度图的频谱图(由fft2实现)未加噪且未滤波的恢复图4.2未加噪但低通滤波如图8是对使用fft2_m变换后的频谱进行低通滤波后恢复的灰度图，由该图可以看出，图片变得模糊柔和，即低通滤波器有效的实现了低通滤波，证明该低通滤波器是有效的。图8低通滤波后的恢复图4.3加噪且低通滤波我们对原二维图片加入高斯噪声（图9），并使用巴特沃斯低通滤波器进行滤波（图11为加入噪声并低通滤波后的频谱图），这样再进行逆傅里叶变换以返回原图时可以有更加明显的对比。图10为图片加入高斯噪声后经过快速傅里叶变换后的频谱图。图9加入高斯噪声后的灰度图未加噪但低通滤波的恢复图加入高斯噪声后的灰度图图10加入高斯噪声后的频谱图图11加入噪声并低通滤波后的频谱图加入高斯噪声后的频谱图巴特沃斯低通滤波后的频谱图图12为加噪并经过低通滤波后，使用逆傅里叶变换（iff2_m）恢复的原始二维图片，由于我们计算时以二为基，对图片矩阵进行了适当的扩大处理，所以在恢复图的边缘加上一些黑色边缘，而且从图12中可以看出，图片变得模糊柔和，即低通滤波器有效的实现了低通滤波，而且与图8相比该图片多了一些模糊的斑点，这是因为加入了高斯噪声。图10加入噪声且低通滤波后恢复的灰度图综合以上三组实验对图片的分析与比较，可以证明算法函数ff22_m和ifft2_m以及低通滤波器都是正确的，它们可以有效的实现快速傅里叶变换、快速逆傅里叶变换以及低通滤波，即完成了本次课程设计的要求。5参考文献【1】《数字信号处理》（第二版），丁玉美、高西全编著【2】《matlab7.x》，周建兴、岂兴明等著加噪低通滤波后恢复的灰度图6附件%%%%%fft2_m.m%DIT-FFT快速傅里叶变换%输入x为灰度图矩阵function[m]=fft2_m(x)x=double(x);%进行数据类型转换，MATLAB不支持图像的无符号整型的计算X=length(x(:,1));%计算矩阵x的行数Y=length(x(1,:));%计算矩阵x的列数M2=nextpow2(Y);%当Y=1时，M2=0;N2=2^M2;M1=nextpow2(X);%当X=1时，M1=0N1=2^M1;m=zeros(N1,N2);%构造一个全零矩阵，用来存放fft2_m的结果%将矩阵x扩展为N1*N2的矩阵fori=1:Xm(i,:)=[x(i,:),zeros(1,N2-Y)];%若该矩阵的行数小于N2或列数小N1，对其尾部补零end%由于m是全零阵，所以只需对行处理即可%对该矩阵的每一行进行倒序fori=1:N1m(i,:)=invert_sequence(m(i,:));endtemp=0;%设置变量，控制行列的变换W=exp(-j*2*pi/N2);%旋转因子%对逆序后矩阵的每一行进行蝶形运算m=fft2_DieXing(m,M2,N2,W,temp);%对该矩阵的每一列进行倒序fori=1:N2m(:,i)=invert_sequence(m(:,i));endtemp=1;W=exp(-j*2*pi/N1);%旋转因子%对逆序后矩阵的每一列进行蝶形运算m=fft2_DieXing(m,M1,N1,W,temp);%%%%%invert_sequence.m%完成DIT-FFT变换中序列的逆序%该函数在fft2_m函数中被调用%设计参考《数字信号处理》（第二版）丁玉美P104%输入x为一维行矩阵function[x]=invert_sequence(x)N=length(x);%计算输入序列的长度，此处和fft2_m一起使用，此值为偶%n=nextpow(x);%N=2^n;M=N/2;%M是M位二进制数最高位的权值j=M+1;%由于matlab的数组下标从1开始，所以后移一个下标，%即下标j处的值为M位二进制的权值，x(j)=M;N1=N-2;%进行倒序时第一个数和最后一个数的位置一定不动fori=2:N1+1%仅将第二个到倒数第二个数进行倒序ifij%如果该下标小于最高位的权值，直接进行数值交换temp=x(i);x(i)=x(j);x(j)=temp;endK=M;while(j-1=K)%由于matlab的数组下标从1开始。如果下标超过该权值j=j-K;%首先将该位由1变为0K=round(K/2);%次高位的权值endj=j+K;end%%%%%fft2_DieXing.m%对矩阵m进行蝶形运算%M为蝶形运算的最大级数%W为旋转因子%temp为变量，用来控制行、列的变换functionm=fft2_DieXing(m,M,N,W,temp)forL=1:M%M为蝶形运算的最大级数B=2^(L-1);forJ=0:B-1%控制旋转因子的系数p=J*2^(M-L);WN=W^p;fork=J+1:2^L:Nkb=k+B;iftemp==0%进行行计算WX=m(:,kb)*WN;m(:,kb)=m(:,k)-WX;%先减后加，可以不用构造新的矩阵来存储%由于该句修改了m(:,kb)的值，而下句需要该原始值%所以要提前计算WX=m(:,kb)*WNm(:,k)=m(:,k)+WX;elseiftemp==1%进行列计算WX=m(kb,:)*WN;m(kb,:)=m(k,:)-WX;m(k,:)=m(k,:)+WX;endendendendend%%%%%ifft2_m.m%实现DIT-IFFT%配合fft2_m使用%对fft2_m变换过的矩阵进行逆傅里叶function[m]=ifft2_m(m)[N1,N2]=size(m);%求出矩阵的行列值M1=nextpow2(m(:,1));%计算列的蝶形计算的最大级数M2=nextpow2(m(1,:));%计算行的蝶形计算的最大级数%对矩阵的每一列进行逆傅里叶%对该矩阵的每一列进行倒序排列fori=1:N2m(:,i)=invert_sequence(m(:,i));endtemp=1;%对该矩阵的每一列进行傅里叶变换W=exp(j*2*pi/N1);%计算逆变换时的旋转因子m=fft2_DieXing(m,M1,N1,W,temp);fori=1:N2m(:,i)=m(:,i)/N2;%对结果乘以1/N2end%对矩阵的每一行进行逆傅里叶%对该矩阵的每一行进行倒序排列fori=1:N1m(i,:)=invert_sequence(m(i,:));endtemp=0;%对该矩阵的每一行进行快速傅里叶变换W=exp(j*2*pi/N2);m=fft2_DieXing(m,M2,N2,W,temp);fori=1:N1m(i,:)=m(i,:)/N1;end%%%%%butter.m%巴特沃斯滤波器functionf=butter(f)%二阶巴特沃斯低通滤波器[N1,N2]=size(f);n=2;d0=80;n1=fix(N1/2);n2=fix(N2/2);fori=1:N1forj=1:N2d=sqrt((i-n1)^2+(j-n2)^2);h=1/(1+0.414*(d/d0)^(2*n));f(i,j)=h*f(i,j);endend