冈萨雷斯数字图像处理4

第4章频域处理一、背景知识频域滤波，就是对图像做傅里叶变换后进行的处理

频域滤波在图像增强、图像复原、图像数据压缩等过程中都起着重要作用频域滤波包括低通滤波、高通滤波和高频强调滤波一、二维离散傅里叶变换

令f(x,y)表示一幅大小为MXN的图像，其中

x=0,1,2,…,M-1,y=0,1,2,…,N-1其中u=0,1,2,…,M-1和v=0,1,2,…,N-1

频域系统是由F(u,v)所构成的坐标系统，其中u和v是频率变量利用欧拉公式可以手工计算傅里叶变换一、二维离散傅里叶变换离散傅里叶逆变换由下式给出：其中x=0,1,2,…,M-1和y=0,1,2,…,N-1F(u,v)在这里称为傅里叶系数Matlab中F(1,1)=F(0,0)一、二维离散傅里叶变换在原点处的频率值F(0,0)称为直流分量傅里叶变换的频谱定义为傅里叶变换的相位角定义为：极坐标表示方法：功率谱一、二维离散傅里叶变换数字图像的傅里叶变换关于原点共轭对称傅里叶频谱也关于原点对称周期性

傅里叶逆变换得到的图像也是周期性的，一般只取一个周期一、二维离散傅里叶变换

由于傅里叶变换的周期性和对称性，在(M,N)矩形内无法得到一个完整的周期，因此可以在计算

F(u,v)之前先对f(x,y)乘以(-1)x+y来将变换的中心移到(M/2,N/2)

在Matlab中，数组下标从1开始，则变换后的中心为(M/2+1,N/2)可以用下式来计算频率矩形中心[floor(M/2)+1,

floor(N/2)+1]二、Matlab中的二维DFT

DFT在matlab中使用快速傅里叶变换FFT来实现，工具箱中对应的函数为fft2，语法：F

=

fft2(f)结果若需要填充到PXQ大小，则语法为：F=fft2(f,

P,

Q)傅里叶频谱S

=

abs(F)变换之后的原点移动到频谱中心：FC=fftshift(F),逆操作F=ifftshift(FC)二、Matlab中的二维DFT显示频谱:FC

=

fftshift(F)imshow(abs(FC),

[])对数变换可以拓展显示范围

S2=log(1+abs(FC)) imshow(S2,[])二、Matlab中的二维DFT傅里叶逆变换：

f=ifft2(F)图像处理中，逆变换结果一般只取实部:

f=real(ifft2(F))扩展形式：f

=

real(ifft2(F,

P,

Q))三、频域滤波理论基础：卷积定理

式(1)表明空间函数的卷积可以通过频域函数的乘积得到

式(2)表明频域函数的卷积可以通过空间函数的乘积得到三、频域滤波

频域滤波一般采用式(1)，其中H(u,v)称为传递函数，或称为频域滤波器，滤波器的一般形式：低通滤波器高通滤波器三、频域滤波

为避免折叠误差的干扰，在做频域滤波前要对输入的图像和滤波器进行扩充补零的操作

设f(x,y)的大小为AXB，h(x,y)的大小为CXD，则扩充后的函数大小为PXQ，其中：P>=A+C-1Q>=B+D-1

自定义函数paddedsize()用来计算P、Q的最小偶数值，以满足快速傅里叶变换FFT的计算需要三、频域滤波

为避免折叠误差的干扰，在做频域滤波前要对输入的图像和滤波器进行扩充补零的操作

设f(x,y)的大小为AXB，h(x,y)的大小为CXD，则扩充后的函数大小为PXQ，其中：P>=A+C-1Q>=B+D-1

自定义函数paddedsize()用来计算P、Q的最小偶数值，以满足快速傅里叶变换FFT的计算需要三、频域滤波函数paddedsize()function

PQ

=

paddedsize(AB,

CD,

PARAM)if

nargin

==

1PQ

=

2*AB;elseif

nargin

==

2

&

~ischar(CD)PQ

=

AB

+

CD

-

1;PQ

=

2

*

ceil(PQ

/

2);elseif

nargin

==

2m

=

max(AB);三、频域滤波P

=

2^nextpow2(2*m);PQ

=

[P,

P];elseif

nargin

==

3m

=

max([AB

CD]);P

=

2^nextpow2(2*m);PQ

=

[P,

P];elseerror("Wrong

number

ofinputs.")end三、频域滤波不使用填充的滤波解释：三、频域滤波使用填充的频域滤波代码实例：

PQ=paddedsize(size(f));

Fp=fft2(f,PQ(1),PQ(2));Hp

=

lpfilter("gaussian",PQ(1),PQ(2),2*sig)Gp

=

Hp.*Fp;gp

=

real(ifft2(Gp));gpc=gp(1:size(f,1),

1:size(f,2));对应的空间滤波代码h

=

fspecial(‘gaussian’,

15,

7);gs=imfilter(f,h);三、频域滤波使用填充的频域滤波解释：四、频域滤波的基本步骤

使用函数paddedsize获得填充参数

PQ=paddedsize(size(f))

得到使用填充的傅里叶变换

F=fft2(f,PQ(1),PQ(2))生成大小为PQ(1)XPQ(2)的滤波函数H将变换乘以滤波函数

G=H.*F获得G的傅里叶逆变换的实部

g=real(ifft2(G))将左上角的矩形修剪为原始大小：

g=g(1:size(f,1),1:size(f,2))四、频域滤波的基本步骤频域滤波的基本步骤四、频域滤波的基本步骤

频域滤波函数dftfiltfunction

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));五、从空间滤波器获得频域滤波器

函数freqz2用于计算FIR滤波器的频率响应，可以用来获得频域滤波器H

=

freqz2(h,

R,

C)

其中h是一个二维空间滤波器，H是相应的二维域滤波器，R为H的行数，C为H的列数可以用三维透视图的形式输出H： freqz2(h,R,C)h

=

[1

0

-1;

2

0

-2;1

0

-1];Freqz2(h);五、从空间滤波器获得频域滤波器例4.2

空间滤波与频域滤波的比较F

=

fft2(f);

%计算傅里叶变换S=fftshift(log(1+abs(F));%对数频谱移到中心S=gscale(S);%将频谱缩放到[0~255]h=special(‘sobel’);%h是sobel滤波器

PQ=paddedsize(size(f));%计算填充尺寸

H=freqz2(h,PQ(1),PQ(2));%频域滤波器H1=ifftshift(H);%重排滤波器,原点位于左上角gs=imfilter(double(f),h);%空间域滤波gf=dftfilt(f,H);%频域滤波imshow(gs,

[]);

figure

imshow(gf,

[]);五、在频率域直接生成滤波器频域滤波器的网格数组function[U,V]=dftuv(M,N) u=0:(M-1);%u是0到M-1的数组

v=0:(N-1);%v是0到N-1的数组idx=find(u>M/2);%idx是大于M/2的索引u(idx)

=

u(idx)-M;idy=find(v>N/2);%idy是大于N/2的索引v(idy)

=

v(idy)-N;[V,U]=meshgrid(v,u);%[V,U]是生成的网格五、在频率域直接生成滤波器则C=[01R=[0001101];2

其中函数meshgrid用来生成网格数组，语法：[C,

R]

=

meshgrid(c,r);c和r是输入的行向量，C和R是输出的矩阵C和R的维数为length(c)*length(r)其中C的行是c的副本，R的列是r的副本如c=[0,1];r=[0,1,2];12];五、在频率域直接生成滤波器频域滤波器的网格数组[U,V]=dftuv(8,5);D=U.^2+V.^2;D=

01441125524588591013131016

17

20

20

179

10

13

13

104

5

8

8

51

2

5

5

2五、低通频域滤波器理想低通滤波器(ILP)具有传递函数：

n阶巴特沃兹低通滤波器(BLPF)，截止频率为D0高斯低通滤波器(GLPF)的传递函数为：五、低通频域滤波器例4.4低通滤波PQ=paddedsize(size(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*(D0^2)));%H是频域滤波器g=dftfilt(f,H);%频域滤波

figure,imshow(fftshift(H),[]);

figure,imshow(g,[]);五、低通频域滤波器低通频域滤波器生成函数funciton

[H,

D]=lpfilter(type,

M,

N,

D0,n)[U,V]

=

dftuv(M,

N);D

=

sqrt(U.^2+V.^2);switch

typecase‘ideal’H=double(D<=D0);Case“btw”if nargin

==

4n

=

1;end五、低通频域滤波器低通频域滤波器生成函数H=1./(1+D./D0).^(2*n));case

‘gaussian’H=

exp(-D.^2)./(2*(D0^2)));otherwiseerror(‘unkown

filter

type’);end六、绘图函数绘制线框图函数mesh绘制表面图函数surf绘制阴影图函数shading

interp七、基本的高通滤波器

给定低通滤波器，可以获得相应的高通滤波器：例4.6高通滤波器

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

offgrid

off七、基本的高通滤波器例4.7高通滤波PQ

=

paddedsize(size(f));D0

=

0.05*PQ(1);H=hpfilter(‘gaussian’,

PQ(1),

PQ(2),D0);G=dftfilter(f,

H);Figure,

imshow(g,

[]);八、高频强调滤波