基于matlab的濾波器設(shè)計(jì)

.z光電圖像課程設(shè)計(jì)報(bào)告書(shū)課題名稱基于matlab的濾波器設(shè)計(jì)一、設(shè)計(jì)任務(wù)及要求：設(shè)計(jì)任務(wù)：如何用matlab來(lái)實(shí)現(xiàn)濾波。要求：1.選圖像，并輸入；2.采用butter命令，設(shè)計(jì)一個(gè)巴特沃斯低通濾波器，以此實(shí)現(xiàn)信號(hào)的濾波；3.在同一面上顯示濾波前后的所有的圖像和程序；4.使用MATLAB分別實(shí)現(xiàn)IIR濾波器和FIR濾波器；指導(dǎo)教師簽名：年月日二、指導(dǎo)教師評(píng)語(yǔ)：指導(dǎo)教師簽名：年月日三、成績(jī)驗(yàn)收蓋章年月日?qǐng)D像復(fù)原的MATLAB實(shí)現(xiàn)1課程設(shè)計(jì)目的〔1〕了解基于matlab的濾波器處理及其根本操作；〔2〕學(xué)習(xí)MATLAB在濾波器中的使用；〔3〕提高學(xué)習(xí)與解決問(wèn)題的能力。2課程設(shè)計(jì)根本內(nèi)容2.1濾波器的根本原理設(shè)計(jì)數(shù)字濾波器的任務(wù)就是尋求一個(gè)因果穩(wěn)定的線性時(shí)不變系統(tǒng)，并使系統(tǒng)函數(shù)H〔z〕具有指定的頻率特性。數(shù)字濾波器從實(shí)現(xiàn)的網(wǎng)絡(luò)構(gòu)造或者從單位沖激響應(yīng)分類，可以分成無(wú)限單位沖激響應(yīng)〔IIR〕數(shù)字濾波器和有限長(zhǎng)單位沖激響應(yīng)〔FIR〕數(shù)字濾波器。數(shù)字濾波器頻率響應(yīng)的三個(gè)要素：幅度平方響應(yīng)相位響應(yīng)群時(shí)延響應(yīng)IIR數(shù)字濾波器：IIR數(shù)字濾波器的系統(tǒng)函數(shù)為有理分?jǐn)?shù)，即IIR數(shù)字濾波器的逼近問(wèn)題就是求解濾波器的系數(shù)和，使得在規(guī)定的物理意義上逼近所要求的特性的問(wèn)題。如果是在s平面上逼近，就得到模擬濾波器，如果是在z平面上逼近，則得到數(shù)字濾波器。FIR數(shù)字濾波器：設(shè)FIR的單位脈沖響應(yīng)h〔n〕為實(shí)數(shù)，長(zhǎng)度為N，則其z變換和頻率響應(yīng)分別為按頻域采樣定理FIR數(shù)字濾波器的傳輸函數(shù)H(z)和單位脈沖響應(yīng)h〔z〕可由它的N歌頻域采值H(k)唯一確定。MATLAB中提供了幾個(gè)函數(shù)，分別用于實(shí)現(xiàn)IIR濾波器和FIR濾波器。卷積函數(shù)conv，調(diào)用格式為，c=conv〔a，b〕該格式可以計(jì)算兩向量a和b的卷積，可以直接用于對(duì)有限長(zhǎng)信號(hào)采用FIR濾波器的濾波。函數(shù)filter的調(diào)用格式為，y=filter〔b，a，*〕該格式采用數(shù)字濾波器對(duì)數(shù)據(jù)進(jìn)展濾波，既可以用于IIR濾波器，也可以用于FIR濾波器。其中向量b和a分別表示系統(tǒng)函數(shù)的分子，分母多項(xiàng)式的系數(shù)，假設(shè)a=1，此時(shí)表示FIR濾波器，否則就是IIR濾波器。該函數(shù)就是利用給出的向量b和a，對(duì)*中的數(shù)據(jù)進(jìn)展濾波，結(jié)果放入向量y。函數(shù)fftfilt的調(diào)用格式為，y=fftfilt〔b，*〕該格式是利用基于FFT的重疊相加法對(duì)數(shù)據(jù)進(jìn)展濾波，這種頻域?yàn)V波技術(shù)只對(duì)FIR濾波器有效。該函數(shù)是通過(guò)向量b描述的濾波器對(duì)*數(shù)據(jù)進(jìn)展濾波。關(guān)于用butter函數(shù)求系統(tǒng)函數(shù)分子與分母的幾種形式。[b，a]=butter〔N,wc，‘high’〕：設(shè)計(jì)N階高通濾波器，wc為它的3dB邊緣頻率，以Hz為單位。[b，a]=butter〔N,wc〕：當(dāng)wc為具有兩個(gè)元素的矢量wc=[w1，w2]時(shí)，它設(shè)計(jì)2N階帶通濾波器，3dB通帶，以Hz為單位。[b，a]=butter〔N,wc，‘stop’〕：假設(shè)wc=[w1，w2]，則它設(shè)計(jì)2N階帶阻濾波器，3dB通帶，以Hz為單位。如果在這個(gè)函數(shù)輸入變?cè)淖詈?，加一個(gè)變?cè)皊〞，設(shè)計(jì)的是模擬濾波器。這里不作討論。為了設(shè)計(jì)任意的選項(xiàng)巴特沃斯濾波器，必須知道階數(shù)N和3dB邊緣頻率矢量wc。這可以直接利用信號(hào)處理工具箱中的buttord函數(shù)來(lái)計(jì)算。如果濾波器指標(biāo)，，和，則調(diào)用格式為[N,wc]=buttord〔wp，ws，Rp，As〕對(duì)于不同類型的濾波器，參數(shù)wp和ws有一些限制：對(duì)于低通濾波器，wp<ws;對(duì)于高通濾波器，wp>ws;對(duì)于帶通濾波器，wp和ws分別為具有兩個(gè)元素的矢量，wp=[wp1,wp2]和ws=[ws1,ws2]，并且ws1<wp1<wp2<ws2;對(duì)于帶阻濾波器wp1<ws1<ws2<wp2。2.2濾波器與MATLAB簡(jiǎn)介濾波器例如：在這里為了說(shuō)明如何用matlab來(lái)實(shí)現(xiàn)濾波，特舉出一個(gè)簡(jiǎn)單的函數(shù)信號(hào)濾波實(shí)例(對(duì)信號(hào)*(n)=sin(n/4)+5cos(n/2)進(jìn)展濾波，信號(hào)長(zhǎng)度為500點(diǎn))，從中了解濾波的實(shí)現(xiàn)過(guò)程。程序如下：Wn=0.2*pi;N=5;[b,a]=butter(N,Wn/pi);N=0:1:499;*=sin(pi*n/4)+5*cos(pi*n/2);*=fft(*,4096);Subplot(2,2,1);plot(*);title(‘濾波前信號(hào)的波形’);Subplot(2,2,2);plot(*);title(‘濾波前信號(hào)的頻譜’);y=filter(b,a,*);Y=fft(y,4096);Subplot(2,2,3);plot(y);title(‘濾波后信號(hào)的波形’);Subplot(2,2,4);plot(Y);title(‘濾波后信號(hào)的頻譜’);結(jié)果如圖：在這里，是采用了butter命令，設(shè)計(jì)出一個(gè)巴特沃斯低通濾波器，從頻譜圖中可以很明顯的看出來(lái)。下面，也就是本課題的主要內(nèi)容，也都是運(yùn)用到了butter函數(shù)，以便容易的得到系統(tǒng)函數(shù)的分子與分母系數(shù)，最終以此來(lái)實(shí)現(xiàn)信號(hào)的濾波。N階高通濾波器的設(shè)計(jì)〔在這里，以5階為例，其中wc為其3dB邊緣頻率，以Hz為但位〕，程序設(shè)計(jì)如下：*=wavread(‘ding.wav’);Sound(*);N=5;wc=0.3;[b,a]=butter(N,wc,’high’);*=fft(*);Subplot(3,2,1);plot(*);title(‘濾波前信號(hào)的波形’);Subplot(3,2,2);plot(*);title(‘濾波前信號(hào)的頻譜’);y=filter(b,a,*);Y=fft(y);Subplot(3,2,3);plot(y);title(‘IIR濾波后信號(hào)的波形’);Subplot(3,2,4);plot(Y);title(‘IIR濾波后信號(hào)的頻譜’);z=fftfilt(b,*);Z=fft(z);Subplot(3,2,5);plot(z);title(‘FIR濾波后信號(hào)的波形’);Subplot(3,2,6);plot(Z);title(‘FIR濾波后信號(hào)的頻譜’);得到結(jié)果如圖：〔3〕2N階帶通濾波器的設(shè)計(jì)〔在這里，以10階為例，其中wc為其3dB邊緣頻率，以Hz為單位，在wc=[w1,w2],w1wcw2〕,程序設(shè)計(jì)如下：*=wavread(‘ding.wav’);Sound(*);N=5;wc=[0.3,0.6];[b,a]=butter(N,wc);*=fft(*);Subplot(3,2,1);plot(*);title(‘濾波前信號(hào)的波形’);Subplot(3,2,2);plot(*);title(‘濾波前信號(hào)的頻譜’);y=filter(b,a,*);Y=fft(y);Subplot(3,2,2);plot(*);title(‘濾波前信號(hào)的頻譜’);y=filter(b,a,*);Y=fft(y);Subplot(3,2,3);plot(y);title(‘IIR濾波后信號(hào)的波形’);Subplot(3,2,4);plot(Y);title(‘IIR濾波后信號(hào)的頻譜’);z=fftfilt(b,*);Z=fft(z);Subplot(3,2,5);plot(z);title(‘FIR濾波后信號(hào)的波形’);Subplot(3,2,6);plot(Z);title(‘FIR濾波后信號(hào)的頻譜’);得到的結(jié)果如圖：〔4〕關(guān)于MATLABMATLAB對(duì)于技術(shù)計(jì)算來(lái)說(shuō)是一種高性能的語(yǔ)言。它以易于應(yīng)用的環(huán)境集成了計(jì)算、可視化和編程，在該環(huán)境下，問(wèn)題及其解以我們熟悉的數(shù)學(xué)表示法來(lái)表示。典型的應(yīng)用包括如下方面：數(shù)學(xué)和計(jì)算/算法開(kāi)發(fā)/數(shù)據(jù)獲取/建模、模擬和原型設(shè)計(jì)/數(shù)據(jù)分析、研究和可視化/科學(xué)和工程圖形/應(yīng)用開(kāi)發(fā)，包括圖像用戶界面構(gòu)建。圖2約束最小二乘濾波恢復(fù)比照?qǐng)D如圖2所示，利用振鈴抑制恢復(fù)圖像是3幾種中恢復(fù)效果最好的，其他幾種方法也可以恢復(fù)但是比擬模糊，效果不是很明顯。使用Lucy-Richardson算法的迭代非線性復(fù)原Lucy-Richardson(LR)算法假設(shè)圖像服從Poission分布，采用最大似然法進(jìn)展估計(jì)，是一種基于貝葉斯分析的迭代算法。Lucy-Richardson算法是目前世界上應(yīng)用最廣泛的函數(shù)恢復(fù)技術(shù)之一，它是一種迭代非線性復(fù)原算法，能夠按照泊松噪聲統(tǒng)計(jì)標(biāo)準(zhǔn)求出與給定的點(diǎn)擴(kuò)散函數(shù)PSF(PointSpreadFunction)卷積后最有可能成為輸入模糊圖像的圖像。當(dāng)PSF但圖像噪聲信息未知時(shí)，也可以使用這個(gè)函數(shù)進(jìn)展有效的恢復(fù)。其最優(yōu)估計(jì)以最大似然準(zhǔn)則作為標(biāo)準(zhǔn)，即要使概率密度函數(shù)最大，推導(dǎo)出的迭代式為：(1)其中，和分別為卷積運(yùn)算和相關(guān)運(yùn)算；為迭代次數(shù)，可以令進(jìn)展迭代，可以證明，當(dāng)噪聲可以忽略當(dāng)不斷增大時(shí)會(huì)依概率收斂于，從而恢復(fù)出原始圖像。當(dāng)噪聲不可忽略時(shí)，把式(1.2)代入式(7.1)可得到：(2)從上式可看出，假設(shè)噪聲不可忽略，則以上過(guò)程的收斂性將難以保證，即LR存在放大噪聲的缺陷。因此，處理噪聲項(xiàng)是LR算法應(yīng)用于低信噪比圖像復(fù)原的關(guān)鍵。是未退化圖像的估計(jì)，這個(gè)算法的迭代本質(zhì)是顯而易見(jiàn)的。它的非線性本質(zhì)是載方程的右邊用來(lái)除產(chǎn)生的。就像大多數(shù)非線性方法一樣，關(guān)于什么時(shí)候停頓L-R算法通常是很難答復(fù)。通常，接下來(lái)的處理途徑是對(duì)于給定的應(yīng)用，在我們獲得滿意的結(jié)果時(shí)，觀察輸出并終止算法。在IPT中，L-R算法是由名為deconvlucy的函數(shù)完成的，該函數(shù)的語(yǔ)法為fr=deconvlucy(g,PSF,NUMIT,DAMPAR,WEIGHT)其中，fr代表復(fù)原的圖像，g代表退化的圖像，PSF是點(diǎn)擴(kuò)散函數(shù)，NUMIT為迭代次數(shù)〔默認(rèn)為10次〕，DAMPAR和WEIGHT定義如下。DAMPAR是一個(gè)標(biāo)量，它指定了結(jié)果圖像與原圖像g之間的偏離閾值。當(dāng)像素偏離原值的*圍在DAMPAR之內(nèi)時(shí)，就不用再迭代。這既抑制了這些像素上的噪聲，又保存了必要的圖像細(xì)節(jié)。默認(rèn)值為0〔無(wú)衰減〕。WEIGHT是一個(gè)與g同樣大小的數(shù)組，它為每一個(gè)像素分配一個(gè)權(quán)重來(lái)反映其重量。當(dāng)用一個(gè)指定的PSF來(lái)模擬模糊時(shí)，WEIGHT可以從計(jì)算像素中剔除那些來(lái)自圖像邊界的像素點(diǎn)，因此，PSF造成的模糊是不同的。假設(shè)PSF的大小為，則在WEIGHT中用到的零邊界的寬度是ceil〔n/2〕。默認(rèn)值是同輸入圖像g同等大小的一個(gè)單位數(shù)組。假設(shè)復(fù)原圖像呈現(xiàn)出由算法中所用的離散傅里葉變換所引入的振鈴，則在調(diào)用函數(shù)deconvlucy之前，要利用函數(shù)edgetaper。下面通過(guò)程序事例說(shuō)明Lucy-Richardson算法恢復(fù)圖像：I=imread('G:\1.jpg');PSF=fspecial('gaussian',5,5);Blurred=imfilter(I,PSF,'symmetric','conv');V=.003;BN=imnoise(Blurred,'gaussian',0,V);luc=deconvlucy(BN,PSF,5);figuresubplot(2,2,1);imshow(I);title('原始圖像');subplot(2,2,2);imshow(Blurred);title('模糊后的圖像');subplot(2,2,3);imshow(BN);title('加噪后的圖像');subplot(2,2,4);imshow(luc);title('恢復(fù)后的圖像');圖3Lucy-Richardson算法恢復(fù)圖像比照?qǐng)D用Lucy-Richardson算法可以較好的恢復(fù)圖像[9]。但由于迭代產(chǎn)生的噪聲痕跡是最大化可能性數(shù)據(jù)逼近法的常見(jiàn)問(wèn)題，在低信噪比條件下，恢復(fù)圖像可能會(huì)出現(xiàn)一些斑點(diǎn)，這些斑點(diǎn)并不代表圖像的真實(shí)構(gòu)造，只不過(guò)是恢復(fù)圖像過(guò)于逼近噪聲所產(chǎn)生的結(jié)果。另外本方法存在一些較嚴(yán)重的缺陷，問(wèn)題一是噪聲放大問(wèn)題，這也是這類方法〔利用最大似然性求解〕的通病。問(wèn)題二是對(duì)于恢復(fù)圖像中的不同局部，分別執(zhí)行多少迭代才適宜的問(wèn)題。因?yàn)閳D像**噪比高的局部可能需要數(shù)百次迭代才能獲得滿意的結(jié)果；而另一些光滑的對(duì)象可能只需很少次數(shù)即可到達(dá)滿意的結(jié)果，所以適中選擇迭代次數(shù)對(duì)圖像恢復(fù)也很重要。這兩個(gè)問(wèn)題假設(shè)得不到解決，將會(huì)對(duì)最終結(jié)果產(chǎn)生不利影響從圖3我們可以看到，經(jīng)屢次迭代，尤其是在低信噪比情況下，重建圖像可能會(huì)出現(xiàn)一些斑點(diǎn)，這些斑點(diǎn)并不代表圖像的真實(shí)構(gòu)造，是輸出圖像過(guò)于逼近噪聲所產(chǎn)生的結(jié)果。用Lucy-Richardson算法可以較好的恢復(fù)圖像。但由于迭代產(chǎn)生的噪聲痕跡是最大化可能性數(shù)據(jù)逼近法的常見(jiàn)問(wèn)題，在低信噪比條件下，恢復(fù)圖像可能會(huì)出現(xiàn)一些斑點(diǎn)，這些斑點(diǎn)并不代表圖像的真實(shí)構(gòu)造，只不過(guò)是恢復(fù)圖像過(guò)于逼近噪聲所產(chǎn)生的結(jié)總結(jié)本文介紹了圖像退化的原因并且簡(jiǎn)要介紹了當(dāng)前主流的圖像復(fù)原方法，并通過(guò)對(duì)各種復(fù)原方法的仿真，了解了各種方法的優(yōu)劣性，為我們?cè)趯?shí)際生活提供依據(jù)。但是無(wú)論是哪一種方法都有所局限性，我們應(yīng)該努力致力于研究新型的優(yōu)秀的圖像復(fù)原方法，來(lái)獲得更好的圖像復(fù)原效果。同時(shí)，我們知道總而言之，圖像恢復(fù)的算法很多，無(wú)論何