MATLAB--数字图像处理 频域图像分析
一、實驗名稱
頻域圖像分析
二、實驗目的
1.熟悉MATLAB軟件的使用。
2.掌握頻域圖像分析的原理及數學運算。
三、實驗內容
1.自選一幅圖像,并對其分別添加一定強度的周期噪聲和高斯噪聲,然后分別采用高斯模板、中值濾波的時域方法以及傅里葉變換和小波變換的頻率濾波方法對該含噪圖像進行去噪處理,并基于PSNR值和視覺效果這兩個指標來比較這四種濾波方法對兩種不同噪聲的去噪能力。
2.編寫一個程序,要求實現下列算法:首先將閣像分割為8x8的子圖像,對每個予圖像進行FFT.對每個了圖像中的64個系數。按照每個系數的方差來排序后,舍去小的變換系數,只保留16個系數,實現4: I的圖像壓縮。
3.給定一幅行和列都為2的整數次幕圖像,用Haar小波基函數對其進行二維小波變換,試著將最低尺度近似分量置零再反變換,結果是什么?如果把垂直方向的細節分量置零,反變換后結果又是什么呢?試解釋一下原因。
4.基于小波變換對圖像進行不同壓縮比的壓縮。在同壓縮比情況下,對于基于小波變換和基于傅里葉變換的壓縮結果,比較=二者保留原圖像能里百分比情況。
四、實驗儀器與設備
Win10 64位電腦
MATLAB R2017a
五、實驗原理
1.傅里葉變換
從純粹的數學意義上看,傅里葉變換是將一個函數轉換為一系列周期函數來處理的。從物理效果看,傅里葉變換是將圖像從空間域轉換到頻率域,其逆變換是將圖像從頻率域轉換到空間域。換句話說,傅里葉變換的物理意義是將圖像的灰度分布函數變換為圖像的頻率分布函數。
傅里葉逆變換是將圖像的頻率分布函數變換為灰度分布函數傅里葉變換以前,圖像(未壓縮的位圖)是由對在連續空間(現實空間)上的采樣得到一系列點的集合,通常用一個二維矩陣表示空間上各點,記為z=f(x,y)。又因空間是三維的,圖像是二維的,因此空間中物體在另一個維度上的關系就必須由梯度來表示,這樣我們才能通過觀察圖像得知物體在三維空間中的對應關系。
2.小波變換
小波變換是時間(空間)頻率的局部化分析,它通過伸縮平移運算對信號(函數)逐步進行多尺度細化,最終達到高頻處時間細分,低頻處頻率細分,能自動適應時頻信號分析的要求,從而可聚焦到信號的任意細節帶噪聲信號經過預處理,然后利用小波變換把信號分解到各尺度中,在每一尺度下把屬于噪聲的小波系數去掉,保留并增強屬于信號的小波系數,最后再經過小波逆變換回復檢測信號。
小波變換在去除噪聲時可提取并保存對視覺起主要作用的邊緣信息,而傳統的基于傅里葉變換去除噪聲的方法在去除噪聲和邊沿保持上存在著矛盾,因為傅里葉變換方法在時域不能局部化,難以檢測到局域突變信號,在去除噪聲的同時,也損失了圖像邊沿信息。由此可知,與傅里葉變換去除噪聲的方法相比較,小波變換法去除噪聲具有明顯的性能優勢。
3.PSNR算法
peak的中文意思是頂點。而ratio的意思是比率或比列的。整個意思就是到達噪音比率的頂點信號,psnr一般是用于最大值信號和背景噪音之間的一個工程項目。通常在經過影像壓縮之后,輸出的影像都會在某種程度與原始影像不同。為了衡量經過處理后的影像品質,我們通常會參考PSNR值來衡量某個處理程序能否令人滿意。它是原圖像與被處理圖像之間的均方誤差相對于(2n-1)2的對數值(信號最大值的平方,n是每個采樣值的比特數),它的單位是dB。
六、實驗過程及代碼
t=imread('a1.jpg');%添加高斯噪聲 t1=imnoise(t,'gaussian',0,0.02); [m,n]=size(t); t2=t;%添加周期噪聲 for i=1:m for j=1:n t2(i,j)=t(i,j)+30*sin(30*i)+30*sin(30*j); end end imshow(t1),title('加入高斯噪聲后') figure,imshow(t2),title('加入周期噪聲后');%進行高斯濾波、中值濾波 t3_1=t; t3_1(:,:,1)=medfilt2(t1(:,:,1),[3 3]); t3_1(:,:,2)=medfilt2(t1(:,:,2),[3 3]); t3_1(:,:,3)=medfilt2(t1(:,:,3),[3 3]); figure,imshow(t3_1),title('對高斯噪聲進行中值濾波','fontsize',16) h=fspecial('gaussian',5,3);%確定濾波方式為高斯濾波 5是模板大小 3是方差 t3_2=imfilter(t2,h);%濾波操作 figure,imshow(t3_2),title('對周期噪聲進行高斯濾波','fontsize',16)t3_1=t; t3_1(:,:,1)=medfilt2(t2(:,:,1),[3 3]); t3_1(:,:,2)=medfilt2(t2(:,:,2),[3 3]); t3_1(:,:,3)=medfilt2(t2(:,:,3),[3 3]); figure,imshow(t3_1),title('對高斯噪聲進行中值濾波','fontsize',16) h=fspecial('gaussian',5,3);%確定濾波方式為高斯濾波 5是模板大小 3是方差 t3_2=imfilter(t1,h);%濾波操作 figure,imshow(t3_2),title('對周期噪聲進行高斯濾波','fontsize',16)%對周期噪聲傅里葉變換濾波 t_f=rgb2gray(t2); %將圖像灰度化 t_f=fft2(double(t_f));%利用fft2()函數將圖像從時域空間轉換到頻域空間 t_f=fftshift(t_f);%將零頻平移到中心位置 [m,n]=size(t_f); m_min=round(m/2); n_min=round(n/2); t_rf=t_f; d_0=100;%設置閾值 for i=1:m for j=1:n d_1=sqrt((i-m_min)^2+(j-n_min)^2); if(d_1>d_0) x=0; else x=1; end t_rf(i,j)=x*t_f(i,j); end end t_rf=ifftshift(t_rf); t_rf=uint8(real(ifft2(t_rf))); figure,imshow(t_rf),title('對周期噪聲進行傅里葉變換濾波后')%對高斯噪聲進行傅里葉變換濾波 t_f=rgb2gray(t1); t_f=fft2(double(t_f)); t_f=fftshift(t_f); [m,n]=size(t_f); m_min=round(m/2); n_min=round(n/2); t_rf=t_f; d_0=100; for i=1:m for j=1:n d_1=sqrt((i-m_min)^2+(j-n_min)^2); if(d_1>d_0) x=0; else x=1; end t_rf(i,j)=x*t_f(i,j); end end t_rf=ifftshift(t_rf); t_rf=uint8(real(ifft2(t_rf))); figure,imshow(t_rf),title('對高斯噪聲進行傅里葉變換濾波后')%對周期噪聲進行小波變換濾波 [c,s]=wavedec2(t2,2,'sym4'); a1=wrcoef2('a',c,s,'sym4'); figure,imshow(uint8(a1)); title('對周期噪聲進行小波變換濾波') %對高斯噪聲進行小波變換濾波 [c,s]=wavedec2(t1,2,'sym4'); a1=wrcoef2('a',c,s,'sym4'); figure, imshow(uint8(a1)); title('對高斯噪聲進行小波變換濾波')SNRP算法
function PSNR = psnr(f1, f2) %計算兩幅圖像的峰值信噪比 k = 8; %k為圖像是表示地個像素點所用的二進制位數,即位深。 fmax = 2.^k - 1; a = fmax.^2; MSE =(double(im2uint8(f1)) -double( im2uint8(f2))).^2; b = mean(mean(MSE)); PSNR = 10*log10(a/b);對圖像進行4:1的壓縮t=imread('a6.jpg'); t=rgb2gray(t);%灰度化 [k,p]=size(t); t=double(t)/255;%歸一化 便于計算 %顯示原圖 imshow(t),title('原圖','fontsize',16); %利用blkproc 進行分塊 并對每一塊進行fft操作 t_fft=blkproc(t,[8 8],'fft2(x)'); %利用im2col進行優化操作 便于計算 t_block=im2col(t_fft,[8 8],'distinct'); [t_change,ix]=sort(t_block);%對每一塊圖像進行排序 [m,n]=size(t_block); nums=48; %對后48位系數清零 for i=1:nt_block(ix(1:nums),i)=0; end t_rchange=col2im(t_block,[8 8],[k p],'distinct'); t_ifft=blkproc(t_rchange,[8 8],'ifft2(x)');%對每一塊進行傅里葉反變換 figure, imshow(t_ifft),title('4:1壓縮后','fontsize',16);haar基函數進行小波變換
%實現一維Haar變換 function [h]=D1Haar(f,N) J=log2(N); Y=zeros(J,N); f1=zeros(1,N); H=zeros(1,N); q=f';%調用倒序函數 f1=Reverse(q,N); for A=1:NY(1,A)=f1(A); end%第一層的迭代 for C=1:N/2Y(2,C)=Y(1,C)+Y(1,C+N/2);Y(2,C+N/2)=Y(1,C)-Y(1,C+N/2); end%余下層的迭代 M=0; for B=2:JK=2*N/(2^B); F=zeros(1,K);for C=1:K/2Y(B+1,C)=Y(B,C)+Y(B,C+K/2);Y(B+1,C+K/2)=Y(B,C)-Y(B,C+K/2);endfor C=K+1:2*KF(1,C-K)=Y(B,C);endZ=Reverse(F,K);for D=1:KY(B+1,D+K)=Z(D);endfor C=2*K+1:NY(B+1,C)=Y(B,C); end end%系數修正 H(1,1)=(1/N)*Y(J+1,1); H(1,2)=(1/N)*Y(J+1,2); for a=2:Jb=2^(a-1);c=b^(1/2);for d=b+1:2*bH(1,d)=(c/N)*Y(J+1,d);end end h=H';%二維哈爾函數的正變換 function [y]=D2Haar(x,M,N) %先進行按列變換 for i=1:Nz=zeros(M,1);for c=1:Mz(c,1)=x(c,i);endz1=D1Haar(z,M);%調用一維哈爾變換函數D1Haarfor c=1:Mx1(c,i)=z1(c,1);end end%再進行按行變換 x2=x1';%將變換好的矩陣進行轉置 for j=1:Mz=zeros(N,1);for c=1:Nz(c,1)=x2(c,j);endz1=D1Haar(z,N);%調用一維哈爾變換函數D1Haarfor c=1:Ny1(c,j)=z1(c,1);end end y=y1';%再最后將變換好的矩陣作轉置%對矩陣作一次哈爾變換 function [y]=D2Har(x,M,N) for i=1:Nz=zeros(M,1);for c=1:Mz(c,1)=x(c,i);endz1=D1Haar(z,M);for c=1:Mx1(c,i)=z1(c,1);end end y=x1;%哈爾函數的逆變換 function [y]=ID2Haar(X,M,N) har1=zeros(M); I1=eye(M);%產生M維單位陣 har1=D2Har(I1,M,M);%調用D2Har產生M維哈爾矩陣 har2=M*har1'; I2=eye(N);%產生N維單位陣 har3=zeros(N); har3=D2Har(I2,N,N);%調用D2Har產生M維哈爾矩陣 har4=N*har3; y=har2*X*har4;%哈爾逆變換%求倒敘函數 function [F]=Reverse(F,M) N=log10(M)/log10(2); Y=zeros(1,M); for x=0:M-1A=dec2bin(x,N);%十進制轉二進制B=fliplr(A);%二進制倒敘C=bin2dec(B);%二進制轉十進制Y(x+1)=F(C+1); end for x=0:M-1F(x+1)=Y(x+1); end %對圖片進行哈爾正變化,并對其進行恢復 clear; M=256; N=256;%讀取圖象 X=imread('C:\Users\LiCongliang\Desktop\數字圖像處理\數字圖像處理第七次作業\tea.jpg'); subplot(2,2,1); imshow(X); title('原圖像');%縮放原圖象 x=imresize(X,[M,N]);%將原圖象縮放成分辨率為256*256的圖象 subplot(2,2,2); imshow(x); title('縮放圖象');%對縮放圖象進行Haar變換 y=D2Haar(x,M,N); subplot(2,2,3); imshow(y); title('變換圖象');%對變換后的圖象進行Harr逆變換,恢復原圖象 z=ID2Haar(y,M,N); %Z=imresize(z,[480,360]); subplot(2,2,4); imshow(z,[0,255]); title('恢復圖象');小波變換進行圖像壓縮 X=imread('a5.jpg'); X=rgb2gray(X); subplot(221); imshow(X); title('原始圖像'); %對圖像用小波進行層小波分解 [c,s]=wavedec2(X,2,'haar'); %提取小波分解結構中的一層的低頻系數和高頻系數 cal=appcoef2(c,s,'haar',1); ch1=detcoef2('h',c,s,1); %水平方向 cv1=detcoef2('v',c,s,1); %垂直方向 cd1=detcoef2('d',c,s,1); %斜線方向 %各頻率成份重構 a1=wrcoef2('a',c,s,'haar',1); h1=wrcoef2('h',c,s,'haar',1); v1=wrcoef2('v',c,s,'haar',1); d1=wrcoef2('d',c,s,'haar',1); c1=[a1,h1;v1,d1]; subplot(222),imshow(c1,[]); title ('分解后低頻和高頻信息'); %進行圖像壓縮 %保留小波分解第一層低頻信息 %首先對第一層信息進行量化編碼 ca1=appcoef2(c,s,'haar',1); ca1=wcodemat(ca1,440,'mat',0); %改變圖像高度并顯示 ca1=0.5*ca1; subplot(223);imshow(cal,[]); title('第一次壓縮圖像'); %保留小波分解第二層低頻信息進行壓縮 ca2=appcoef2(c,s,'haar',2); %首先對第二層信息進行量化編碼 ca2=wcodemat(ca2,440,'mat',0); %改變圖像高度并顯示 ca2=0.25*ca2; subplot(224);imshow(ca2,[]); title('第二次壓縮圖像');七、實驗結果與分析
圖 1原圖
1.加入周期噪聲、高斯噪聲
2.對添加了高斯噪聲和周期噪聲的圖像進行高斯濾波
PSNR值
1.對高斯噪聲進行高斯濾波后 23.0287
2.對周期噪聲進行高斯濾波后 23.4837
2.中值濾波
PSNR值:
1.對高斯噪聲進行中值濾波 23.9931
2.對周期噪聲進行中值濾波 24.3134
3.傅里葉變換濾波
PSNR值:
1.對添加了高斯噪聲的圖像進行傅里葉變換濾波 20.4922
2.對添加了周期噪聲的圖像進行傅里葉變換濾波 18.9736
4.小波變換濾波
PSNR值:
1.對添加了高斯噪聲的圖像進行小波變換濾波 23.4712
2.對添加了周期噪聲的圖像進行小波變換濾波 24.4525
分析:
對于高斯噪聲,高斯濾波和傅里葉變換濾波聲的除噪效果較好,中值濾波效果較差,小波變換濾波的處理效果也比較好
對于周期噪聲,中值濾波和高斯濾波效果不是很好,傅里葉變換變換濾波對噪聲的去處效果比較好,對于原圖像損壞不大,小波變換對原圖的損壞較大,但是圖片可以看出噪聲也去除的比較好。
5.圖像壓縮(4:1壓縮) 原圖-左 壓縮后-右
分析:
圖像壓縮算法就是先將一副圖像分成很多小塊,然后分別對這些小塊進行變換,這里采用的是傅里葉變換,然后過濾掉冗余的像素點,然后再利用反變換得到壓縮后的圖像即可。
小波變換
1.定義
小波變換是時間(空間)頻率的局部化分析,它通過伸縮平移運算對信號(函數)逐步進行多尺度細化,最終達到高頻處時間細分,低頻處頻率細分,能自動適應時頻信號分析的要求,從而可聚焦到信號的任意細節帶噪聲信號經過預處理,然后利用小波變換把信號分解到各尺度中,在每一尺度下把屬于噪聲的小波系數去掉,保留并增強屬于信號的小波系數,最后再經過小波逆變換回復檢測信號。
2.優點
小波變換在去除噪聲時可提取并保存對視覺起主要作用的邊緣信息,而傳統的基于傅里葉變換去除噪聲的方法在去除噪聲和邊沿保持上存在著矛盾,因為傅里葉變換方法在時域不能局部化,難以檢測到局域突變信號,在去除噪聲的同時,也損失了圖像邊沿信息。由此可知,與傅里葉變換去除噪聲的方法相比較,小波變換法去除噪聲具有明顯的性能優勢。
Haar基函數進行小波變換
圖 2原圖
圖 3 haar變換
圖 4 haar反變換后
圖 5 最低分量近似置零
圖 6 垂直分量置零
小波變換進行圖像壓縮與傅里葉變換壓縮對比
1.壓縮比 1:2(左-小波壓縮 右-傅里葉壓縮)
2.壓縮比 1:4(左-小波壓縮 右-傅里葉壓縮)
八、實驗總結及心得體會
通過這次實驗,學到了很多。特別是在傅里葉變換和小波變換等方面,開始的時候連傅里葉變換的基礎基礎也不懂,后來在csdn上看了一篇講解傅里葉變換的文章,豁然開朗,傅里葉變換居然可以將一個時域信號轉化到頻域,而且自己還對與i有了更加深刻的理解。雖然傅里葉變換可以把信號從時域轉換到頻域,但是頻域與時域的對應關系卻無法一一對應,所以誕生了小波變換。小波變換的特別之處就是可以把一個時域上的信息轉換為時域-頻域一一對應,這對應特殊信號的提取是有很好的效果,在一定程度上比傅里葉變換更厲害。但是在傅里葉、小波等基礎概念知識方面,自己還是涉獵的比較少,原理的論證公式太復雜了。
更多
獲取更多資料、代碼,微信公眾號:海轟Pro
回復 海轟 即可
總結
以上是生活随笔為你收集整理的MATLAB--数字图像处理 频域图像分析的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 用友nc java启动不了_用友NC开发
- 下一篇: 使用Wireshark抓包分析TCP协议