第5章图像傅里叶变换的MATLAB实现
计算机图像处理中，所谓图像变换就是为达到图像处理的某种目的而使用的一种数学技巧。图像函数经过变换后处理起来较变换前更加简单和方便，由于这种变换是对图像函数而言的，所以称为图像变换。现在研究的图像变换基本都是正交变换，通过正交变换可以减少图像数据的相关性，获取图像的整体特点，有利于用较少的数据量来表示原始图像，这对图像的分析、存储以及传输都是非常有意义的。
5.1傅里叶变换的物理意义
从纯粹的数学意义上看，傅里叶变换是将一个图像函数转换为一系列周期函数来处理； 从物理效果看，傅里叶变换是将图像从空间域转换到频域，其逆变换是将图像从频域转换到空间域。换言之，傅里叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅里叶逆变换是将图像的频率分布函数变换为灰度分布函数。实际上，对图像进行二维傅里叶变换得到的频谱图，就是图像梯度的分布图。傅里叶频谱图上看到的明暗不一的亮点，实际上是图像上某一点与邻域点差异的强弱，即梯度的大小，也即该点的频率大小。如果频谱图中暗的点数更多，那么实际图像是比较柔和的； 反之，如果频谱图中亮的点数多，那么实际图像一定是尖锐的，边界分明且边界两边像素差异较大。
5.2傅里叶变换的定义
傅里叶变换分为一维连续傅里叶变换、一维离散傅里叶变换、二维连续傅里叶变换、二维离散傅里叶变换等。
5.2.1一维连续傅里叶变换
假设函数f(x)为实变量，且在(-∞,+∞)内绝对可积，则f(x)的傅里叶变换定义如下：
F(u)=∫+∞-∞f(x)e-2jπuxdx
假设F(u)可积，求f(x)的傅里叶变换定义为：
f(x)=∫+∞-∞F(u)e2jπuxdu
在积分区间内，f(x)必须满足只含有限个第一类间断点、有限个极值点和绝对可积的条件，并且F(u)也是可积的。正、反傅里叶变换称为傅里叶变换对，并且是可逆的。正、反傅里叶变换的唯一区别是幂的符号。F(u)为一个复函数，由实部和虚部构成，如式(51)所示。
F(u)=R(u)+jI(u)(51)
由于F(u)为复函数，根据复数的特点，可以知道复数的模和实部、虚部之间的关系，如式(52)所示； 复数在实平面上的向量角度和实部、虚部之间的关系，如式(53)所示。
F(u)=［R(u)2+I(u)2］(52)
θ(u)=arctanI(u)R(u)(53)
其中，F(u)称为f(x)的振幅谱或傅里叶谱； F(u)称为f(x)的幅值谱； θ(u)称为f(x)的相位谱； E(u)=F2(u)，E(u)称为f(x)的能量谱。
5.2.2一维离散傅里叶变换
对于有限长序列f(x)(x=0,1,…,N-1)，定义一维离散傅里叶变换对如下：
F(u)=DFT［f(x)］=∑N-1x=0f(x)Wux,u=0,1,…,N-1(54)
f(x)=IDFT［F(u)］=1N∑N-1u=0F(u)W-ux,x=0,1,…,N-1(55)
其中，W=e-j2πN，称为变换核。由式(54)可见，给定序列f(x)，可以求出其傅里叶谱F(u)； 反之由傅里叶谱F(u)也可以求出f(x)。离散傅里叶变换对可以简记为：
f(x)F(u)(56)
离散傅里叶变换的矩阵形式为：
F(0)
F(1)
F(N-1)=W0W0W0…W0
W0W1×1W2×1…W(N-1)×1
W0W1×(N-1)W2×(N-1)…W(N-1)×(N-1)f(0)
f(1)
f(N-1)(57)
f(0)
f(1)
f(N-1)=W0W0W0…W0
W0W-1×1W-2×1…W-(N-1)×1
W0W-1×(N-1)W-2×(N-1)…W-(N-1)×(N-1)F(0)
F(1)
F(N-1)(58)
5.2.3二维连续傅里叶变换
从一维傅里叶变换很容易推广到二维傅里叶变换。
如果假设f(x,y)为实变量，并且E(u,v)可积，则存在以下傅里叶变换对，其中，u、v为频率变量：
F(u,v)=∫+∞-∞∫+∞-∞f(x,y)e-j2π(ux+vy)dxdy(59)
其逆变换为：
f(x,y)=∫+∞-∞∫+∞-∞F(u,v)ej2π(ux+vy)dudv(510)
与一维傅里叶变换一样，二维傅里叶变换可写为如下形式。
振幅谱为：
F(u,v)=［R2(u,v)+I2(u,v)］(511)
相位谱为：
θ(u)=arctanI(u,v)R(u,v)
能量谱为：
p(u,v)=F(u,v)2=［R2(u,v)+I2(u,v)］(512)
振幅谱表明了各正弦分量出现了多少次，而相位谱信息表明了各正弦分量在图像中出现的位置。对于整幅图像来说，只要各正弦分量保持原相位，幅值就不那么重要了。所以大多数实用滤波器都只能影响幅值，而几乎不改变其相位信息。
5.2.4二维离散傅里叶变换
一幅静止的数字图像可以看成二维数据阵列，因此，数字图像处理主要是二维数据处理。一维的DFT和FFT是二维离散信号处理的基础。
将一维离散傅里叶变换推广到二维，则二维离散傅里叶变换对被定义为：
F［f(x,y)］=F(u,v)=1MN∑M-1x=0∑N-1y=0f(x,y)e-j2πuxM+vyN(513)
F-1［F(u,v)］=f(x,y)=∑M-1x=0∑N-1y=0F(u,v)ej2πuxM+vyN(514)
式中，u,x=0,1,2,…,M-1； v,y=0,1,2,…,N-1； x,y为时域变量； u,v为频域变量。
同一维离散傅里叶变换一样，系数1/MN可以在正变换或逆变换中； 也可以在正变换和逆变量前分别乘以1/MN，只要两式系数的乘积等于1/MN即可。
二维离散函数的复数形式、指数形式、振幅、相角、能量谱的表示类似于二维连续函数相应的表达式。
图51矩形函数图
下面通过一个矩形函数来帮助读者加深对二维傅里叶变换的理解。函数f(m,n)只在矩形中心区域有值，取值为1，其他区域取值为0，为了简单起见，将f(m,n)显示为连续形式，矩形函数图如图51所示。
图52显示了其二维离散傅里叶变换后的振幅谱图，其中最大值是F(0,0)，是f(m,n)所有元素的和。从图52中可以看出： 高频部分水平方向的能量比垂直方向的能量更高，这是因为水平方向为窄脉冲，垂直方向为宽脉冲，而窄脉冲比宽脉冲含有更多的高频成分。
另一种显示二维傅里叶变换的方法是将log2F(u,v)作为像素值，使用不同颜色来表示像素值的大小，如图53所示。
图52矩形函数的二维傅里叶变换振幅谱
图53傅里叶变换幅度的对数显示
5.3二维离散傅里叶变换的性质
离散傅里叶变换建立了函数在空间域与频域之间的转换关系，使在空间域难以显示的特征在频域中能够十分清楚地显示出来。在数字图像处理中，经常需要利用这种转换关系和转换规律。下面介绍二维离散傅里叶变换的基本性质。
1. 可分离性
如果图像函数f(x,y)的傅里叶变换为F(u,v)，图像函数g(x,y)的傅里叶变换为G(u,v)，则图像函数h(x,y)=f(x,y)·g(x,y)，它的傅里叶变换H(u,v)=F(u,v)·G(u,v)。
2. 线性
如果图像函数f1(x,y)的傅里叶变换为F1(u,v)，图像函数f2(x,y)的傅里叶变换函数为F2(u,v)，则af1(x,y)+bf2(x,y)的傅里叶变换为aF1(u,v)+bF2(u,v)。
3. 共轭对称性
如果f(x,y)是实函数，则它的傅里叶变换具有共轭对称性：
F(u,v)=F*(-u,-v)
F(u,v)=F(-u,-v)
其中，F*(u,v)是F(u,v)的复共轭。
4. 位移性
如果图像函数f(x,y)的傅里叶变换为F(u,v)，则f(x-x0,y-y0)的傅里叶变换为F(u,v)e-j2π(ux0+vy0)/N，f(x,y)ej2π(u0x+v0y)/N的傅里叶变换为F(u-u0,v-v0)。
5. 尺度变换性
如果图像函数f(x,y)的傅里叶变换为F(u,v)，则图像函数f(ax,by)的傅里叶变换为1|ab|Fua,vb。
6. 周期性
傅里叶变换和逆变换均以N为周期，即：
F(u,v)=F(u+N,v)=F(u,v+N)=F(u+N,v+N)
傅里叶变换的周期性表明，尽管F(u,v)对无穷多个u和v的值重复出现，但只需根据任意周期内的N个值就可以从F(u,v)得到f(x,y)。也就是说，只需一个周期内的变换就可以将F(u,v)完全确定。这一性质对于f(x,y)在空间域中也同样成立。
7. 旋转不变性
如果引入极坐标，使
x=rcosθ
y=rcosθu=ωcosφ
v=ωcosφ
则f(x,y)和F(u,v)分别表示为f(r,θ)，F(ω,φ)。
在极坐标中，存在以下的变换对：
f(r,θ+θ0)F(ω,φ+θ0)
该式表明，如果f(x,y)在空域旋转θ0角度，则相应的傅里叶变换F(u,v)在频域上也旋转同一角度θ0。
8. 卷积性
如果图像函数f(x,y)的傅里叶变换为F(u,v)，图像函数g(x,y)的傅里叶变换为G(u,v)，则图像函数h1(x,y)=f(x,y)*g(x,y)，它的傅里叶变换H1(u,v)=F(u,v)·G(u,v)； 图像函数h2(x,y)=f(x,y)·g(x,y)，它的傅里叶变换H2(u,v)=F(u,v)*G(u,v)。
5.4傅里叶变换的实现
在MATLAB中，通过fft函数进行一维离散傅里叶变换，通过ifft函数进行一维离散傅里叶逆变换。这两个函数用法可通过MATLAB帮助文档了解。MATLAB同时提供了fft2函数进行二维离散傅里叶变换，fft函数与fft2函数的关系为fft2(X)=fft(fft(X).').'。fft2函数与ifft2函数的调用格式为：
Y = fft2(X)： 返回二维离散傅里叶变换，结果Y和X的大小相同。其等价于变换形式fft(fft(X).').'。
Y = fft2(X,m,n)： 在变换前，把X截短或者添加0成m×n的数组，返回结果大小为m×n。
Y = ifft2(X)： 运用快速傅里叶逆变换(IFFT)算法，计算矩阵X的二维离散傅里叶逆变换值Y。Y与X的维数相同。
Y = ifft2(X,m,n)： 计算矩阵X的二维离散傅里叶逆变换矩阵Y。在变换前先将X补零到m×n矩阵。如果m或n比X的维数小，则将X截短。Y的维数为m×n。
y = ifft2(…, 'symmetric')： 强制认为矩阵X为共轭对称矩阵计算矩阵X的二维离散傅里叶逆变换值Y。
y = ifft2(…, 'nonsymmetric')： 不强制认为矩阵X为共轭对称矩阵X的二维离散傅里叶逆变换值Y。
【例51】实现图像的傅里叶变换。
>> clear all;
I=imread('cameraman.tif'); %导入图像
subplot(131);imshow(I);
title('原始图像');
J=fft2(I);%图像傅里叶变换
subplot(132);imshow(J);
title('傅里叶变换后图像');
K=ifft2(J)/255; %傅里叶逆变换
subplot(133);imshow(K);
title('傅里叶逆变换后图像')
运行程序，效果如图54所示。
图54图像的傅里叶变换
在MATLAB中，可以通过fftshift函数将变换后的坐标原点移到频谱图窗口中央，坐标原点是低频，向外是高频。fftshift函数的调用格式为：
Y = fftshift(X)： 把fft函数、fft2函数和fftn函数输出的结果的零频率部分移到数组的中间。对于观察傅里叶变换频谱中间零频率部分十分有效。对于向量，fftshift(X)把X左右部分交换一下； 对于矩阵，fftshift(X)把X的第一、第三象限和第二、第四象限交换； 对于高维数组，fftshift(X)在每维交换X的半空间。
Y = fftshift(X,dim)： 把fftshift操作应用到dim维上。
【例52】图像变亮后进行傅里叶变换。
>> clear all;
I=imread('peppers.png');
J=rgb2gray(I); %将彩色图像转换为灰度图像
J=J*exp(1);%变亮
J(find(J>255))=255;
K=fft2(J);%傅里叶变换
K=fftshift(K);%平移
L=abs(K/256);
figure;
subplot(121);imshow(J);
title('变亮后的图像');
subplot(122);imshow(uint8(L));%频谱图
title('频谱图');
运行程序，效果如图55所示。
图55灰度图像变亮后进行傅里叶变换
5.5傅里叶变换的应用
通过傅里叶变换将图像从时域转换到频域，对其进行相应处理，例如滤波和增强等； 再通过傅里叶变换将图像从频域转换到时域。
5.5.1在图像特征定义中的应用
傅里叶变换可用于与卷积密切相关的运算(correlation)。数字图像处理中的相关运算通常用于匹配模板，可用于对某些模板对应的特征进行定位。
【例53】假如希望在图像text.tif中定位字母a，如图56(a)所示，可以采用下面的方法定位。
解析： 将包含字母a的图像与图像text.png进行相关运算，也就是对字母a的图像和图像text.png进行傅里叶变换，然后利用快速卷积的方法，计算字母a和图像text.png的卷积，提取卷积运算的峰值，即得到在图像text.png中对应字母a的定位结果。
>> clear all;
bw=imread('text.png');
a=bw(32:45,88:98);
subplot(1,2,1),imshow(bw);
title('原始图像');
subplot(1,2,2),imshow(a);
title('模板图像')
运行程序，效果如图56所示。
将模板a和text.png图像进行相关运算，就是先分别对其进行快速傅里叶变换，然后利用快速卷积的方法，计算模板和text.png的卷积。如图57所示，提取卷积运算结果的最大值，即图57右图所示的白色亮点，即得到图像text.png中字母a的定位结果。
>> clear all;
bw=imread('text.png');
a=bw(32:45,88:98);%从图像中提取字母a
C=real(ifft2(fft2(bw).*fft2(rot90(a,2),256,256)));
subplot(121),imshow(C,［］);
title('模板与卷积')
max(C(:))
thresh=60%设定门限
subplot(122),imshow(C>thresh)
title('字母a定位')
图56在图形中定位字母a
运行程序，输出如下，效果如图57所示。
ans =
68
thresh =
60
图57字母a的识别效果
5.5.2在滤波器中的应用
巴特沃斯低通滤波器的传递函数为：
H(u,v)=11+［D(u,v)/D0］2n
其中，D0为截止频率，D(u,v)=u2+v2。由于进行了中心化，频率的中心为M2,N2，因此D(u,v)=u-M22+v-N2212。参数n为巴特沃斯滤波器的阶数，n越大滤波器的形状越陡峭。
巴特沃斯高通滤波器的传递函数为：
H(u,v)=11+［D0/D(u,v)］2n
其参数的意义和巴特沃斯低通滤波器相同。
【例54】对图像进行巴特沃斯低通滤波器。
>> clear all;
I=imread('cameraman.tif');
I=im2double(I);
J=fftshift(fft2(I));%傅里叶变换和平移
［x,y］=meshgrid(-128:127,-128:127); %产生离散数据
z=sqrt(x.^2+y.^2);
D1=10;D2=35; %滤波器的截止
n=6; %滤波器的阶数
H1=1./(1+(z/D1).^(2*n)); %滤波器
H2=1./(1+(z/D2).^(2*n));
K1=J.*H1;
K2=J.*H2;
L1=ifft2(ifftshift(K1));%傅里叶逆变换
L2=ifft2(ifftshift(K2));
subplot(131);imshow(I);
title('原始图像');
subplot(132);imshow(L1); %显示载频频率为10Hz
title('巴特沃斯低通滤波器(10Hz)');
subplot(133);imshow(L2); %载频频率为35Hz
title('巴特沃斯低通滤波器(35Hz)');
运行程序，效果如图58所示。
图58图像的巴特沃斯低通滤波效果
在程序中读入灰度图像，接着对图像进行二维离散傅里叶变换和平移，然后设计巴特沃斯低通滤波器，在频域对图像进行滤波，最后进行二维离散傅里叶逆变换。
【例55】对图像进行巴特沃斯高通滤波器。
>> clear all;
I=imread('cameraman.tif');
I=im2double(I);
J=fftshift(fft2(I));%傅里叶变换和平移
［x,y］=meshgrid(-128:127,-128:127); %产生离散数据
z=sqrt(x.^2+y.^2);
D1=10;D2=35; %滤波器的截止
n1=4; n2=8%滤波器的阶数
H1=1./(1+(D1./z).^(2*n1)); %滤波器
H2=1./(1+(D2./z).^(2*n2));
K1=J.*H1;
K2=J.*H2;
L1=ifft2(ifftshift(K1));%傅里叶逆变换
L2=ifft2(ifftshift(K2));
subplot(131);imshow(I);
title('原始图像');
subplot(132);imshow(L1); %显示载频频率为10Hz
title('巴特沃斯高通滤波器（10Hz）');
subplot(133);imshow(L2); %载频频率为35Hz
title('巴特沃斯高通滤波器（35Hz）');
运行程序，效果如图59所示。
图59图像的巴特沃斯高通滤波效果
在程序中读入灰度图像，接着对图像进行二维离散傅里叶变换和平移，然后设计巴特沃斯高通滤波器，通过频域的相乘来进行滤波，最后进行二维离散傅里叶逆变换。
展开
