實驗4DFT及頻域濾波

MACROBUTTONMTEditEquationSection2SEQMTEqn\r\hSEQMTSec\r1\hSEQMTChap\r1\h實驗名稱：數(shù)字圖像DFT及頻域濾波實驗目的了解數(shù)字圖像各種正交變換的概念、原理和用途。深熟練掌握數(shù)字圖像的DFT/DCT的原理、方法和實現(xiàn)流程，熟悉兩種變換的性質(zhì)，并能對數(shù)字圖像DFT及DCT的結(jié)果進行必要解釋。熟悉熟悉和掌握利用MATLAB工具進行數(shù)字圖像FFT及DCT的基本步驟、MATLAB函數(shù)使用及具體變換的處理流程，并能根據(jù)需要進行必要的頻譜分析和可視化顯示。熟悉利用空域濾波器構(gòu)建對應的頻域濾波器的方法和關鍵步驟。熟悉和掌握幾種典型的頻域低通濾波器及高通濾波器的原理、特性和作用。搞清空域圖像處理與頻域圖像處理的異同，包括處理流程、各自的優(yōu)勢等。掌握頻域濾波的基本原理和基本流程，并能編寫出相應的程序代碼。實驗原理FFT實驗傅立葉變換：非周期函數(shù)表示為正弦和/或余弦乘以加權(quán)函數(shù)的積分。一維連續(xù)Fourier變換對函數(shù)f(x)進行傅立葉變換得到F(u) 逆變換，即將F(u)變換到f(x)為 一維離散Fourier變換正變換(DFT) 逆變換(IDFT) 用幅值和相位表示傅立葉變換 用途：頻域濾波，圖像增強，頻譜分析。實驗步驟 模型圖像的FFT實驗圖SEQ圖\*ARABIC1圖SEQ圖\*ARABIC1（a）（b）對變換結(jié)果做頻譜中心化處理，并分別顯示出其2D頻譜圖以及對應的3D頻譜圖。對以上兩幅原始圖像FFT后的頻譜圖進行分析，可以得出什么樣的結(jié)論或驗證了DFT的什么性質(zhì)。圖SEQ圖\*ARABIC2實際圖像的FFT實驗讀任意讀取一幅灰度圖像，對其進行FFT變換，分別畫出變換前的原始圖像及其FFT后頻譜的2D平面圖，要求變換結(jié)果做頻譜中心化處理。畫出對應的3D中心化頻譜。圖SEQ圖\*ARABIC3數(shù)字圖像的頻域濾波處理設定截止頻率D0=100，試分別構(gòu)建256×256的頻域理想低通濾波器（ILPF）和頻域理想高通濾波器（IHPF）。編寫MATLAB代碼，分別畫出它們的頻域濾波器響應3D圖及其對應的2D投影平面圖。任意讀取一幅數(shù)字圖像。編寫MATLAB代碼，分別利用理想低通濾波器（ILPF）和理想高通濾波器（IHPF）對其進行頻域濾波處理。具體可參考如下步驟：消除折疊現(xiàn)象的填充：P=2M，Q=2N，其中M、N為原圖像的尺寸；原圖f(x,y)的傅立葉變換，產(chǎn)生F(u,v)；頻譜中心化處理（當然也可以不變換，則H(u,v)要改變）；用濾波器函數(shù)H(u,v)乘以F(u,v)；傅立葉反變換；取實數(shù)部分，絕對值很小的虛數(shù)部分是浮點運算存在誤差造成的；空間域中心還原變換（反中心化）。分別顯示ILPF和IHPF兩種濾波器濾波前的原始圖像、濾波器頻域響應（2D）及濾波結(jié)果圖像。對濾波結(jié)果進行必要分析，如振鈴現(xiàn)象、圖像模糊、圖像變暗等說明原因，以及盡可能克服這些現(xiàn)象的有效措施等。圖SEQ圖\*ARABIC4實驗結(jié)果及分析圖SEQ圖圖SEQ圖\*ARABIC5圖圖SEQ圖\*ARABIC6圖圖SEQ圖\*ARABIC7分析I1的白色矩形面積小與I，可以認為I1是I的縮小。 而傅立葉有如下的變換性質(zhì)： 從圖（5）和圖（6）可以驗證該性質(zhì)。圖圖SEQ圖\*ARABIC8圖圖SEQ圖\*ARABIC9圖圖SEQ圖\*ARABIC10圖SEQ圖SEQ圖\*ARABIC11圖SEQ圖\*ARABIC12振鈴現(xiàn)象：已知頻域上的方波對應著空域的sinc波，頻域：，傅立葉反變換到時域有：，頻域相乘相當于時域作卷積。圖像經(jīng)過理想低通濾波器后，時域上相當于原始圖像與sinc函數(shù)卷積，由于sinc函數(shù)振蕩，則卷積后圖像也會振蕩；或者說由于sinc函數(shù)有兩個負邊帶，卷積后圖像信號兩側(cè)出現(xiàn)“過沖現(xiàn)象”，而且能量不集中，即產(chǎn)生振鈴效應?？梢允褂酶咚沟屯V波解決該現(xiàn)象。圖像模糊：由于理想低通濾波器，濾去高通成分，削弱了圖像的細節(jié)，從而使圖像變得模糊圖像變暗：由于高通濾波器會濾去低通成分，只突出圖像細節(jié)，所以圖像變暗，為改善這種情況，將高通濾波器的轉(zhuǎn)移函數(shù)加上一常數(shù)量以便引入一些低頻分量，也就是高頻提升濾波器。六.實驗心得體會和建議心得體會：通過這次實驗使我明白了可以利用傅立葉變換將空域濾波轉(zhuǎn)換在頻域濾波，更加符合人們的認知，而且便于在實際中實現(xiàn)。建議：可以讓大家編寫相關代碼解決振鈴現(xiàn)象，圖像變暗。七.程序源代碼模型圖像的FFT實驗fs=1000;I=zeros(256,256);%創(chuàng)建I矩陣I(98:158,98:158)=1;I1=zeros(256,256);%創(chuàng)建I1矩陣I1(118:138,118:138)=1;subplot(121),imshow(I);title('I');%分別在同一幅圖中畫出I和I1矩陣subplot(122),imshow(I1);title('I1');F=fft2(I);F_1=abs(fftshift(F));%分別對I和I1矩陣做二維傅立葉變換，并對結(jié)果中心化處理F1=fft2(I1);F1_1=abs(fftshift(F1));figure,%分別顯示2D頻譜圖subplot(121),imshow(log(1+F_1),[]);title('The2DSpectrumofI');subplot(122),imshow(log(1+F1_1),[]);title('The2DSpectrumofI1');figure,%分別顯示3D頻譜圖subplot(121),mesh(F_1);title('The3DSpectrumofI');xlabel('fx'),ylabel('fy'),zlabel('amplitude');subplot(122),mesh(F1_1);title('The3DSpectrumofI1');xlabel('fx'),ylabel('fy'),zlabel('amplitude');實際圖像的FFT實驗clc,clear,closeall;I=imread('barbara.jpg');I=double(I);I1=fft2(I);I2=abs(fftshift(I1));subplot(121),imshow(I,[]);title('原圖')subplot(122),imshow(log(1+I2),[]);title('頻譜圖')[M,N]=size(I2);m=71;n=71;Fx=(0:m-1)-floor(m/2);Fy=(0:n-1)-floor(n/2);figure[x,y]=meshgrid(Fx,Fy);surfl(x,y,I2(M/2-floor(m/2):M/2+floor(m/2),N/2-floor(n/2):N/2+floor(n/2))/max(I2(:)));title('3D頻域圖'),xlabel('fx'),ylabel('fy'),zlabel('振幅')數(shù)字圖像的頻域濾波處理M=256;N=256;foru=1:2*Mforv=1:2*ND(u,v)=sqrt((u-M)^2+(v-N)^2);endendd0=100;foru=1:2*Mforv=1:2*NifD(u,v)>d0H1(u,v)=0;%H1低通濾波H2(u,v)=1;%H2高通濾波elseH1(u,v)=1;H2(u,v)=0;endendendfigure,%分別顯示ILPF和IHPF的2D投影平面圖subplot(121),imshow(H1,[]);title('2DILPF頻譜圖')subplot(122),imshow(H2,[]);title('2DBLPF頻譜圖')figure,%分別顯示ILPF和IHPF的3D圖subplot(121),mesh(H1);title('3DILPF頻譜圖')xlabel('u'),ylabel('v'),zlabel('振幅')subplot(122),mesh(H2);title('3DBLPF頻譜圖')xlabel('u'),ylabel('v'),zlabel('振幅')I=imread('cameraman.tif');%讀入要處理的圖像I=double(I);Ip=zeros(2*M,2*N);%預分配內(nèi)存Ip(1:M,1:N)=I;%加0處理F=fft2(Ip);%2DfftF0=fftshift(F);%頻域中心化處理F1=F0.*H1;F2=F0.*H2;Fa=real(ifft2(ifftshift(F1)));%傅立葉反變換取實部Fb=real(ifft2(ifftshift(F2)));figure,subplot(231),imshow(I,[]);title('原圖');subplot(232),imshow(log(1+F1),[]);title('中心化頻率譜');subplot(233),imshow(Fa(1:M,1:N),[]);title('ILPF濾波結(jié)果');subplot(235),imshow(log(1+F2),[]);title('中心化頻率譜');subplot(236),imshow(Fb(1:M,1:N),[])