基于matlab的圖像處理濾波器設計

1、 數字信號處理 課程設計任務書20112012學年 第一學期第 15 周 19 周 題目基于matlab的圖像處理濾波器設計內容及要求：1、設計一個低通FIR濾波器和一個低通IIR濾波器；2、分析比較上述兩種濾波器的優(yōu)劣；3、分析上述兩個模型的幅頻特性、相頻特性、相延遲、群延遲。進度安排： 1、任務分配、查閱資料 2天 2、方案論證 3天 3、分析、設計、調試程序 5天 4、書寫、整理實驗報告和小結 3天成員組成：09044106蘇青文 08陳舒龍（組長） 09方雪松 指導時間： 指導地點： F 618 任務下達2011年12 月 16 日任務完成2012年1月 6日考核方式1.評閱 2.答辯

2、 3.實際操作 4.其它 指導教師汪傳忠系（部）主任王長坤注：1、此表一組一表二份，課程設計小組組長一份；任課教師授課時自帶一份備查。2、課程設計結束后與“課程設計小結”、“學生成績單”一并交院教務存檔。目錄摘要3一、課程設計目的及要求4二、課程設計內容及任務安排4 三、設計原理及設計方法53.1 FIR數字濾波器5 3.2 IIR數字濾波器7四、與設計相關的知識8五、設計過程14 5.1圖像加噪處理及功率密度譜對比14 5.2 FIR濾波器的設計22 (1)布萊克曼窗.22 (2)海明窗26 （3）漢寧窗.28 （4）多爾夫-切比雪夫窗.32 （5）巴特利特窗.35 5.3 IIR 數字濾波

3、器設計.38 （1）巴特沃茲.38 （2）切比雪夫1.41 （3）切比雪夫2 .43 （4）橢圓濾波.45六、FIR和IIR的比較.47七、個人設計總結49附錄50摘要MATLAB軟件具有使用簡單、方便，易編程，語言簡練，函數庫可任意擴充，采用全新數據類型和面向對象編程技術等特點，具有數值分析容易、語法結構簡單、矩陣數值運算高效、圖形功能完備、圖像處理方便等強大的功能，因此已被廣泛應用在教學、科研和工程設計的各個領域。傳統(tǒng)的濾波器分析與設計均使用繁瑣的公式計算，改變參數后需要重新運算，從而在分析與設計濾波器尤其是高階濾波器時工作量特別大，利用MATLAB信號處理工具箱可以快速而有效地實現濾波器

4、的分析、設計與仿真。本課設介紹基于MATLAB環(huán)境的IIR數字濾波器和FIR數字濾波器的設計方法，原理和步驟。IIR數字濾波器的常用設計方法有脈沖響應不變法和雙線性變換法等，FIR數字濾波器的設計方法只要有窗函數設計法，頻率抽樣設計法和最優(yōu)化設計法等，最后利用設計的濾波器對含噪聲的圖像進行濾波處理并驗證所設計濾波器的性能，同時比較兩種設計方法的最優(yōu)性。關鍵字：數字信號處理 Matlab FIR濾波器 IIR濾波器 圖像加噪 圖像濾波 FIR與IIR濾波器比較 1、 課程設計目的及要求1.1目的 綜合運用數字信號處理的理論知識進行頻譜分析和濾波器設計FIR和IIR濾波器，通過理論推導得出相應結論

5、，再利用 MATLAB 作為編程工具進行計算機實現，從而加深對所學知識的理解，建立概念。 1、鞏固和加深理解數字信號處理基本理論知識，提高綜合運用這些知識進行程序設計的能力。2、能根據需要查閱相關手冊和資料，提高獨立分析問題和解決問題的能力。3、能正確使用儀器設備，掌握測試原理，掌握綜合性程序設計的基本過程。4、掌握時域離散時間信號和時域離散系統(tǒng)基本概念和基本性質。5、掌握時域離散信號、系統(tǒng)的頻域分析工具。6、掌握基于Matlab的數字信號處理領域的計算機輔助分析與設計的基本方法。7、培養(yǎng)團隊協(xié)作精神、創(chuàng)新意識、嚴肅認真的治學態(tài)度和嚴謹求實的工作作風1.2基本要求 1、 學會 MATLAB 的

6、使用，掌握 MATLAB 的程序設計方法。 2、掌握在 Windows 環(huán)境下數字圖像處理的方法。 3、 掌握數字信號處理的基本概念、基本理論和基本方法。 4、 掌握 MATLAB 設計 FIR 和 IIR 數字濾波器的方法。 5、 學會用 MATLAB 對數字圖像進行分析和處理二、課程設計內容及任務安排 每人用自己的一張照片，并對照片進行加噪聲處理，通過設計濾波器濾除噪聲。每組三人，其中一人設計FIR濾波器進行噪聲處理，一人設計IIR濾波器濾除噪聲，另外一人分析FIR和IIR濾波器的優(yōu)劣比較。其中在進行濾波器設計中，要求每種模型都要求設計，FIR濾波器設計包括布萊克曼窗、海明窗、漢寧窗、多爾

7、夫切比雪夫窗、巴特利特窗，IIR濾波器設計包括巴特沃茲、切比雪夫1和2，橢圓模型。并且分析每個模型的幅頻特性，相頻特性，相延時，群延時。最后得出結論。具體實施方案：1、蘇青文主要負責設計IIR無限脈沖響應數字濾波器。2、方雪松主要負責設計FIR有限脈沖響應數字濾波器。3、組長陳舒龍主要負責圖像處理、整理報告及FIR和IIR兩種數字濾波器的比較。三、設計原理及設計方法數字濾波器設計的一個重要步驟就是確定一個可實現的傳遞函數G(z)來逼近指定的頻率響應。在G(z)確定下來之后，接下來要設計用一種合適的濾波器結構來實現它。濾波器設計過程的第一步是估計傳輸函數階數。然后根據各種模型用Matlab函數得

8、出傳輸函數。3.1 FIR數字濾波器3.1.1 基本介紹FIR是有限沖激響應（Finite Impulse Response）的簡稱。 由線性系統(tǒng)理論可知，在某種適度條件下，輸入到線性系統(tǒng)的一個沖擊完全可以表征系統(tǒng)。當我們處理有限的離散數據時，線形系統(tǒng)的響應（包括對沖擊的響應）也是有限的。若線性系統(tǒng)僅是一個空間濾波器，則通過簡單地觀察它對沖擊的響應，我們就可以完全確定該濾波器。通過這種方式確定的濾波器稱為有限沖擊響應（FIR)濾波器. 3.1.2 設計方法 FIR濾波器設計的最簡單方法是，對期望的理想頻率響應進行離散時間傅里葉逆變換，得到雙無限長沖擊響應，對所得的無限系數簡單地截尾為有限系數。

9、濾波器類型確定之后，接下來是估計濾波器階數N，matlab里提供了階數估計的函數，我們可以直接用，階數估計之后就是確定傳輸函數num和den，每種模型都有這樣的函數求出傳輸函數。本設計用了5個模型加窗函數：布萊克曼窗、海明窗、漢寧窗、多爾夫切比雪夫窗、巴特利特窗，來設計FIR濾波器。幾種窗函數的性能比較窗函數 主瓣寬度ML相對旁瓣級別（dB）最小阻帶衰減（dB）過渡帶帶寬矩形窗 4/(2M+1)13.3 20.9 0.92/M巴特利特窗 4/(M+1)26.5漢寧窗 8/(2M+1)31.543.93.11/M海明窗8/(2M+1)42.754.5 3.32/M布拉克曼窗 12/(2M+1)5

10、8.175.3 5.56/M3.2 IIR數字濾波器3.2.1 IIR數字濾波器的特點 a.封閉函數IIR數字濾波器的系統(tǒng)函數可以寫成封閉函數的形式。 b.IIR數字濾波器采用遞歸型結構IIR數字濾波器采用遞歸型結構，即結構上帶有反饋環(huán)路。IIR濾波器運算結構通常由延時、乘以系數和相加等基本運算組成，可以組合成直接型、正準型、級聯(lián)型、并聯(lián)型四種結構形式，都具有反饋回路。由于運算中的舍入處理，使誤差不斷累積，有時會產生微弱的寄生振蕩。 c.借助成熟的模擬濾波器的成果IIR數字濾波器在設計上可以借助成熟的模擬濾波器的成果，如巴特沃斯、契比雪夫和橢圓濾波器等，有現成的設計數據或圖表可查，其設計工作量

11、比較小，對計算工具的要求不高。在設計一個IIR數字濾波器時，我們根據指標先寫出模擬濾波器的公式，然后通過一定的變換，將模擬濾波器的公式轉換成數字濾波器的公式。 d.需加相位校準網絡IIR數字濾波器的相位特性不好控制，對相位要求較高時，需加相位校準網絡3.2.2 IIR數字濾波器的設計方法IIR數字濾波器的最通用的方法是借助于模擬濾波器的設計方法，基于原型模擬傳輸函數的雙線性變換。模擬傳輸函數通常有四種模型：巴特沃茲、切比雪夫1型和2型，橢圓模型。模擬濾波器設計已經有了相當成熟的技術和方法，有完整的設計公式，還有比較完整的圖表可以查詢，因此設計數字濾波器可以充分利用這些豐富的資源來進行。 對于I

12、IR數字濾波器的設計具體步驟如下： (1)按照一定的規(guī)則將給出的數字濾波器的技術指標轉換為模擬低通濾波器的技術指標。 (2)根據轉換后的技術指標設計模擬低通濾波器G(s)(G(s)是低通濾波器的傳遞函數)。 (3)再按照一定的規(guī)則將G(s)轉換成H(z)(H(z)是數字濾波器的傳遞函數)。若設計的數字濾波器是低通的，上述的過程可以結束，若設計的是高通、帶通或者是帶阻濾波器，那么還需要下面的步驟：將高通、帶通或帶阻數字濾波器的技術指標轉換為低通模擬濾波器的技術指標，然后設計出低通G(s)，再將G(s)轉換為H(z)。 本設計應用matlab進行設計，對模型直接應用Matlab提供的函數來估計階數

13、，然后確定傳輸函數。本設計設計了四種模型：巴特沃茲、切比雪夫1型和2型，橢圓模型。四、與設計相關的背景知識1、功率譜密度 功率譜密度（PSD），它定義了信號或者時間序列的功率如何隨頻率分布。這里功率可能是實際物理上的功率，或者更經常便于表示抽象的信號被定義為信號數值的平方，也就是當信號的負載為1歐姆(ohm)時的實際功率。由于平均值不為零的信號不是平方可積的，所以在這種情況下就沒有傅里葉變換。維納-辛欽定理（Wiener-Khinchin theorem）提供了一個簡單的替換方法，如果信號可以看作是平穩(wěn)隨機過程，那么功率譜密度就是信號自相關函數的傅里葉變換。 信號的功率譜密度當且僅當信號是廣義

14、的平穩(wěn)過程的時候才存在。如果信號不是平穩(wěn)過程，那么自相關函數一定是兩個變量的函數，這樣就不存在功率譜密度，但是可以使用類似的技術估計時變譜密度。功率譜的概念是針對功率有限信號的，所表現的是單位頻帶內信號功率隨頻率的變化情況。保留了頻譜的幅度信息，但是丟掉了相位信息，所以頻譜不同的信號其功率譜是可能相同的。 有兩點需要注意：（1）. 功率譜是隨機過程的統(tǒng)計平均概念，平穩(wěn)隨機過程的功率譜是一個確定函數；而頻譜是隨機過程樣本的Fourier變換，對于一個隨機過程而言，頻譜也是一個“隨機過程”。（隨機的頻域序列）（2）. 功率概念和幅度概念的差別。此外，只能對寬平穩(wěn)的各態(tài)歷經的二階矩過程談功率譜，其存

15、在性取決于二階矩是否存在并且二階矩的Fourier變換收斂；而頻譜的存在性僅僅取決于該隨機過程的該樣本的Fourier變換是否收斂。功率譜密度譜是一種概率統(tǒng)計方法，是對隨機變量均方值的量度。一般用于隨機振動分析，連續(xù)瞬態(tài)響應只能通過概率分布函數進行描述，即出現某水平響應所對應的概率。 功率譜密度的定義是單位頻帶內的“功率”（均方值） 功率譜密度是結構在隨機動態(tài)載荷激勵下響應的統(tǒng)計結果，是一條功率譜密度值頻率值的關系曲線，其中功率譜密度可以是位移功率譜密度、速度功率譜密度、加速度功率譜密度、力功率譜密度等形式。數學上，功率譜密度值頻率值的關系曲線下的面積就是方差，即響應標準偏差的平方值。2、ui

16、getfile命令功能：創(chuàng)建標準的對話框并通過交互式操作取得文件名uigetfile：顯示一個模態(tài)對話框，對話框列出了當前目錄下的文件和目錄，用于可以選擇一個將要打開的文件名。如果文件名是有效的且該文件存在，則當用戶點擊Open時函數uigetfile返回該文件名。若不存在，uigetfile顯示一個控制返回對話框值的錯誤提示信息，此時用戶可以輸入另外的文件名或點擊Cancel按鈕。如果用戶點擊Cancel按鈕或關閉對話框，函數uigetfile將返回0。FileName,PathName,FilterIndex = uigetfile(FilterSpec)：創(chuàng)建并顯示對話框,列出當前目錄下

17、的文件。FilterSpec決定文件的初始顯示，它可以為一個文件全名或者包含通配符*，如*.doc表示列出所有的文檔文件。如果FilterSpec是一個包含文件名的字符串，則該文件名所在的文件名域將被選中并顯示出來，且該文件的擴展名將作為過濾作用；如果FilterSpec是一個包含文件路徑的字符串，例如.,.或/,例如./*.m表示列出當前目錄上的所有M文件；如果FilterSpec是一個單元數組，則其第一列為擴展名列表，第二列為描述列表；如果FilterSpec沒有指定路徑對象，則uigetfile將使用默認路徑下的所有文件類型。此函數的用法為FileName,PathName,Filter

18、Index= uigetfile(FilterSpec,DialogTitle,DefaultName)FileName：返回的文件名PathName：返回的文件的路徑名FilterIndex：選擇的文件類型FilterSpec：文件類型設置DialogTitle：打開對話框的標題DefaultName：默認指向的文件名3、Imread 讀取圖像文件 函數功能：該函數用于讀取圖片文件中的數據。在matlab的命令窗口中輸入doc imread或者help imread即可獲得該函數的幫助信息。 調用格式：A = imread(filename,fmt) X,map = imread(filen

19、ame,fmt) . = imread(filename) . = imread(URL,.) . = imread(.,idx) (CUR, GIF, ICO, and TIFF only) . = imread(.,'PixelRegion',ROWS, COLS) (TIFF only) . = imread(.,'frames',idx) (GIF only) . = imread(.,ref) (HDF only) . = imread(.,'BackgroundColor',BG) (PNG only) A,map,alpha = i

20、mread(.) (ICO, CUR, and PNG only4、Imwrite 圖像數據保存 函數功能：將圖像數據寫入到圖像文件中， 存儲在磁盤上。在matlab命令窗口中鍵入doc imwrite或help imwrite可以獲得更多關于該函數的幫助信息。 調用格式： imwrite(A,filename,fmt) A是圖像數據， filename是目標圖像名字， fmt是要生成的圖片的格式。圖像格式有：bmp（1-bit、8-bit和24-bit）、gif（8-bit）、hdf、jpg（或jpeg）（8-bit、12-bit和16-bit）、jp2或jpx、pbm、pcx（8-bit）

21、、pgm、png、pnm、ppm、ras、tif（或tiff）、xwd。各種格式支持的圖像位數不一樣， 比如bmp格式不支持16-bit，而png格式支持， 又如gif只支持8-bit格式。5、Imnoise 加噪聲處理函數 函數imnoise來使用噪聲污染一幅圖像，該函數的基本語法為： g=imnoise(f,type,parameters) f為是輸入圖像。函數imnoise在給圖像添加噪聲之前，將它轉換為范圍0,1內的double類圖像。指定噪聲參數時必須考慮到這一點。 g=imnoise(f,'gaussian',m,var)將均值M,方差為var的高斯噪聲加到圖像f上

22、，默認值為均值是0，方差是0.01的噪聲。 g=imnoise（f,'localvar',V)將均值為0，局部方差為V的高斯噪聲添加到圖像f上，其中V是與f大小相同的一個數組，它包含了每一個點的理想方差值。 g=imnoise(f,'localvar',image_intensity,var)將均值為0的高斯噪聲添加到圖像f中，其中噪聲的局部方差var是圖像f的亮度值的函數。參量image_intensity和var是大小相同的向量，plot(image_intensity,var)繪制出噪聲方差和圖像亮度的函數關系。向量image_intensity必須包含范

23、圍在0,1內的歸一化亮度值。 g=imnoise(f,'salt&pepper',d)用椒鹽噪聲污染圖像f，其中d是噪聲密度（即包括噪聲值的圖像區(qū)域的百分比）。因此，大約有d*numel（f）個像素受到影響。默認的噪聲密度為0.05。 g=imnoise(f,'speckle',var)用方程g=f+n*f)將乘性噪聲添加到圖像f上，其中n是均值為0，方差為var的均勻分布的隨機噪聲，var的默認值是0.04。 g=imnoise(f,'poisson')從數據中生成泊松噪聲，而不是將人工的噪聲添加到數據中，為了遵守泊松統(tǒng)計，unit8和

24、unit16類圖像的亮度必須和光子的數量相符合。當每個像素的光子數量大于65535時，就要使用雙精度圖像。亮度值在0到1之間變化，并且對應于光子的數量除以10e12。五、設計過程5.1 圖像加噪處理及功率密度譜對比(1)陳舒龍 加噪前加噪后 (2)蘇青文 加噪前加噪后(3)方雪松 加噪前加噪后5.2 FIR濾波器1，布萊克曼'blackman'經多次調試下面參數濾波效果最好。Fp=2200; Fs=3600; Rp=0.6; Rs=15; FT=20000; num,den=FIR(Fp,Fs,Rp,Rs,FT,'blackman'); 2，海明窗'ha

25、mming'經多次調試下面參數濾波效果最好。Fp=2400;Fs=3200;Rp=0.6; Rs=15; FT=20000; num,den=FIR(Fp,Fs,Rp,Rs,FT,'hamming'); 3，漢寧窗 經多次調試下面參數濾波效果最好。 Fp=2400; Fs=2800; Rp=0.6; Rs=15; FT=20000; num,den=FIR(Fp,Fs,Rp,Rs,FT,'hanning'); 4，多爾夫切比雪夫窗 經多次調試下面參數濾波效果最好。Fp=2300;Fs=3100;Rp=0.6; Rs=15; FT=20000; num,d

26、en=FIR(Fp,Fs,Rp,Rs,FT,'cheby'); 5，巴特利特窗'battlet' 經多次調試下面參數濾波效果最好。 Fp=1960; Fs=5200; Rp=0.6; Rs=15; FT=20000; num,den=FIR(Fp,Fs,Rp,Rs,FT,'battlet'); 5.3 IIR數字濾波器設計關鍵是確定Wp，Ws，Rp，Rs參數用 for nn=1:1 Fp=1940+nn*100;Fs=5300+nn*500;循環(huán)5次，Fs每次先500甚至更大增大，觀察圖像濾波效果，選擇兩張圖像較好的進行大定位，然后以在兩個值之間

27、進行小范圍循環(huán)5次，觀察濾波效果，多次嘗試，最終選擇一個最好的參數。1、巴特沃茲濾波Wp=0.5; Ws=0.6;Rp=1.5; Rs=36;此參數效果最好。 num den=iirlp(Wp,Ws,'butter'); 2. 切比雪夫1 Wp=0.3; Ws=0.6;Rp=1.5; Rs=38;此參數效果最好。 num den=iirlp(Wp,Ws,' cheby1');3.切比雪夫2型Wp=0.5; Ws=0.6;Rp=1.5; Rs=38;此參數效果最好。num den=iirlp(Wp,Ws,' cheby2');4、橢圓濾波Wp=0.

28、5; Ws=0.6;Rp=1.5; Rs=8;此參數效果最好num den=iirlp(Wp,Ws,'ellip');比較結果分析利用巴特沃茲、切比雪夫1、切比雪夫2、橢圓濾波器都可以進行低通、高通、帶通、帶阻濾波器的設計，但是各有特點。巴特沃茲濾波器通帶內的幅頻響應曲線能得到最大限度的平滑，但犧牲了截止頻率的坡度；切比雪夫1型濾波器通帶內等波紋，阻帶內單調；切比雪夫2型濾波器通帶內單調，阻帶內等波紋；橢圓濾波器阻帶內和通帶內都是等波紋的，但下降的坡度更大，而且可以以更低的階數實現和其他兩類濾波器一樣的性能指標。六、FIR和IIR比較a.從性能上進行比較 從性能上來說，IIR濾

29、波器傳輸函數的極點可位于單位圓內的任何地方，因此可用較低的階數獲得高的選擇性，所用的存貯單元少，所以經濟而效率高。但是這個高效率是以相位的非線性為代價的。選擇性越好，則相位非線性越嚴重。相反，FIR濾波器卻可以得到嚴格的線性相位，然而由于FIR濾波器傳輸函數的極點固定在原點，所以只能用較高的階數達到高的選擇性；對于同樣的濾波器設計指標，FIR濾波器所要求的階數可以比IIR濾波器高510倍，結果，成本較高，信號延時也較大；如果按相同的選擇性和相同的線性要求來說，則IIR濾波器就必須加全通網絡進行相位較正，同樣要大增加濾波器的節(jié)數和復雜性。b.從結構上看IIR濾波器必須采用遞歸結構，極點位置必須在

30、單位圓內，否則系統(tǒng)將不穩(wěn)定。另外，在這種結構中，由于運算過程中對序列的舍入處理，這種有限字長效應有時會引入寄生振蕩。相反，FIR濾波器主要采用非遞歸結構，不論在理論上還是在實際的有限精度運算中都不存在穩(wěn)定性問題，運算誤差也較小。此外，FIR濾波器可以采用快速付里葉變換算法，在相同階數的條件下，運算速度可以快得多。兩者的優(yōu)缺點如下：FIR有限沖激響應濾波器原理：無反饋（無遞歸），極點在原點，卷積（或FFT卷積，部分也可用遞歸計算）優(yōu)點：線性相位，頻譜特性好缺點：速度較慢（由于階數多延時長）應用場合：需線性相位的濾波（如圖像信號等），線性相位是指不同頻率的信號經過濾波器其延時一樣。FIR的應用靈活

31、，易于適應某些特殊的應用。IIR無限沖激響應濾波器原理：有反饋（有遞歸），極點可位于單位圓內任何地方，遞歸計算優(yōu)點：階數少（延時少），借助成熟的模擬濾波器的成果（如巴特沃斯、契比雪夫和橢圓濾波器等）缺點：非線性相位，易出現震蕩（不穩(wěn)定，因為遞歸計算）應用場合：不需線性相位的濾波（如音頻信號等）總之，可以將IIR濾波器和FIR濾波器的性能指標以及優(yōu)缺點歸納如下：IIR 數字濾波器FIR 數字濾波器（1）相位一般是非線性的（1）相位可以做到嚴格線性（2）不一定穩(wěn)定（2）一定是穩(wěn)定的（3）對頻率分量的選擇性好（零極點可同時起作用）（3）選擇性差（4）相同性能下階次較低（4）相同性能下階次高（5）有噪

32、聲反饋，噪聲大（5）噪聲?。?）設計有封閉形式的公式，一次完成（6）沒有封閉形式的設計公式，須靠經驗與反復調試（7）對計算手段的要求較低（7）一般需用計算機計算七、個人設計總結1.陳舒龍本次課設我作為組長給咱組員分配任務，不過由于蘇青文同學家里有事提前回家了，他之前負責做的IIR濾波器設計不太準確，所以我就做了較大修改。通過對巴特沃茲、切比雪夫1型和2型、橢圓型濾波器的設計，從圖形直觀分析驗證理論知識。比較這四種濾波器發(fā)現：巴特沃茲濾波器有最寬的過渡帶，有單調遞減的增益響應。兩種類型的切比雪夫濾波器有相等寬度的過渡帶，該寬度小于巴特沃茲但大于橢圓。切比雪夫1型在過渡帶中有一個比切比雪夫2型稍微

33、快一些的下降。在通帶中，切比雪夫2型的幅度響應與巴特沃茲的幾乎一樣。橢圓有最窄的過渡帶并具有等波紋通帶和等波紋阻帶響應。 另外我的本職是分析比較IIR和FIR，其中，FIR濾波器的最大好處是穩(wěn)定、線性相位和廣泛的適用范圍，而它的最大缺點是階數高，從而帶來時延大、存儲單元多等問題。但IIR與FIR濾波器各有所長，所以在實際應用時應該從多方面考慮來加以選擇。從使用要求上來看，在對相位要求不敏感的場合，如語言通信等，選用IIR較為合適，這樣可以充分發(fā)揮其經濟高效的特點；對于圖像信號處理，數據傳輸等以波形攜帶信息的系統(tǒng)，則對線性相位要求較高。如果有條件，采用FIR濾波器較好。當然，在實際應用中可能還要

34、考慮更多方面的因素。不論IIR和FIR，階數越高，信號延遲越大；同時在IIR濾波器中，階數越高，系數的精度要求越高，否則很容易造成有限字長的誤差使極點移到單位圓外。因此在階數選擇上是綜合考慮的。2.方雪松通過本次的課程設計，讓我平時在課本上學習的知識得以運用。在這次課設中，我主要任務是負責FIR濾波器的設計，雖然設計過程中遇到許多問題，但通過與同學的交流以及在網上查閱相關資料都逐漸得以解答。通過對FIR濾波器的設計讓我認識到它可進行準確的線性相位設計且結構具有穩(wěn)定的量化濾波器系數，FIR濾波器可以得到嚴格的線性相位，然而由于FIR濾波器傳輸函數的極點固定在原點，所以只能用較高的階數達到高的選擇

35、性，成本較高，信號延時也較大。經過本次課程設計后，讓我對以前在課堂上學習的知識有了更深層次的理解以及FIR濾波器的設計和基本應用，最后和其他兩位組員合作共同完成了本次課程設計。附錄：%圖像濾波function DE_NOISE_IMAGEfilename, pathname = uigetfile( .'*.jpg','MATLAB Files (*.jpg)' '*.m', 'M-files (*.m)' . '*.fig','Figures (*.fig)' . '*.mat',

36、'MAT-files (*.mat)' . '*.mdl','Models (*.mdl)' . '*.*', 'All Files (*.*)', . 'Pick a file');q,qq=size(pathname);S=''q1,qq1=size(filename);RNAME=pathname;RNAME(qq+1)=S(1);RNAME(qq+2:qq+1+qq1)=filename(1:qq1);sh=imread(RNAME,'jpg');%選被噪聲污

37、染的圖象；M,N,O=size(sh);%sh=double(sh)/255;R=rand*0.025;if R<0.015 R=0.015;enddh1=sh;for o=1:3 dh1(:,:,o)=imnoise(dh1(:,:,o),'gaussian',0.015);enddh=dh1;dhL=dh1;dhI=dh1;%參數設計% Wp=0.5; Ws=0.6; Rp=1; Rs=40; %'iir_type butter','cheby1','cheby2','ellip'.%fir-type &

38、#39;blackman'',''hamming'',''hanning'',''cheby'',''kaiser''%num den=firlp(Wp,Ws,'blackman');num den=iirlp(Wp,Ws,'ellip');%從此處選擇濾波器類型for o=1:3 A(1:M,1:N)=dh(1:M,1:N,o); AL=A; AI=A; %用設計的濾波器對圖象濾波 %按列處理； for k=1:N x

39、=AL(:,k); AL(:,k)=filter(num,den,x);%設計出的IIR濾波器的參數num,den end dhL(1:M,1:N,o)=AL(1:M,1:N); %行列終和 %DH(1:M,1:N,o)=(dhI(1:M,1:N,o)+dhL(1:M,1:N,o)/2; DH(1:M,1:N,o)=dhL(1:M,1:N,o); %有約束濾波 endfigure;%圖片subplot(1,3,1);imshow(sh);title('原始圖像');subplot(1,3,2);imshow(dh1);title('加噪圖像');subplot(

40、1,3,3);imshow(DH);title('去噪圖像');figure;%功率譜subplot(3,1,1);PGJ_INPUT(sh);grid on;axis(0,1,-50,40);title('原始圖像功率譜');subplot(3,1,2);PGJ_INPUT(dh1);grid on;axis(0,1,-50,40);title('加噪功率譜');subplot(3,1,3);PGJ_INPUT(DH);grid on;%axis(0,1,-100,60);title('去噪功率譜');end%FIR濾波器設計f

41、unction num den=firlp(Wp,Ws,windows)%N num den=firlp(windows,ftype)%函數名：FIR低通濾波器的設計.%參數：窗函數windows,濾波器響應類型 ftype.%返回：濾波器的階數N和系數num den.%濾波器指標 Rp=0.1;Rs=25; dp=1-10(-Rp/20);ds=10(-Rs/20); %估計階數N Wn,beta,ftype=kaiserord(Wp Ws,1 0,dp ds);N=N+18;switch lower(windows) case 'blackman' taper=blackm

42、an(N+1); case 'hamming' N=N-10; taper=hamming(N+1); case 'hanning' N=N-13; taper=hanning(N+1); case 'cheby' N=N+6; taper=chebwin(N+1); case 'kaiser' N=N-15; taper=kaiser(N+1,beta); case 'battlet' N=N+120; taper=bartlett(N+1); otherwise error('Unknown windo

43、ws of ''%s'' ,please select one from ''blackman'',''hamming'',''hanning'',''cheby'',''kaiser''.',windows)endnum=fir1(N,Wn,ftype,taper);den=1; freqRespone1(num,den,windows,Wp,Ws);%頻率響應(subf)%function f

44、reqRespone(num,den,FT,methods,Wp,Ws)function freqRespone(num,den,methods,Wp,Ws)h wa=freqz(num,den,512);%濾波器頻率響應figure,subplot(2,2,1);plot(wa,(abs(h);gridylabel('振幅(dB)');xlabel('頻率(Hz)');title('頻率響應','(IIR-',methods,'逼近','Wp=',num2str(Wp),' Hz,'

45、;,'Ws=',num2str(Ws),' Hz.)');subplot(2,2,2);plot(wa,unwrap(angle(h)*180/pi);gridxlabel('頻率(Hz)');ylabel('相位(degree)');subplot(2,2,3);plot(wa,grpdelay(h);xlabel('頻率(Hz)');ylabel('群延時(grpdelay)');subplot(2,2,4);plot(wa,phasedelay(h);xlabel('頻率(Hz)

46、9;);ylabel('相延時(phasedelay)')figure,zplane(num,den);title('零極點圖','(IIR-',methods,'逼近)');%頻率響應%function freqRespone1(num,den,windows,Wp,Ws)h,wa=freqz(num,den,512);%濾波器頻率響應figure;subplot(2,2,1);xlabel('階數(N)'),ylabel('h(N)');title('沖激響應','(N=

47、',. num2str(length(num),',FIR-',windows,'窗','Wp=',num2str(Wp),' Hz,','Ws=',num2str(Ws),' Hz.)');plot(wa,abs(h);gridxlabel('頻率(Hz)');ylabel('振幅(dB)');subplot(2,2,2);plot(wa,unwrap(angle(h)*180/pi);gridxlabel('頻率(Hz)');ylabel(

48、'相位(degree)');subplot(2,2,3);plot(wa,grpdelay(h);xlabel('頻率(Hz)');ylabel('群延時(grpdelay)');subplot(2,2,4);plot(wa,phasedelay(h);xlabel('頻率(Hz)');ylabel('相延時(phasedelay)')figure,zplane(num,den);%title('零極點圖'(FIR-',windows,'窗')%雙向性變換法設計IIR濾波器%

49、function num den=iirlp(Wp,Ws,methods)%N num den=iirlp(methods,ftype)%函數名：IIR低通濾波器的設計.%參數：濾波器的逼近類型 methods：'butter','cheby1','cheby2','ellip'.%返回：濾波器的階數N和傳輸函數num den.ftype='low'%濾波器設計%濾波器指標Rp=1;Rs=20;%設計switch lower(methods) case 'butter'%巴特沃斯逼近 N,Wn=buttord(Wp,Ws,Rp,Rs); num den=butter(N,Wn,ftype); freqRespone(num,den,methods,Wp,Ws); case 'cheby1'%切比