数字图像处理学4

1、 离散傅里叶变换（离散傅里叶变换（DFTDFT） 傅里叶变换是建立以时间为自变量的“信号”与以频率为自变量的“频谱函数”之间的某种变换关系。对于有限长序列的离散傅里叶变换的表达式如下：第第4章章 频率域图像增强频率域图像增强4.1 二维离散傅里叶变换二维离散傅里叶变换 离散傅里叶反变换式为： 二维离散傅里叶变换的定义二维离散傅里叶变换的定义 令f(x,y)表示一幅大小为MN的图像，其中x=0,1,2, ,M-1和y=0,1,2, N-1。f 的二维傅立叶变换可表示为F(u,v)，如下式所示：其中u=0,1,2, ,M-1 和 v=0,1,2, ,N-11010)(2),(),(MxNyNvyM

2、uxjeyxfvuF 频域系统是由F(u,v)所张成的坐标系，其中u和v用作频率变量。由u=0,1,2, ,M-1 和v=0,1,2, ,N-1定义的MN矩形区域常称为频率矩形频率矩形。频率矩形的大小与输入图像的大小相同。 在频域原点处变换的值F(0,0)称为傅里叶变换的直流分量。二维离散傅里叶逆变换的表达式如下：1010)(2),(1),(MuNvNvyMuxjevuFMNyxf 傅里叶频谱傅里叶频谱 令R(u,v)和I(u,v)分别表示F(u,v)的实部和虚部，则傅里叶频谱定义为： 2/122),(),(),(vuIvuRvuF),(),(),(22vuIvuRvuP),(/ ),(arc

3、tan),(vuRvuIvu 傅里叶变换的性质傅里叶变换的性质 共轭对称性： 周期性： ),(),(*vuFvuF),(),(),(),(NvMuFNvuFvMuFvuF 傅里叶频谱的平移傅里叶频谱的平移 对于图像f(x,y), f(x,y)的平移不影响傅里叶变换的幅值，不管图像延X轴平移，沿Y轴平移，还是两个方向上都有平移，对应的幅值谱和原始图像是完全一样的。根据此性质，当对图像进行傅里叶变换后，如果原点不在图像的中心，则可以先进行平移，省掉幅值谱的原点不在图像的中心带来的麻烦。函数函数fft2 F=fft2(f) 返回一个大小仍为MN的傅立叶变换。 使用傅立叶变换进行滤波时，需要对输入数据

4、进行零填充。 F=fft2(f,P,Q) 使用该语法， fft2将使用所要求的0的个数对输入图像进行填充，以便结果函数的大小为PQ。4.2 二维离散傅里叶变换的二维离散傅里叶变换的MATLAB计算计算 频谱函数频谱函数abs s=abs(F) 返回傅里叶变换的频谱。 F=fft2(f); S=abs(F); imshow(S, ) 将频谱的原点移动到频率矩形的中心： F=fft2(f); Fc=fftshift(F); S=abs(Fc); imshow(S, ) 改变频谱显示的动态范围： F=fft2(f); Fc=fftshift(F); S=abs(Fc); S2=log(1+S); i

5、mshow(S2, ) 傅里叶逆变换傅里叶逆变换 f=ifft2(F) f1=real(f) 频域滤波与线性空间滤波频域滤波与线性空间滤波 两个空间函数的卷积可以通过计算两个傅立两个空间函数的卷积可以通过计算两个傅立叶变换函数的乘积的逆变换得到。相反地，两叶变换函数的乘积的逆变换得到。相反地，两个空间函数的卷积的傅立叶变换恰好等于两个个空间函数的卷积的傅立叶变换恰好等于两个函数的傅立叶变换的乘积。函数的傅立叶变换的乘积。 f(x,y)*h(x,y) H(u,v) F(u,v) f(x,y)h(x,y) H(u,v)*F(u,v)4.3 频域滤波频域滤波 频域滤波的原理频域滤波的原理 线性空间滤

6、波由图像f(x,y)与滤波掩膜h(x,y)卷积完成的。根据卷积定理，在频域中可以通过让F(u,v)乘以H(u,v)来得到相同的结果，即空间滤波器的傅立叶变换， H(u,v)称为滤波传递函数。 频域滤波的目的是选择一个滤波器的传递函数，以便按照指定的方式修改F(u,v)。 折叠误差干扰折叠误差干扰 使用DFT进行滤波操作，则图像及其变换都将自动地看成是周期性的。若周期关于函数的非零部分的持续时间很靠近，则对周期函数执行卷积运算会导致相邻周期之间的干扰，称为折叠误差的干扰。可以通过使用零填充零填充来填充函数的方法来避免。 零填充零填充 假设函数f(x,y) 和和h(x,y)的大小分别为的大小分别为

7、A x B和 C x D。通过对f 和 h补零，构造两个大小均为 P x Q 的扩充的函数。通过选择 P A+C-1 Q B+D-1以避免折叠误差。 通过使用用户自定义函数paddedsize计算PQ。 自定义paddedsize函数：计算用于FFT的最小填充尺寸 function PQ=paddedsize(AB,CD,PARAM) if nargin=1 PQ=2*AB; elseif nargin=2&ischar(CD) PQ=AB+CD-1; PQ=2*ceil(PQ/2); elseif nargin=2 m=max(AB); P=2nextpow2(2*m); PQ=P

8、P; elseif nargin=3 m=max(AB CD); P=2nextpow2(2*m); PQ=P P; else error(wrong number of inputs.) end 使用填充和不使用填充的滤波效果使用填充和不使用填充的滤波效果 M,N=size(f); F=fft2(f); sig=10; H=lpfilter(gaussian,M,N,sig); G=H.*F; g=real(ifft2(G); imshow(g, ) PQ=paddedsize(size(f); Fp=fft2(f,PQ(1),PQ(2); Hp=lpfilter(gaussian,PQ(1

9、),PQ(2),2*sig); Gp=Hp.*Fp; gp=real(ifft2(Gp); gpc=gp(1:size(f,1),1:size(f,2); imshow(gp, ) 线性空间滤波线性空间滤波 h=fspecial(gaussian,15,7); gs=imfilter(f,h); imshow(gs, )空间滤波结果空间滤波结果DFT滤波的基本步骤滤波的基本步骤 1.使用paddedsize 获得填充函数： PQ= paddedsize(size(f); 2.得到使用填充的傅立叶变换 F=fft2(f,PQ(1),PQ(2); 3.生成大小为PQ(1) x PQ(2)的滤波函数

10、H 4.将变换乘以滤波函数： G=H.*F; 5.获得G傅立叶逆变换的实部： g=real(ifft2(G); 6.将左上部的矩形修剪为原始大小： g=g(1:size(f,1),1:size(f,2); 用于滤波的用于滤波的M函数函数 以下函数可接收输入图像和一个滤波器，可处理所有的滤波细节并输出经滤波和剪切后的图像。 function g=dftfilt(f,H) F=fft2(f,size(H,1),size(H,2); g=real(ifft2(H.*F); g=g(1:size(f,1),1:size(f,2); 如何将空间滤波器转换为等价的频域滤波器如何将空间滤波器转换为等价的频域

11、滤波器 函数freqz2：计算二维频率响应 H=freqz2(h,R,C) h是一个二维空间滤波器，H是相应的二维频域滤波器，R为H的行数，C为 H 的列数。通常取 R=PQ(1),C=PQ(2)。4.4 从空间滤波器获得频域滤波器从空间滤波器获得频域滤波器 空间域滤波与频域滤波的比较空间域滤波与频域滤波的比较 F=fft2(double(f); S=fftshift(log(1+abs(F); S1=gscale(S); imshow(S1, ) 使用函数使用函数fspecial来生成空间滤波器来生成空间滤波器 h=fspecial(sobel) freqz2(h)sobel对应的频域滤波器

12、图像显示频域滤波器 PQ=paddedsize(size(f); H=freqz2(h,PQ(1),PQ(2); H1=ifftshift(H); imshow(abs(H); figure,imshow(abs(H1) 空间滤波 gs=imfilter(double(f),h); imshow(gs) 频域滤波 gf=dftfilt(f,H1); figure,imshow(gf) 通过创建一副二值图像，可以更清楚地看到边缘： figure,imshow(abs(gs)0.2*abs(max(gs(:) figure,imshow(abs(gs)0.2*abs(max(gf(:)通过计算两者

13、的差来显示两种滤波的差别： d=abs(gs-gf); max(d(:) min(d(:) 建立用于实现频域滤波器的网格数组建立用于实现频域滤波器的网格数组 在频域滤波中，需要计算在频率矩形中矩形频率矩形中矩形原点到任意点的距离原点到任意点的距离，下面通过自定义M函数dftuv来实现提供距离计算的网格数组网格数组。 function U,V=dftuvM,N u=0:(M-1); v=0:(N-1); idx=find(uM/2); u(idx)=u(idx)-M;4.5 在频域中直接生成滤波器在频域中直接生成滤波器 idy=find(vN/2); v(idy)=v(idy)-N; V,U=m

14、eshgrid(v,u); 在大小为在大小为85的矩阵中计算矩阵原点到每一的矩阵中计算矩阵原点到每一点的距离平方点的距离平方 U,V=dftuv8,5; D=U.2+V2; fftshift(D) 低通滤波器（低通滤波器（LPF） 理想低通滤波器的传递函数：理想低通滤波器的传递函数： n阶巴特沃兹低通滤波器的传递函数：阶巴特沃兹低通滤波器的传递函数：00),(),(01),(DvuDDvuDvuHnDvuDvuH20/ ),(11),( 高斯低通滤波器的传递函数：高斯低通滤波器的传递函数：2022/ ),(),(DvuDevuH 低通滤波举例低通滤波举例 PQ=paddedsize(size(

15、f); U,V=dftuv(PQ(1), PQ(2); D0=0.05*PQ(2); F=fft2(f, PQ(1), PQ(2); H=exp(-(U.2+V.2)/(2*(D02); g=dftfilt(f,H); figure,imshow(fftshift(H), ) figure,imshow(log(1+abs(fftshift(F), ) figure,imshow(g, ) 自定义自定义M函数生成所有低通滤波器的传递函数函数生成所有低通滤波器的传递函数 function H,D=lpfilter(type,M,N,D0,n) U, V = dftuv(M, N); D = sq

16、rt(U.2 + V.2); switch type case ideal H = double(D H=fftshift(lpfilter(gaussian,500,500,50); mesh(H(1:10:500, 1:10:500) axis(0 50 0 50 0 1) colormap(0 0 0) axis off grid off view(-25,30) view(-25, 0) 绘制表面图：绘制表面图： H=fftshift(lpfilter(gaussian,500,500,50); surf(H(1:10:500, 1:10:500) axis(0 50 0 50 0 1

17、) colormap(gray) axis off grid off shading interp 高通滤波器（高通滤波器（HPF） 高通滤波器的传递函数：高通滤波器的传递函数： 自定义自定义M函数生成所有高通滤波器的传递函函数生成所有高通滤波器的传递函数数 function H=hpfilter(type,M,N,D0,n),(1),(vuHvuHlphp 高通滤波器的可视化高通滤波器的可视化 H=fftshift(hpfilter(ideal,500,500,50); mesh(H(1:10:500, 1:10:500) axis(0 50 0 50 0 1) colormap(0 0 0) axis off grid off figure, imshow(H, ) 高通滤波举例高通滤波举例 PQ=paddedsize(size(f); D0=0.05*PQ(1); H=hpfilter(gaus