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

基于二维图像的FFT算法实现1 摘要FFT算法基本分为两大类：时域抽取法FFT(Decimation-In-Time FFT，简称DIT-FFT)和频域抽取法FFT(Decimation-In-Frequency FFT，简称DIF-FFT)。本文选取时域抽取法，即DIT-FFT算法，利用matlab编程实现基于二维图像的FFT算法，并选取二维图片，对该图片进行加噪和滤波处理，最后使用逆傅里叶变换恢复原始图片，从而检验该算法的有效性。 2 算法描述设序列的长度为， 且满足，为自然数按的奇偶把分解为两个点的子序列， ， 则的DFT为 由于所以 ，其中和分别是的点DFT，即由于和均以为周期，且，所以又可以表示为 (4.1) (4.2)这样，就将点DFT分解为两个点的DFT和4.1式以及4.2式的运算。4.1式和4.2式的运算可用图1所示的流图符号表示，根据其形状称其为蝶形运算。A+BCAA-BCBC图1由于，仍然是偶数，故可以对点DFT再做进一步分解。由DIT-FFT算法的分解过程可知，时，其运算流图应有级蝶形，每一级都由个蝶形运算构成。通过一次分解，可以使运算量减少近一半，从而，DIT-FFT算法比直接计算DFT的运算次数大大减少，可以使运算效率极大提高。 而逆傅里叶只要将DITA-DFT运算式中的系数改变为，最后乘以，就是DIT-IDFT的运算公式。3 程序流程开始送入灰度图矩阵，输出结束NY构造全零矩阵，将送入倒序，图2 DIT-FFT运算和程序框图NYY图3 倒序程序框图4 实验结果实验中，我们分三组对图片进行不同情况的处理，从而检验该算法的有效性。4.1 未加噪且未滤波如图4 为原始的二维图片，图5 为经过快速傅里叶变换后的频谱图，通过与matlab自带的基二快速傅里叶函数fft2变换后的频谱图比较（图6），可以看出fft2_m（自己编写的快速傅里叶变换函数）函数已经可以很好的实现基二的快速傅里叶变换。图4 原始灰度图图5 原灰度图的频谱图图6 fft2_m实现的原图的频谱图与fft2实现的原图的频谱图如图7是由ifft2_m（自己编写的快速逆傅里叶变换函数）函数恢复的原始灰度图，由于我们计算时以二为基，对图片矩阵进行了适当的扩大处理，所以在恢复图的边缘加上一些黑色边缘，通过与原始灰度图对比可知，该算法很好的恢复了原始二维图片，证明该算法是有效的。图7 ifft2_m的恢复图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 = 2M2;M1 = nextpow2(X);%当X=1 时，M1=0N1 = 2M1;m = zeros(N1,N2);%构造一个全零矩阵，用来存放fft2_m 的结果%将矩阵 x 扩展为N1*N2的矩阵for i = 1 : X m(i,:) = x(i,:), zeros(1,N2 - Y);%若该矩阵的行数小于N2或列数小N1，对其尾部补零end %由于 m 是全零阵，所以只需对行处理即可%对该矩阵的每一行进行倒序for i = 1 : N1 m(i,:) = invert_sequence(m(i,:);endtemp = 0;%设置变量，控制行列的变换W = exp(-j*2*pi/N2);%旋转因子%对逆序后矩阵的每一行进行蝶形运算m = fft2_DieXing(m, M2, N2, W, temp);%对该矩阵的每一列进行倒序for i = 1 : N2 m(:,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 = 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 开始。如果下标超过该权值 j = j - K;%首先将该位由 1 变为 0 K = round(K/2);%次高位的权值 end j = 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 : N kb = k + B; if temp = 0%进行行计算 WX = m(:,kb)*WN; m(:,kb) = m(:,k) - WX;%先减后加，可以不用构造新的矩阵来存储 %由于该句修改了 m(:,kb) 的值，而下句需要该原始值 %所以要提前计算 WX = m(:,kb)*WN m(:,k) = m(:,k) + WX; else if temp = 1%进行列计算 WX = m(kb,:)*WN; m(kb,:) = m(k,:) - WX; m(k,:) = m(k,:) + WX; end end end endend%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,:);%计算行的蝶形计算的最大级数%对矩阵的每一列进行逆傅里叶%对该矩阵的每一列进行倒序排列for i = 1 : N2 m(:,i) = invert_sequence(m(:,i);endtemp = 1;%对该矩阵的每一列进行傅里叶变换W = exp(j*2*pi/N1);%计算逆变换时的旋转因子m = fft2_DieXing(m, M1, N1, W, temp);for i = 1 : N2 m(:,i) = m(:,i) / N2;%对结果乘以 1 / N2end%对矩阵的每一行进行逆傅里叶%对该矩阵的每一行进行倒序排列for i = 1 : N1 m(i,:) = invert_sequence(m(i,:);endtemp = 0;%对该矩阵的每一行进行快速傅里叶变换W = exp(j*2*pi/N2);m = fft2_DieXing(m, M2, N2, W, temp);for i = 1 : N1 m(i,:) = m(i,:) / N1;end%butter.m%巴特沃斯滤波器funct