岡薩雷斯數(shù)字圖像處理4

第4章頻域處理一、背景知識頻域?yàn)V波，就是對圖像做傅里葉變換后進(jìn)行的處理

頻域?yàn)V波在圖像增強(qiáng)、圖像復(fù)原、圖像數(shù)據(jù)壓縮等過程中都起著重要作用頻域?yàn)V波包括低通濾波、高通濾波和高頻強(qiáng)調(diào)濾波一、二維離散傅里葉變換

令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

頻域系統(tǒng)是由F(u,v)所構(gòu)成的坐標(biāo)系統(tǒng)，其中u和v是頻率變量利用歐拉公式可以手工計(jì)算傅里葉變換一、二維離散傅里葉變換離散傅里葉逆變換由下式給出：其中x=0,1,2,…,M-1和y=0,1,2,…,N-1F(u,v)在這里稱為傅里葉系數(shù)Matlab中F(1,1)=F(0,0)一、二維離散傅里葉變換在原點(diǎn)處的頻率值F(0,0)稱為直流分量傅里葉變換的頻譜定義為傅里葉變換的相位角定義為：極坐標(biāo)表示方法：功率譜一、二維離散傅里葉變換數(shù)字圖像的傅里葉變換關(guān)于原點(diǎn)共軛對稱傅里葉頻譜也關(guān)于原點(diǎn)對稱周期性

傅里葉逆變換得到的圖像也是周期性的，一般只取一個周期一、二維離散傅里葉變換

由于傅里葉變換的周期性和對稱性，在(M,N)矩形內(nèi)無法得到一個完整的周期，因此可以在計(jì)算

F(u,v)之前先對f(x,y)乘以(-1)x+y來將變換的中心移到(M/2,N/2)

在Matlab中，數(shù)組下標(biāo)從1開始，則變換后的中心為(M/2+1,N/2)可以用下式來計(jì)算頻率矩形中心[floor(M/2)+1,

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

DFT在matlab中使用快速傅里葉變換FFT來實(shí)現(xiàn)，工具箱中對應(yīng)的函數(shù)為fft2，語法：F

=

fft2(f)結(jié)果若需要填充到PXQ大小，則語法為：F=fft2(f,

P,

Q)傅里葉頻譜S

=

abs(F)變換之后的原點(diǎn)移動到頻譜中心：FC=fftshift(F),逆操作F=ifftshift(FC)二、Matlab中的二維DFT顯示頻譜:FC

=

fftshift(F)imshow(abs(FC),

[])對數(shù)變換可以拓展顯示范圍

S2=log(1+abs(FC)) imshow(S2,[])二、Matlab中的二維DFT傅里葉逆變換：

f=ifft2(F)圖像處理中，逆變換結(jié)果一般只取實(shí)部:

f=real(ifft2(F))擴(kuò)展形式：f

=

real(ifft2(F,

P,

Q))三、頻域?yàn)V波理論基礎(chǔ)：卷積定理

式(1)表明空間函數(shù)的卷積可以通過頻域函數(shù)的乘積得到

式(2)表明頻域函數(shù)的卷積可以通過空間函數(shù)的乘積得到三、頻域?yàn)V波

頻域?yàn)V波一般采用式(1)，其中H(u,v)稱為傳遞函數(shù)，或稱為頻域?yàn)V波器，濾波器的一般形式：低通濾波器高通濾波器三、頻域?yàn)V波

為避免折疊誤差的干擾，在做頻域?yàn)V波前要對輸入的圖像和濾波器進(jìn)行擴(kuò)充補(bǔ)零的操作

設(shè)f(x,y)的大小為AXB，h(x,y)的大小為CXD，則擴(kuò)充后的函數(shù)大小為PXQ，其中：P>=A+C-1Q>=B+D-1

自定義函數(shù)paddedsize()用來計(jì)算P、Q的最小偶數(shù)值，以滿足快速傅里葉變換FFT的計(jì)算需要三、頻域?yàn)V波

為避免折疊誤差的干擾，在做頻域?yàn)V波前要對輸入的圖像和濾波器進(jìn)行擴(kuò)充補(bǔ)零的操作

設(shè)f(x,y)的大小為AXB，h(x,y)的大小為CXD，則擴(kuò)充后的函數(shù)大小為PXQ，其中：P>=A+C-1Q>=B+D-1

自定義函數(shù)paddedsize()用來計(jì)算P、Q的最小偶數(shù)值，以滿足快速傅里葉變換FFT的計(jì)算需要三、頻域?yàn)V波函數(shù)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);三、頻域?yàn)V波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三、頻域?yàn)V波不使用填充的濾波解釋：三、頻域?yàn)V波使用填充的頻域?yàn)V波代碼實(shí)例：

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));對應(yīng)的空間濾波代碼h

=

fspecial(‘gaussian’,

15,

7);gs=imfilter(f,h);三、頻域?yàn)V波使用填充的頻域?yàn)V波解釋：四、頻域?yàn)V波的基本步驟

使用函數(shù)paddedsize獲得填充參數(shù)

PQ=paddedsize(size(f))

得到使用填充的傅里葉變換

F=fft2(f,PQ(1),PQ(2))生成大小為PQ(1)XPQ(2)的濾波函數(shù)H將變換乘以濾波函數(shù)

G=H.*F獲得G的傅里葉逆變換的實(shí)部

g=real(ifft2(G))將左上角的矩形修剪為原始大?。?/p>

g=g(1:size(f,1),1:size(f,2))四、頻域?yàn)V波的基本步驟頻域?yàn)V波的基本步驟四、頻域?yàn)V波的基本步驟

頻域?yàn)V波函數(shù)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));五、從空間濾波器獲得頻域?yàn)V波器

函數(shù)freqz2用于計(jì)算FIR濾波器的頻率響應(yīng)，可以用來獲得頻域?yàn)V波器H

=

freqz2(h,

R,

C)

其中h是一個二維空間濾波器，H是相應(yīng)的二維域?yàn)V波器，R為H的行數(shù)，C為H的列數(shù)可以用三維透視圖的形式輸出H： freqz2(h,R,C)h

=

[1

0

-1;

2

0

-2;1

0

-1];Freqz2(h);五、從空間濾波器獲得頻域?yàn)V波器例4.2

空間濾波與頻域?yàn)V波的比較F

=

fft2(f);

%計(jì)算傅里葉變換S=fftshift(log(1+abs(F));%對數(shù)頻譜移到中心S=gscale(S);%將頻譜縮放到[0~255]h=special(‘sobel’);%h是sobel濾波器

PQ=paddedsize(size(f));%計(jì)算填充尺寸

H=freqz2(h,PQ(1),PQ(2));%頻域?yàn)V波器H1=ifftshift(H);%重排濾波器,原點(diǎn)位于左上角gs=imfilter(double(f),h);%空間域?yàn)V波gf=dftfilt(f,H);%頻域?yàn)V波imshow(gs,

[]);

figure

imshow(gf,

[]);五、在頻率域直接生成濾波器頻域?yàn)V波器的網(wǎng)格數(shù)組function[U,V]=dftuv(M,N) u=0:(M-1);%u是0到M-1的數(shù)組

v=0:(N-1);%v是0到N-1的數(shù)組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]是生成的網(wǎng)格五、在頻率域直接生成濾波器則C=[01R=[0001101];2

其中函數(shù)meshgrid用來生成網(wǎng)格數(shù)組，語法：[C,

R]

=

meshgrid(c,r);c和r是輸入的行向量，C和R是輸出的矩陣C和R的維數(shù)為length(c)*length(r)其中C的行是c的副本，R的列是r的副本如c=[0,1];r=[0,1,2];12];五、在頻率域直接生成濾波器頻域?yàn)V波器的網(wǎng)格數(shù)組[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五、低通頻域?yàn)V波器理想低通濾波器(ILP)具有傳遞函數(shù)：

n階巴特沃茲低通濾波器(BLPF)，截止頻率為D0高斯低通濾波器(GLPF)的傳遞函數(shù)為：五、低通頻域?yàn)V波器例4.4低通濾波PQ=paddedsize(size(f));%填充大小

[U,V]=dftuv(PQ(1),PQ(2));%生成網(wǎng)格

D0=0.05*PQ(2);%截止頻率F=fft2(f,PQ(1),PQ(2));%傅里葉變換

H=exp(-(U.^2+V.^2)./(2*(D0^2)));%H是頻域?yàn)V波器g=dftfilt(f,H);%頻域?yàn)V波

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

figure,imshow(g,[]);五、低通頻域?yàn)V波器低通頻域?yàn)V波器生成函數(shù)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五、低通頻域?yàn)V波器低通頻域?yàn)V波器生成函數(shù)H=1./(1+D./D0).^(2*n));case

‘gaussian’H=

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

filter

type’);end六、繪圖函數(shù)繪制線框圖函數(shù)mesh繪制表面圖函數(shù)surf繪制陰影圖函數(shù)shading

interp七、基本的高通濾波器

給定低通濾波器，可以獲得相應(yīng)的高通濾波器：例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,

[]);八、高頻強(qiáng)調(diào)濾波