PS:大三上學(xué)期學(xué)的數(shù)字信號(hào)處理，Matlab大實(shí)驗(yàn)可以自選題，想到老師上課說的IIR、FIR的區(qū)別，借助圖像觀察兩種濾波器的區(qū)別。當(dāng)然，現(xiàn)在大家使用的圖像處理算法是現(xiàn)代濾波器，與經(jīng)典濾波器分析問題的角度不同，但本質(zhì)上還是對(duì)圖像的濾波。本文為我基于Matlab語言實(shí)現(xiàn)的IIR、FIR濾波。
一、實(shí)驗(yàn)任務(wù)
1. ?利用?Matlab?實(shí)現(xiàn)簡單的圖像操作
2.?為圖像加上噪聲，并用Matlab制作FIR、IIR帶限濾波器，觀察處理效果
3.?利用?Simulink?搭建簡單的圖像處理工程
4.?探究?Matlab?和?ModelSim?結(jié)合的FPGA圖像處理仿真平臺(tái)
二、主要實(shí)驗(yàn)儀器及材料?
Windows 10操作系統(tǒng)、Matlab2018a?
三、基于Matlab的圖像濾波
1．圖像的讀入
實(shí)驗(yàn)內(nèi)容：
通過Matlab讀入彩色圖像，并獲得其長，寬等基本信息。
實(shí)驗(yàn)步驟：
在Matlab工程目錄下，保存一幅彩色圖片，命名為`image.jpg`。
?
使用imread函數(shù)可讀入圖像，并返回一個(gè)對(duì)象；使用size函數(shù)可返回函數(shù)的長寬信息；由于圖像要通過濾波器，使用rgb三通道需要分別設(shè)計(jì)三種濾波器，為實(shí)驗(yàn)過程增加了不必要的工作量，所以本次實(shí)驗(yàn)處理灰色圖像，使用rgb2gray函數(shù)可將rgb圖像轉(zhuǎn)化成灰度圖。
%%?圖像的讀取以及轉(zhuǎn)換
rawImg=imread('image.jpg');?%讀取jpg圖像
grayImg=rgb2gray(rawImg);???%生成灰度圖像
[row,col]=size(grayImg);????????%求圖像長寬
生成的灰度圖像如下圖所示：
2．添加噪聲
實(shí)驗(yàn)內(nèi)容：
為原圖像添加噪聲，為驗(yàn)證帶限濾波器和非帶限濾波器（中值濾波、卡爾曼
濾波）的特點(diǎn)、分別加入以下兩類噪聲。
1）多個(gè)單頻率的正弦波噪聲疊加。頻率分別為：350Hz、400Hz、450Hz；?
2）高斯白噪聲；‘
實(shí)驗(yàn)步驟：
首先規(guī)定掃描頻率和掃描時(shí)間間隔
。再分別生成三種頻率的正弦波，其長度需要與圖像的像素點(diǎn)數(shù)匹配，并疊加成噪聲：
????fz1=350;fz2=400;fz3=450;%?三個(gè)噪聲頻率
????noise=0.4*sin(2*pi*fz1*n*T)+...
????????0.7*sin(2*pi*fz2*n*T)+...
????????0.5*sin(2*pi*fz3*n*T);%?噪聲序列
????為了將正弦波噪聲疊加到圖像上，首先需要將二維圖像的數(shù)據(jù)作歸一化處理，然后映射到一維。
normImgMartix=im2double(grayImg);???%?圖像數(shù)據(jù)進(jìn)行歸一化
rawMartix=zeros(1,row*col);????????????????????%?初始化一維矩陣
for?i=1:row
????for?j=1:col
????????rawMartix(col*(i-1)+j)=normImgMartix(i,j);
????end
end??????????%將M*N維矩陣變成1維矩陣
再進(jìn)行灰度圖像和噪聲信號(hào)的疊加，并將一維變換到二維。
????rawMartixWithNoise=rawMartix+noise;%?加入噪聲的序列
????noiseMartix=zeros(row,col);
????%?一維變M*N矩陣
????for?i=1:row
????????for?j=1:col
????????????noiseMartix(i,j)=rawMartixWithNoise(col*(i-1)+j);
????????end
????end
加入單頻混疊干擾的圖像如下圖所示：
再然后，再原始圖像上加入高斯白噪聲，Matlab內(nèi)置了噪聲函數(shù)，調(diào)用語法如下。
whiteNoiseImg=imnoise(grayImg);?%?加高斯噪聲，給simulink用
加入高斯白噪聲的圖像如下圖所示：
?
3．設(shè)計(jì)數(shù)字濾波器
實(shí)驗(yàn)內(nèi)容：
設(shè)計(jì)IIR、FIR數(shù)字低通濾波器。
低通濾波器性能指標(biāo)見代碼。
觀察兩種濾波器對(duì)加入了單頻混疊噪聲的圖像的濾除效果，比較兩者的優(yōu)缺點(diǎn)。再觀察兩種低通濾波器對(duì)加入了高斯白噪聲的圖像的屢出效果，分析原因。
實(shí)驗(yàn)步驟：
首先，設(shè)計(jì)IIR低通濾波器，用直接設(shè)計(jì)數(shù)字濾波器法設(shè)計(jì)巴特沃斯低通濾波器。再用hamming窗設(shè)計(jì)FIR低通濾波器，相關(guān)設(shè)計(jì)代碼如下：
%%?設(shè)計(jì)IIR濾波器并分析相關(guān)指標(biāo)
wp=250*2/fs;ws=300*2/fs;Rp=3;Rs=20;
[Nm,Wc]=buttord(wp,ws,Rp,Rs);
[b,a]=butter(Nm,Wc);?????????
H=freqz(b,a,f*2*pi/fs);
mag=abs(H);pha=angle(H);
mag1=20*log((mag+eps)/max(mag));
%%?設(shè)計(jì)FIR濾波器并分析相關(guān)指標(biāo)
wc=280*2/fs;?
%6dB截止頻率280kHz
fx=[0 wc wc 1];
m=[1 1 0 0];
%理想頻幅響應(yīng)
b1=fir2(40,fx,m,hamming(41));
H1=freqz(b1,1,f*2*pi/fs);
mag2=abs(H1);pha1=angle(H1);
mag3=20*log((mag2+eps)/max(mag2));
觀察設(shè)計(jì)出的兩種濾波器的特性。
?
????分別用IIR、FIR濾波器濾除圖像中的單頻混疊噪聲，并且為了便于觀察圖像的頻譜變化，做一次中心變換。
rawMartixWithNoiseWithIIR=filter(b,a,rawMartixWithNoise);
rawMartixWithNoiseWithIIRFFT=fft(rawMartixWithNoiseWithIIR);
rawMartixWithNoiseWithIIRFFTShift=fftshift(rawMartixWithNoiseWithIIRFFT);
其中，b，a為設(shè)計(jì)出的濾波器的分子分母系數(shù)，將其帶入為IIR、FIR濾波器相應(yīng)系數(shù)即可。觀察其濾除后的圖像。
對(duì)上述結(jié)果作簡單分析：經(jīng)過兩個(gè)低通濾波器，高頻干擾被濾除了。但是圖像有部分失真，體現(xiàn)在圖像的左邊有小塊黑點(diǎn)，這些黑點(diǎn)原本是屬于圖像的右側(cè)。原因在于同一個(gè)像素點(diǎn)，經(jīng)過濾波器會(huì)產(chǎn)生一定的延時(shí)，階數(shù)越高，延時(shí)越大，由于相同指標(biāo)下，F(xiàn)IR的階數(shù)比IIR的階數(shù)要高，所以其產(chǎn)生的延時(shí)越大，效果也越差一點(diǎn)。但是由于FIR是線性相位，每一個(gè)點(diǎn)的延時(shí)都是相同的， 所以很好進(jìn)行修正。下面探討修正的方法。
4．對(duì)濾波后的圖像進(jìn)行相位修正
實(shí)驗(yàn)內(nèi)容：
由于信號(hào)通過兩類濾波器像素點(diǎn)會(huì)產(chǎn)生延時(shí)，會(huì)出現(xiàn)圖像部分像素點(diǎn)偏移的
情況，現(xiàn)在探討如何對(duì)這種情況進(jìn)行修正。利用Matlab自帶的grpdelay函數(shù)能夠獲得濾波器的延時(shí)情況。根據(jù)情況，給出解決方案。
實(shí)驗(yàn)步驟：
????求解濾波器延時(shí)調(diào)用代碼格式為：
grd=grpdelay(b,a,f*2*pi/fs);
上述代碼中，b,a為對(duì)應(yīng)濾波器的分子分母系數(shù)。由于結(jié)果過于長，這里不列出。僅給出結(jié)論。對(duì)于IIR濾波器，不同的像素點(diǎn)輸入后的群延時(shí)不同，對(duì)于FIR濾波器，每一個(gè)像素點(diǎn)的延時(shí)都是定值，在本實(shí)驗(yàn)中，該值為20，所以我們只需要將FIR濾波器處理后的圖像向前推移20位即可。編寫代碼如下：
K=round(mean(grd));
rawMartixWithNoiseWithFIRFixed=[rawMartixWithNoiseWithFIR((K+1):L),rawMartixWithNoiseWithFIR(1:K)];
再將圖像轉(zhuǎn)換成二維矩陣顯示，處理后的結(jié)果如下圖所示：
但是對(duì)于IIR濾波器，由于延時(shí)不同，其失真不可挽回，且對(duì)于每個(gè)像素點(diǎn)都要減去相應(yīng)的延時(shí)。即使能夠勉強(qiáng)恢復(fù)圖像，耗費(fèi)的時(shí)間也是很大的。所以可以得出結(jié)論對(duì)于圖像處理來說FIR濾波器比IIR濾波器更適合。
5．比較濾波前后圖像的頻譜?
實(shí)驗(yàn)內(nèi)容：
在一個(gè)窗口里畫出圖像經(jīng)過Matlab處理過程中的頻譜圖，反應(yīng)頻譜變化。
實(shí)驗(yàn)步驟：
????各狀態(tài)的頻譜圖如下圖所示：
可以看到濾波器工作正常，各頻率成分有被很好的濾除。調(diào)用sound函數(shù)回放聲音，感受聲音的變化。
5．高斯白噪聲的濾除
實(shí)驗(yàn)內(nèi)容：
對(duì)于高斯白噪聲，首先嘗試FIR、IIR濾波器的濾除方法，觀察效果。后結(jié)
合Simulink給出更好的解決方案。
實(shí)驗(yàn)步驟：
???首先，嘗試用IIR、FIR濾波器對(duì)摻入噪聲的圖片進(jìn)行濾波，效果如下圖。?
可以發(fā)現(xiàn)，結(jié)果中看不到原始圖像的特征。分析原因：高斯白噪聲的頻帶覆蓋了整個(gè)頻譜，會(huì)將原始圖像的頻譜淹沒。如果此時(shí)仍然使用帶通濾波器很難再將原有的圖像恢復(fù)。在網(wǎng)絡(luò)上查閱資料不難發(fā)現(xiàn)，高斯白噪聲較好的濾除方法是中值濾波。
中值濾波是基于排序統(tǒng)計(jì)理論的一種能有效抑制噪聲的非線性信號(hào)處理技術(shù)，中值濾波的基本原理是把數(shù)字圖像或數(shù)字序列中一點(diǎn)的值用該點(diǎn)的一個(gè)鄰域中各點(diǎn)值的中值代替，讓周圍的像素值接近的真實(shí)值，從而消除孤立的噪聲點(diǎn)。從其算法原理來看，對(duì)于上述人為添加的高頻混疊干擾來說，中值濾波起不到很好效果，但對(duì)于高斯白噪聲這一類隨機(jī)噪聲，能夠起到較好的過濾作用。
接下來，在SimuLink里搭建圖像處理平臺(tái)與中值濾波環(huán)境。簡單介紹一下使用到的模塊：
Image From Workspace： 獲取工作區(qū)的圖像類型，獲取白噪聲。
Median Filter：中值濾波。
Video Viewer：能夠?qū)D像矩陣可視化顯示的模塊。
在SimuLink中搭建工程如下圖所示：
????運(yùn)行后，結(jié)果如下圖所示：
????可見，相對(duì)于帶限濾波器而言，中值濾波對(duì)白噪聲的處理有更好的效果。
對(duì)應(yīng)的代碼可以去我的個(gè)人博客下載：點(diǎn)擊進(jìn)入
