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

1、基于二维图像的 FFT 算法实现1 摘要FFT 算法基本分为两大类：时域抽取法 FFT(Decimation-In-Time FFT，简称 DIT-FFT)和频域抽取法 FFT(Decimation-In-Frequency FFT，简称 DIF-FFT)。本文选取时域抽取法，即 DIT-FFT 算法，利用 matlab 编程实现基于二维图像的 FFT 算法，并选取二维图片，对该图片进行加噪和滤波处理，最后使用逆傅里叶变换恢复原始图片，从而检验该算法的有效性。2 算法描述设序列 的长度为 ， 且满足()xnN， 为自然数2M按 的奇偶把 分解为两个 点的子序列()x/， 1()r0,1.2Nr

2、， 21x则 的 DFT 为()xn()()()knknNNnXkxWx偶 数 奇 数/21/212200()()Nkr krr r由于 222 2jkrNjkrkr krNWeW所以 ，()1()XX0,1.N其中 和 分别是 的 点 DFT，即1()k22xr和 （ /()()DFTX由于 和 均以 为周期，且 ，所以 又可以表示为1()Xk2()/N2NkkW()X(4.1)2(),01,.kX(4.2)()1(22Nk这样，就将 点 DFT 分解为两个 点的 DFT 和 4.1 式以及 4.2 式的运算。4.1 式/和 4.2 式的运算可用图 1 所示的流图符号表示，根据其形状称其为蝶

3、形运算。图 1由于 ， 仍然是偶数，故可以对 点 DFT 再做进一步分解。由 DIT-2MN/ /2NFFT 算法的分解过程可知， 时，其运算流图应有 级蝶形，每一级都由 个2M /2N蝶形运算构成。通过一次分解，可以使运算量减少近一半，从而，DIT-FFT 算法比直接计算 DFT 的运算次数大大减少，可以使运算效率极大提高。而逆傅里叶只要将 DITA-DFT 运算式中的系数 改变为 ，最后乘以 ，就是knNWknN1/DIT-IDFT 的运算公式。AB CA+BCA-BC3 程序流程送入灰度图矩阵 x开始，12MN构造全零矩阵 ，将 送入mx1,L12LB0,J2MLP1,LkJN(:,)W

4、XmkbN输出结束倒序， tep(:,)(:,)Xk0temp(,:)WXkbN,:,:()mYN图 2 DIT-FFT 运算和程序框图/21MNj2,iij()tempxijtKM1j(/2)jKroundjKYY N图 3 倒序程序框图4 实验结果实验中，我们分三组对图片进行不同情况的处理，从而检验该算法的有效性。4.1 未加噪且未滤波如图 4 为原始的二维图片，图 5 为经过快速傅里叶变换后的频谱图，通过与 matlab自带的基二快速傅里叶函数 fft2 变换后的频谱图比较（图 6） ，可以看出 fft2_m（自己编写的快速傅里叶变换函数）函数已经可以很好的实现基二的快速傅里叶变换。且且

5、且且且图 4 原始灰度图且且且且且且且且图 5 原灰度图的频谱图且 且 且 且 且 且 且 且 (且 fft2-m且 且 ) 且 且 且 且 且 且 且 且 (且 fft2且 且 )图 6 fft2_m 实现的原图的频谱图与 fft2 实现的原图的频谱图如图 7 是由 ifft2_m（自己编写的快速逆傅里叶变换函数）函数恢复的原始灰度图，由于我们计算时以二为基，对图片矩阵进行了适当的扩大处理，所以在恢复图的边缘加上一些黑色边缘，通过与原始灰度图对比可知，该算法很好的恢复了原始二维图片，证明该算法是有效的。且且且且且且且且且且且图 7 ifft2_m 的恢复图4.2 未加噪但低通滤波如图 8 是

6、对使用 fft2_m 变换后的频谱进行低通滤波后恢复的灰度图，由该图可以看出，图片变得模糊柔和，即低通滤波器有效的实现了低通滤波，证明该低通滤波器是有效的。且且且且且且且且且且且且图 8 低通滤波后的恢复图4.3 加噪且低通滤波我们对原二维图片加入高斯噪声（图 9） ，并使用巴特沃斯低通滤波器进行滤波（图 11为加入噪声并低通滤波后的频谱图） ，这样再进行逆傅里叶变换以返回原图时可以有更加明显的对比。图 10 为图片加入高斯噪声后经过快速傅里叶变换后的频谱图。且且且且且且且且且且且图 9 加入高斯噪声后的灰度图且且且且且且且且且且且图 10 加入高斯噪声后的频谱图且且且且且且且且且且且且且图

7、11 加入噪声并低通滤波后的频谱图图 12 为加噪并经过低通滤波后，使用逆傅里叶变换（iff2_m）恢复的原始二维图片，由于我们计算时以二为基，对图片矩阵进行了适当的扩大处理，所以在恢复图的边缘加上一些黑色边缘，而且从图 12 中可以看出，图片变得模糊柔和，即低通滤波器有效的实现了低通滤波，而且与图 8 相比该图片多了一些模糊的斑点，这是因为加入了高斯噪声。且且且且且且且且且且且且且图 10 加入噪声且低通滤波后恢复的灰度图综合以上三组实验对图片的分析与比较，可以证明算法函数 ff22_m 和 ifft2_m 以及低通滤波器都是正确的，它们可以有效的实现快速傅里叶变换、快速逆傅里叶变换以及低通

8、滤波，即完成了本次课程设计的要求。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 = 2M2;M1 = nextpow2(X);

9、%当 X=1 时，M1=0N1 = 2M1;m = zeros(N1,N2);%构造一个全零矩阵，用来存放 fft2_m 的结果%将矩阵 x 扩展为 N1*N2 的矩阵for i = 1 : Xm(i,:) = x(i,:), zeros(1,N2 - Y);%若该矩阵的行数小于 N2 或列数小 N1，对其尾部补零end %由于 m 是全零阵，所以只需对行处理即可%对该矩阵的每一行进行倒序for i = 1 : N1m(i,:) = invert_sequence(m(i,:);endtemp = 0;%设置变量，控制行列的变换W = exp(-j*2*pi/N2);%旋转因子%对逆序后矩阵的

10、每一行进行蝶形运算m = fft2_DieXing(m, M2, N2, W, temp);%对该矩阵的每一列进行倒序for i = 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

11、 = invert_sequence(x)N = length(x);%计算输入序列的长度，此处和 fft2_m 一起使用，此值为偶%n = nextpow(x);%N = 2n;M = N/2;%M 是 M 位二进制数最高位的权值j = M + 1;%由于 matlab 的数组下标从 1 开始，所以后移一个下标，%即 下标 j 处的值为 M 位二进制的权值，x(j) = M;N1 = N - 2;%进行倒序时第一个数和最后一个数的位置一定不动for i = 2 : N1 + 1 %仅将第二个到倒数第二个数进行倒序if i = K)%由于 matlab 的数组下标从 1 开始。如果下标超过该权

12、值j = j - K;%首先将该位由 1 变为 0K = round(K/2);%次高位的权值endj = j + K;end%fft2_DieXing.m%对矩阵 m 进行蝶形运算%M 为蝶形运算的最大级数%W 为旋转因子%temp 为变量，用来控制行、列的变换function m = fft2_DieXing(m, M, N, W, temp)for L = 1 : M %M 为蝶形运算的最大级数B = 2(L - 1);for J = 0 : B - 1 %控制旋转因子的系数p = J*2(M - L);WN = Wp;for k = J + 1: 2L : Nkb = k + B;if

13、 temp = 0%进行行计算WX = m(:,kb)*WN;m(:,kb) = m(:,k) - WX;%先减后加，可以不用构造新的矩阵来存储%由于该句修改了 m(:,kb) 的值，而下句需要该原始值%所以要提前计算 WX = m(:,kb)*WNm(:,k) = m(:,k) + WX;else if temp = 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 变换过的矩阵进行逆傅里叶

14、function m = ifft2_m(m)N1,N2 = size(m);%求出矩阵的行列值M1 = nextpow2(m(:,1);%计算列的蝶形计算的最大级数M2 = nextpow2(m(1,:);%计算行的蝶形计算的最大级数%对矩阵的每一列进行逆傅里叶%对该矩阵的每一列进行倒序排列for i = 1 : N2m(:,i) = invert_sequence(m(:,i);endtemp = 1;%对该矩阵的每一列进行傅里叶变换W = exp(j*2*pi/N1);%计算逆变换时的旋转因子m = fft2_DieXing(m, M1, N1, W, temp);for i = 1 :

15、 N2m(:,i) = m(:,i) / N2;%对结果乘以 1 / N2end%对矩阵的每一行进行逆傅里叶%对该矩阵的每一行进行倒序排列for i = 1 : N1m(i,:) = invert_sequence(m(i,:);endtemp = 0;%对该矩阵的每一行进行快速傅里叶变换W = exp(j*2*pi/N2);m = fft2_DieXing(m, M2, N2, W, temp);for i = 1 : N1m(i,:) = m(i,:) / N1;end%butter.m%巴特沃斯滤波器function f = butter(f)%二阶巴特沃斯低通滤波器N1,N2 = size(f);n = 2;d0 = 80;n1 = fix(N1/2);n2 = fix(N2/2);for i = 1 : N1for j = 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