基于頻率抽樣法二維非遞歸濾波器設(shè)計(jì)和簡(jiǎn)單應(yīng)

1、 . . . 一、 開展本課題的意義與工作容：數(shù)字濾波器被廣泛應(yīng)用于語音、圖像、無線電等領(lǐng)域，具有廣闊的發(fā)展空間，特別是在圖像處理中的應(yīng)用非常廣泛。隨著數(shù)字信號(hào)處理系統(tǒng)的大量運(yùn)用與數(shù)字信號(hào)處理技術(shù)的不斷發(fā)展，人們對(duì)數(shù)字濾波器的要求也越來越多樣化。本設(shè)計(jì)主要根據(jù)FIR濾波器性能要求，運(yùn)用最小二乘法設(shè)計(jì)一個(gè)一維優(yōu)化的FIR濾波器，然后通過一定的映射關(guān)系將一維濾波器轉(zhuǎn)換成二維數(shù)字濾波器，實(shí)現(xiàn)對(duì)數(shù)字圖像的濾波，并分析濾波器性能。二、 課題工作的總體安排與進(jìn)度：總體安排：2007-2008學(xué)年第一學(xué)期19-21周確定畢業(yè)設(shè)計(jì)課題；20072008學(xué)年第二學(xué)期1-7周在教師指導(dǎo)下進(jìn)行做畢業(yè)論文（設(shè)計(jì)），第

2、8周，交畢業(yè)設(shè)計(jì)（論文）成果，院組織專門人員根據(jù)規(guī)要求對(duì)畢業(yè)設(shè)計(jì)（論文）進(jìn)行形式審查。第9周，指導(dǎo)教師評(píng)定本科畢業(yè)生的畢業(yè)論文（設(shè)計(jì)），畢業(yè)論文交叉評(píng)閱，學(xué)生答辯資格審定與組織學(xué)生答辯。進(jìn)度：2007-2008學(xué)年第一學(xué)期19-21周確定畢業(yè)設(shè)計(jì)課題，在寒假期間與2007-2008學(xué)年第二學(xué)期1-2周準(zhǔn)備好畢業(yè)設(shè)計(jì)的相關(guān)文獻(xiàn)資料。3-7周進(jìn)行畢業(yè)設(shè)計(jì)與論文的編寫。其中在第6周基本完成預(yù)期工作，論文初稿與仿真程序提交審批。04月22日對(duì)照審批意見對(duì)初稿進(jìn)行修改并在24好最終定稿。三、 課題預(yù)期達(dá)到的效果：設(shè)計(jì)基本達(dá)到預(yù)期的成果，成功設(shè)計(jì)了一個(gè)一維優(yōu)化的FIR濾波器，然后通過一定的映射關(guān)系將一維濾

3、波器轉(zhuǎn)換成二維數(shù)字濾波器，實(shí)現(xiàn)對(duì)數(shù)字圖像的濾波?；陬l率抽樣法二維非遞歸濾波器設(shè)計(jì)與其簡(jiǎn)單應(yīng)用摘要隨著數(shù)字信號(hào)處理系統(tǒng)的大量運(yùn)用與數(shù)字信號(hào)處理技術(shù)的不斷發(fā)展，人們對(duì)數(shù)字濾波器的要求也越來越多樣化。數(shù)字濾波器在數(shù)字信號(hào)處理的各種應(yīng)用中發(fā)揮著十分重要的作用，它是通過對(duì)采樣數(shù)據(jù)信號(hào)進(jìn)行數(shù)學(xué)運(yùn)算處理來達(dá)到頻域?yàn)V波的目的。數(shù)字濾波器是提取有用信息非常重要、非常靈活的方法，是現(xiàn)代信號(hào)處理的重要容，它被廣泛應(yīng)用于語音、圖像、無線電、醫(yī)療器件等領(lǐng)域，具有廣闊的發(fā)展空間，特別是在醫(yī)療圖像處理中的應(yīng)用非常廣泛。本設(shè)計(jì)利用最小二乘優(yōu)化算法設(shè)計(jì)一個(gè)一維FIR數(shù)字濾波器，再根據(jù)頻率抽樣法設(shè)計(jì)二維FIR濾波器，并在Mat

4、lab中進(jìn)行仿真。利用所設(shè)計(jì)出的濾波器對(duì)一幅X光兒童手腕骨醫(yī)學(xué)圖像濾波，為X光圖像的后續(xù)處理作準(zhǔn)備。1 前言1.1 課題研究目的和意義11.2二維數(shù)字濾波器的研究現(xiàn)狀.11.3本課題研究的容22 FIR濾波器設(shè)計(jì)技術(shù)與波波特性32.1 傳統(tǒng)FIR濾波器設(shè)計(jì)方法32.2改進(jìn)的頻率抽樣法43 基于頻率抽樣法二維FIR濾波器的設(shè)計(jì)與圖像濾波實(shí)現(xiàn)93.1設(shè)計(jì)一維FIR濾波器93.2 最小二乘優(yōu)化設(shè)計(jì)103.3 基于頻率抽樣法二維FIR濾波器的程序?qū)崿F(xiàn)123.4基于頻率抽樣法圖像濾波的實(shí)現(xiàn)143.5結(jié)果與展望.163.5.1設(shè)計(jì)結(jié)果與所取得成績(jī)163.5.2 研究展望16參考文獻(xiàn)17致18附錄A 程序代

5、碼.191前言1.1課題研究的目的和意義根據(jù)結(jié)構(gòu)的不同，數(shù)字濾波器分為無限沖擊響應(yīng)（IIR）數(shù)字濾波器和有限沖擊響應(yīng)（FIR）數(shù)字濾波器。其中，由于FIR數(shù)字濾波器總是穩(wěn)定的，而且可以做到嚴(yán)格的線性相位，故而人們對(duì)它的設(shè)計(jì)問題進(jìn)行了大量而深入的研究。1FIR數(shù)字濾波器的設(shè)計(jì)指標(biāo)通常有最小二乘指標(biāo)、Chebyshev指標(biāo)、Lp指標(biāo)與混合指標(biāo)等。因?yàn)樽钚《酥笜?biāo)設(shè)計(jì)簡(jiǎn)單，Chebyshev指標(biāo)物理意義明確，所以此兩種指標(biāo)最為常用。70年代McClellan和Parks設(shè)計(jì)出基于Remez算法的計(jì)算機(jī)程序McClellan-Parks(MP)程序能高效的解決線性相位FIR數(shù)字濾波器數(shù)字濾波器的Che

6、byshev設(shè)計(jì)問題，這也使得Chebyshev誤差指標(biāo)在FIR濾波器的設(shè)計(jì)中被廣泛的采用。2數(shù)字濾波器被廣泛應(yīng)用于語音、圖像、無線電等領(lǐng)域，具有廣闊的發(fā)展空間，特別是在圖像處理中的應(yīng)用非常廣泛。隨著數(shù)字信號(hào)處理系統(tǒng)的大量運(yùn)用與數(shù)字信號(hào)處理技術(shù)的不斷發(fā)展，人們對(duì)數(shù)字濾波器的要求也越來越多樣化。比如在通訊系統(tǒng)中要求盡量減小阻帶能量，在圖像處理中要求保持信號(hào)的直流分量不變，在許多場(chǎng)合下都要求去掉50Hz的電源干擾信號(hào)，而為達(dá)到這些目的就需要在濾波器的設(shè)計(jì)中加入各種約束。1.2二維數(shù)字濾波器的研究現(xiàn)狀濾波器的種類很多，分類方法也不同，如可以從功能上分，也可以從實(shí)現(xiàn)方法上分，或從設(shè)計(jì)方法上來分等等。經(jīng)

7、典濾波器是假定輸入信號(hào)x(n)中的有效信號(hào)和噪聲（或干擾）信號(hào)成分各在不同的頻帶，當(dāng)x(n)通過一個(gè)線性濾波系統(tǒng)后，可以將噪聲信號(hào)成分有效地去除。可是，如果有效信號(hào)和噪聲信號(hào)的頻率帶相互重疊，那么經(jīng)典的濾波器將無能為力。現(xiàn)在的地質(zhì)雷達(dá)信號(hào)處理中的濾波器主要采用經(jīng)典的濾波器進(jìn)行處理。因此有時(shí)濾波效果較好，有時(shí)較差?，F(xiàn)代濾波理論研究的主要容是從含有噪聲的數(shù)據(jù)記錄（又稱為時(shí)間序列）中估計(jì)出信號(hào)的某些特征或信號(hào)本身。一旦信號(hào)被估計(jì)出，那么估計(jì)出的信號(hào)將比原信號(hào)會(huì)有高的信噪比?，F(xiàn)代濾波器把信號(hào)和噪聲都視為隨機(jī)信號(hào)，利用它們的統(tǒng)計(jì)特征（如自相關(guān)函數(shù)、功率譜函數(shù)等等）導(dǎo)出一套最佳的估值算法，然后用硬件和軟件

8、實(shí)現(xiàn)。目前現(xiàn)代濾波器主要有：維納濾波器、卡爾曼濾波器、線性預(yù)測(cè)器、自適應(yīng)濾波器等，很多專家將基于特征分解的頻率估計(jì)與奇異值分解算法都?xì)w入現(xiàn)代濾波器的疇。21.3本課題研究的容本設(shè)計(jì)主要根據(jù)FIR濾波器性能要求，運(yùn)用最小二乘法設(shè)計(jì)一個(gè)一維優(yōu)化的FIR濾波器，然后通過一定的映射關(guān)系將一維濾波器轉(zhuǎn)換成二維數(shù)字濾波器，實(shí)現(xiàn)對(duì)數(shù)字圖像的濾波，并分析濾波器性能。2 FIR濾波器設(shè)計(jì)技術(shù)與濾波特性2.1傳統(tǒng)FIR濾波器設(shè)計(jì)方法 FIR濾波器設(shè)計(jì)最常用的方法是窗函數(shù)設(shè)計(jì)法，其它常用的方法有頻率抽樣法，頻率變換法，切比雪夫逼近法等，以下對(duì)這些常用方法作簡(jiǎn)單介紹。窗函數(shù)法又稱傅里葉級(jí)數(shù)法，是最常用的設(shè)計(jì)方法，它以

9、hd(n)為媒介的時(shí)域設(shè)計(jì)法, 而濾波器指標(biāo)往往是在頻域給出的,為此,要由Hd( ej)算出hd(n) ,加窗后又從h(n)算出hd(n)來檢驗(yàn)。窗函數(shù)設(shè)計(jì)方法首先根據(jù)要求選擇一個(gè)適當(dāng)?shù)睦硐霝V波器,由于理想濾波器的脈沖響應(yīng)是非因果且無限長(zhǎng),用適當(dāng)?shù)拇昂瘮?shù)來截取它的脈沖響應(yīng),從而得到線性相位和因果的FIR濾波器。頻率抽樣法先對(duì)理想頻響Hd( ej)采樣,得到樣值H(k) ,再利用插值公式直接求出系統(tǒng)轉(zhuǎn)換函數(shù)H(z)；或者求出頻響H(ej),以便與理想頻響進(jìn)行比較。在0 ,2區(qū)間上對(duì)Hd(ej)進(jìn)行N 點(diǎn)采樣,等效于時(shí)域以N 為周期延拓。頻率抽樣法可以看做為插值法,它在采樣的i上保證等于；而在非插

10、值點(diǎn)(采樣點(diǎn)) 上,H(ej)是插值函數(shù)的線性疊加。為了提高逼近的質(zhì)量，減少逼近誤差，可以采用人為的擴(kuò)展過渡帶的方法，即在頻率相應(yīng)的過渡帶插入一個(gè)或多個(gè)比較連續(xù)的采樣點(diǎn)，使得過渡帶比較連續(xù)，從而使得通帶和阻帶之間變法比較緩慢，使得設(shè)計(jì)得到的濾波器對(duì)理想濾波器的逼近誤差較小3。切比雪夫方法從數(shù)值逼近的理論來看,對(duì)某個(gè)函數(shù)f(x)的逼近一般有三個(gè)方法10,11：插植法；最小平方逼近法；最佳一致逼近法。所謂插值, 即尋找一n 階多項(xiàng)式(或三角多項(xiàng)式) p(x) ,使它在n + 1 個(gè)點(diǎn)x0, x1 xn處滿足p(xk) = f(xk) , k= 0 ,1 , n而在非插值點(diǎn)上,p(x)是f(xk)

11、的某種組合。當(dāng)然,在非插植點(diǎn)上, p(x)和f(xk)存在一定的誤差。頻率抽樣法可以看作為插值法,它在抽樣點(diǎn)k上保證了= ,而在非抽樣點(diǎn)上, H(ej)是插值函數(shù)S (, k)的線性組合,其權(quán)重是。這種設(shè)計(jì)方法的缺點(diǎn)是通帶和阻帶的邊緣不容易精確的確定。這種設(shè)計(jì)方法是著眼于使整個(gè)區(qū)間 a , b的總誤差為最小,但它并不一定能保證在每個(gè)局部位置誤差都最小。在某些位置上,有可能存在著較大的誤差。實(shí)際上,傅立葉級(jí)數(shù)法就是一種最小平方逼近法。該方法在間斷點(diǎn)處出現(xiàn)了較大的過沖( Gibbs現(xiàn)象) 。為了減少這種過沖和欠沖,采用了加窗口的方法,當(dāng)然,加窗口以后的設(shè)計(jì)方法,已不在是最小平方逼近。2.2 改進(jìn)的

12、頻率抽樣法頻率抽樣法的特點(diǎn)是:(1)在采樣頻率上的逼近誤差為零，也就是理想和實(shí)際響應(yīng)的差為零；(2)其余頻率上的逼近誤差取決于理想響應(yīng)的形狀；理想響應(yīng)的輪廓越陡，則逼近誤差越大；(3)靠近帶的邊緣的誤差大，在帶的誤差小。頻率抽樣法目前有兩種設(shè)計(jì)方法，第一種直接用上面的基本思想，對(duì)逼近誤差不加任何限制；也就是說無論設(shè)計(jì)所得的誤差有多大我們都接受，這種方法叫樸素設(shè)計(jì)法。第二種方法則通過改變過渡帶的樣本值，努力使阻帶中的誤差極小化，以便產(chǎn)生一個(gè)較好的設(shè)計(jì)，這種方法叫最優(yōu)設(shè)計(jì)法。頻率抽樣法設(shè)計(jì)濾波器最大的優(yōu)點(diǎn)是直接從頻率域進(jìn)行設(shè)計(jì)，比較直觀，也適合于設(shè)計(jì)具有任意幅度特性的濾波器。缺點(diǎn)是邊緣頻率不易控制

13、。如果增加采樣點(diǎn)數(shù)N，對(duì)確定邊緣頻率有好處，但N加大會(huì)增加濾波器的成本。因此它適合于窄帶濾波器的設(shè)計(jì)。MATLAB信號(hào)處理工具箱提供了一個(gè)頻率抽樣法的設(shè)計(jì)函數(shù)fir2，它的典型調(diào)用方法為:h=fir2(M，f, A)。其中M是FIR濾波器的階數(shù)(濾波器的長(zhǎng)度為N=M+1)，長(zhǎng)度為N的數(shù)組h為濾波器系數(shù)(或脈沖響應(yīng))。數(shù)組f中包含各邊緣頻率，其單位為pi,0f1。f=1對(duì)應(yīng)于采樣頻率的一半，即奈奎斯特頻率。這些頻率必須以遞增順序排列，從0開始，到1結(jié)束。數(shù)組A為各指定頻率上預(yù)期的幅度響應(yīng)，f與A長(zhǎng)度必須相等，plot(f,A)應(yīng)該給出預(yù)期的濾波器幅頻特性。4現(xiàn)在利用頻率抽樣直接法設(shè)計(jì)一FIR濾波

14、器，低通濾波器技術(shù)指標(biāo)為p=0.2,Rp=0.25dBs=0.3,As=50dB選M=20，以使在p有一個(gè)頻率樣本，也即在k=2：p=0.2；下一個(gè)樣本在s，也即在k=3：s=0.3。這樣在通帶0p有3個(gè)樣本，在阻帶s有7個(gè)樣本。(2-1)由于M=20，并且這是一個(gè)II類線性相位濾波器，有(2-2)現(xiàn)在根據(jù) (2-3)將H(k)集合起來，并由確定脈沖響應(yīng)h(n)，用MATLAB實(shí)現(xiàn)，得到圖2.1。圖2.1 直接頻率抽樣設(shè)計(jì)法由圖2.1可見，最小阻帶衰減大約是16dB，這顯然是不可接受的。如果增加M，那么在過濾帶就一定有一些樣本，而對(duì)這些我們又并不完全知道頻率響應(yīng)。因此，在實(shí)際中很少采用直接設(shè)計(jì)

15、法。為了得到更大的衰減，就必須增大M，并讓過渡帶的樣本作為自由樣本；也就是說，改變它們的值以得到在已給定M下的最大衰減與其過渡帶寬。這個(gè)問題被認(rèn)為是一個(gè)優(yōu)化問題，可以用線性規(guī)劃技術(shù)來解決?，F(xiàn)在我們用最優(yōu)設(shè)計(jì)法來對(duì)上面的低通濾波器?，F(xiàn)選M=40以使有一個(gè)樣本在過濾帶。過濾帶的樣本是在k=5和在k=40-5=35。現(xiàn)用T1表示這些樣本值，且，那么已采樣的振幅響應(yīng)是 (2-4)由于，相位響應(yīng)的樣本是 (2-5)現(xiàn)在能夠變化T1以得到最好的最小阻帶衰減。這將會(huì)引起過渡帶加寬。首先看看當(dāng)T10.5時(shí)得到圖2.2。圖2.2最優(yōu)頻率設(shè)計(jì)法：T1=0.5由圖2.2可見，這個(gè)設(shè)計(jì)的最小阻帶衰減現(xiàn)在是30dB，這

16、就是比直接設(shè)計(jì)的衰減要好一些，但是仍然不在可以接受的50dB的水平上。通過人為地改變T1值可以示得最佳T1值，并且找到接近最優(yōu)解的值T1=0.39。由圖2.3可見，最優(yōu)阻帶衰減是43dB。明顯地要進(jìn)一步增加衰減就必須在過渡帶改變更多的樣本值。圖2.3 最優(yōu)頻率設(shè)計(jì)法：T1=0.39很清楚，通過變化一個(gè)樣本就能得到一種好得多的設(shè)計(jì)，從這點(diǎn)來看這種方法是很優(yōu)越的。實(shí)際上，過渡帶寬一般是很小的，就只包括一個(gè)或兩個(gè)樣本。因此需要優(yōu)化最多也就是兩個(gè)樣本以獲得最大的最小阻帶衰減。在絕對(duì)的意義上這也就是等效于使最大旁瓣幅度最小化，所以這類優(yōu)化問題也稱為最大最小化問題?，F(xiàn)在重新考慮上面的低通濾波器設(shè)計(jì)，用在過

17、渡帶的兩個(gè)樣本來解決它，以求得到一個(gè)更好的阻帶衰減。先M=60使得在過渡帶有兩上樣本。設(shè)這兩個(gè)過渡帶的樣本值是T1和T2，那么，給出為 (2-6)利用T1=0.5925和T2=0.1099這兩個(gè)值，可用MATLAB計(jì)算h(n)，得圖2.4。圖2.4低通濾波器設(shè)計(jì)圖由圖2.4可見，最小阻帶衰減現(xiàn)在是大約是63dB,這就是可以接受的了。53 基于頻率抽樣法二維FIR濾波器的設(shè)計(jì)與圖像濾波實(shí)現(xiàn)3.1 設(shè)計(jì)一維FIR濾波器按頻域采樣定理，F(xiàn)IR數(shù)字濾波器的傳輸函數(shù)H(z)和單位脈沖響應(yīng)h(n)可由它的N個(gè)頻域采樣值H(k)惟一確定，即 (3-1) (3-2)假定要設(shè)計(jì)的FIR數(shù)字濾波器頻率響應(yīng)為，直接

18、從頻域出發(fā)，對(duì)理想頻率響應(yīng)在02之間進(jìn)行N點(diǎn)的等間隔休樣，得到為 (3-3)取H(k)為理想濾波器在各頻率采樣點(diǎn)上的值，即令 (3-4)然后按頻域采樣定理，由H(k)求H(z)或h(n)。為了使設(shè)計(jì)的濾波器為線性相位濾波器，濾波器的頻率響應(yīng)與頻率采樣值必須滿中一定的條件，同時(shí)，在有些情況下，給定的理想濾波器的幅度響應(yīng)只給出了(0)轉(zhuǎn)頻率采樣的幅度值，因而上式中的k只能取0到N/2這（N/2+1）點(diǎn)，另外（N/2-1）點(diǎn)（從N/2+1到N-1）則需根據(jù)線性相位的奇偶對(duì)稱條件求得。對(duì)N為奇數(shù)，設(shè)計(jì)公式為 (3-5) (3-6) (3-7)對(duì)N為偶數(shù)，設(shè)計(jì)公式為 (3-8) (3-9) (3-10)

19、 (3-11)3.2最小二乘優(yōu)化設(shè)計(jì)單位脈沖響應(yīng)為h(n1,n2),n1=0,1,N1-1,n2=0,1,N2-1的二維FIR 濾波器,其頻率響應(yīng)為 (3-12)二維FIR 濾波器的約束最小二乘設(shè)計(jì)是指選擇N1×N2個(gè)脈沖響應(yīng)系數(shù)h(n1,n2) ,使H(ej1,ej2)逼近某個(gè)期望響應(yīng), 在逼近誤差小于預(yù)定值的條件下, 使平方逼近誤差的積分最小.一般情形下, N1×N2個(gè)脈沖響應(yīng)是獨(dú)立的。當(dāng)滿足一定對(duì)稱性時(shí), 獨(dú)立脈沖響應(yīng)系數(shù)有所減少, 且濾波器具有大多數(shù)應(yīng)用場(chǎng)合所要求的線性相位。為描述簡(jiǎn)單, 考慮N1, N2均為奇數(shù)、脈沖響應(yīng)矩形對(duì)稱的情形,即h(n1,n2)=h(N1

20、-1-n1,n2)=h(n1,N2-1-n2)=h(N1-1-n1,N2-1-n2)(3-13)此時(shí),共有數(shù)目為n=(N1+1)(N2+1)/4的獨(dú)立脈沖響應(yīng),即h(n1,n2),n1=0,1,r1,n2=0,1,，r2,其中r1=(N1-1)/2,r2=(N2-1)/2；頻率響應(yīng)變?yōu)?(3-14)其中,相頻響應(yīng)為-(r11+r22),是頻率1,2的線性函數(shù)(線性相位) ,幅頻響應(yīng)P(1,2)等于 (3-15)是獨(dú)立脈沖響應(yīng)的線性函數(shù)。用i|i=1,2,，n來表示獨(dú)立的脈沖響應(yīng),則幅頻響應(yīng)可表示為： (3-16)其中是w1,w2的某個(gè)函數(shù)。當(dāng)i=h(n1,n2),0n1<r1,0n2&l

21、t;r2時(shí) (3-17)當(dāng)i=h(n1,r2),0n1<r1時(shí)(3-18)當(dāng)i=h(r1,n2),0n2<r2時(shí)當(dāng)i=h(r1,r2)時(shí), =1。令=，T，=1,2,nT,則幅頻響應(yīng)可進(jìn)一步化為 (3-19)當(dāng)所設(shè)計(jì)的濾波器就是一個(gè)線性相位濾波器時(shí),只需考慮幅頻響應(yīng)的逼近問題,問題簡(jiǎn)化為:使P(1,2)按如下約束最小二乘準(zhǔn)則逼近期望幅頻響應(yīng)D(1,2)： (3-20) (3-21)其中1與2是頻率區(qū)域(1,2)-p1p,-p2p的子集,(1,2)為預(yù)定的逼近誤差上限函數(shù)。容易看到,P(1,2)關(guān)于1,2矩形對(duì)稱, 因此要求D(1,2)也具有矩形對(duì)稱性,并且只需考慮頻率子區(qū)域(1,2

22、)01p,02p上的逼近問題。此時(shí),逼近問題仍然表示為(1)、(2) ,但其中的1和2是頻率子區(qū)域的子集。一般情況下,式(4-16)的積分需采用數(shù)值方法計(jì)算,式(4-17)的約束用稠密頻率點(diǎn)集上的約束集代替。為此,我們將頻率子區(qū)域離散化為稠密的頻率點(diǎn)集=(1i,2j) =(ip/L,jp/L)|i,j=0,1,LL是一個(gè)充分大的正整數(shù)。將式(4-16)積分替換為1上的求和,將式(4-17)約束替換為2上的約束集。二維線性相位FIR濾波器的約束最小二乘設(shè)計(jì)問題可化為:6 (3-22)(3-23)3.3基于頻率抽樣法二維FIR濾波器的程序?qū)崿F(xiàn)Matlab信號(hào)處理工具箱中的firls函數(shù)用來設(shè)計(jì)具有

23、最小二乘線性相位FIR濾波器。該濾波器可使指定頻段的理想分段線性函數(shù)與濾波器頻響應(yīng)之間的誤差平方和最小7。b=firls(n,f,m)可得到n階FIR濾波器，其幅頻特性與由f和m指定的性能匹配。濾波器系數(shù)矢量b共包含有n+1個(gè)系數(shù)，且滿足對(duì)稱關(guān)系b(k)=b(n+2-k) k=1，n+1 對(duì)應(yīng)濾波器即為型和型線性相位濾波器。b=firls（n,f,m,w）可利用權(quán)值矢量w對(duì)各頻率段進(jìn)行加權(quán)擬合，w的長(zhǎng)度為f和m長(zhǎng)度的一半，用以指定各頻率段的權(quán)值。b=firls(n,f,m,ftype)和b=firls（f,m,w,ftype）可指定濾波器的類型：當(dāng)ftype=hilbert時(shí)，設(shè)計(jì)的濾波器為奇

24、對(duì)稱的線性相位濾波器（型和型），b的系數(shù)滿足以b(k)=-b(n+2-k),k=1,n+1。這類濾波器包括Hilbert變換器。當(dāng)ftype=differentiator時(shí)，也可設(shè)計(jì)出型和型濾波器，但采用了特殊的加權(quán)技術(shù)，對(duì)非零幅值的頻段，將誤差乘以(1+f)2，這樣可使低頻段誤差大大小于高頻段誤差。8fsamp2函數(shù)是使用頻率抽樣法設(shè)計(jì)二維FIR濾波器。h= fsamp2(hd),設(shè)計(jì)一個(gè)頻率響應(yīng)為hd的二維FIR濾波器，返回濾波器系數(shù)到矩陣h,濾波器h有一個(gè)頻率響應(yīng)經(jīng)過點(diǎn)hd。如果hd是在M、N之間，那么H 也在M、N之間。hd是一個(gè)包含理想頻率響應(yīng)采樣點(diǎn)的矩陣，采樣點(diǎn)在x和y頻率軸-1.

25、0 和1.0 之間的等間隔點(diǎn)上，1.0對(duì)應(yīng)采樣頻率的1/2, 或pi弧度。為了結(jié)果的準(zhǔn)確, 使用 FREQSPACE返回頻率點(diǎn)得到hd。h = fsamp2(f1,f2,hd,M N)產(chǎn)生一個(gè)從M到N 的FIR濾波器，使濾波器響應(yīng)hd在f1和f2之間。得到的濾波器盡可能接近理想的頻率響應(yīng)。最好的結(jié)果是必須有至少M(fèi)*N個(gè)理想的頻率點(diǎn)。如果您指定的點(diǎn)少于M*N點(diǎn),fsamp2將出現(xiàn)警告。3下面是利用該函數(shù)和二維FIR濾濾器的頻率采樣設(shè)計(jì)函數(shù)fsamp2進(jìn)行設(shè)計(jì)：設(shè)計(jì)要求：利用頻率抽樣法設(shè)計(jì)一20階的具有最小二乘線性相位的FIR帶阻濾波器。首先必須給出一維的相應(yīng)FIR濾波器。f=0 0.1 0.2

26、0.3 0.4 0.6 0.7 0.8 0.9 1;m=1 1 1 0 0 0 0 1 1 1；產(chǎn)生一個(gè)20階的具有最小二乘線性相位的一維FIR帶阻濾波器。B=firls(20,f,m,Hilbert);h0,w=freqz(b,1,512,2);Plot(f,m,-,w,abs(h0);在一維濾波器函數(shù)取理想帶阻濾波器時(shí)，采用頻率抽樣法產(chǎn)生一個(gè)20階的具有最小二乘線性相位的二維FIR帶阻濾波器h1。h1=fsamp2(f,f,m,5 5);figure(2)freqz2(h1,32 32);圖3.1 一維頻率FIR響應(yīng)圖3.1中虛線代表理想FIR頻率響應(yīng)曲線，實(shí)線代表實(shí)際的FIR頻率響應(yīng)曲線

27、。在一維波器函數(shù)取實(shí)際帶阻濾波器時(shí)，采用頻率抽樣法產(chǎn)生一個(gè)20階的具有最小二乘線性相位的二維FIR帶阻濾波器h2。ff=1/512:1/512:1;h2=fsamp2(ff,ff,h0,5 5);figure(3)freqz2(h2,32 32);圖3.2 理想二維FIR頻率響應(yīng)圖3.3 實(shí)際二維FIR頻率響應(yīng)3.4基于頻率抽樣法圖像濾波的實(shí)現(xiàn)數(shù)字濾波器被廣泛應(yīng)用于語音、圖像、無線電等領(lǐng)域，具有廣闊的發(fā)展空間，特別是在圖像處理中的應(yīng)用非常廣泛，因?yàn)楸救怂鶎W(xué)知識(shí)有限，只能利用濾波器對(duì)圖象進(jìn)行一些簡(jiǎn)單的處理，但我堅(jiān)信應(yīng)用濾波器處理圖象將是未來圖象處理的一大趨勢(shì)并占主導(dǎo)地位。在本設(shè)計(jì)中,將對(duì)自己設(shè)計(jì)

28、的濾波器進(jìn)行一些簡(jiǎn)單的應(yīng)用。在現(xiàn)實(shí)生活當(dāng)中濾波器越來越廣泛地被應(yīng)用于圖象處理領(lǐng)域，因此，在本次設(shè)計(jì)中將把濾波思想應(yīng)用于X線片圖像處理，旨在通過FIR濾波器快速高效的對(duì)X線片圖像進(jìn)行濾波處理，使處理后的圖片更具立體感。圖3.5 圖像濾波流程圖利用上面設(shè)計(jì)的濾波器進(jìn)行圖象分析，導(dǎo)入原始圖像，對(duì)于彩色圖片先要進(jìn)行灰度化10,圖像濾波的流圖如圖3.5所示。圖3.4 原始圖像f=imread(lena256.bmp);利用上面設(shè)計(jì)的FIR濾波器h1進(jìn)行濾波fh1=filter2(h1,f);figure(4)imshow(fh1);利用上面設(shè)計(jì)的FIR濾波器h2進(jìn)行濾波fh2=filter2(h2,f)

29、;figure(5)imshow(fh2);3.5 結(jié)果與展望3.5.1 設(shè)計(jì)結(jié)果與所取得成績(jī)圖3.6 理想二維FIR濾波結(jié)果 圖3.7 實(shí)際二維FIR濾波結(jié)果圖像濾波結(jié)果如圖3.6，圖3.7所示。所完成的工作：(1) 研究了數(shù)字濾波的理論知識(shí)，為系統(tǒng)整體奠定了理論基礎(chǔ)；(2)了解FIR濾波器的一些設(shè)計(jì)方法，利用頻率抽樣法設(shè)計(jì)具有最小二乘線性相位的FIR帶阻濾波器；(3)利用設(shè)計(jì)出的濾波器實(shí)現(xiàn)圖像濾波；(4)研究了MATLAB軟件在數(shù)字信號(hào)處理，尤其是數(shù)字濾波器處理中的應(yīng)用。3.5.2 研究展望在對(duì)醫(yī)學(xué)知識(shí)的深入了解下,上述濾波器的快速處理X線片圖象的特點(diǎn)可用于快速高效的去除人體手腕部X線片的

30、噪聲，從而方便醫(yī)學(xué)人員準(zhǔn)確評(píng)價(jià)骨齡。下面將簡(jiǎn)要介紹骨齡評(píng)價(jià)，以便更好的理解這一應(yīng)用背景。人的骨骼年齡簡(jiǎn)稱為骨齡，是一種與生活中人們常說的“時(shí)間年齡” 不一致的、能準(zhǔn)確反映人體生長(zhǎng)發(fā)育信息的生物學(xué)年齡。骨齡是出生后大多數(shù)正常兒童隨年齡增長(zhǎng)而出現(xiàn)的有規(guī)律的骨骺線解剖變化標(biāo)志，是由于個(gè)體間營(yíng)養(yǎng)發(fā)育的不同以與地區(qū)、種族、性別的不同而表現(xiàn)出的與時(shí)間年齡的不一致。另外，通過研究大量X線片樣本，比較各地區(qū)兒童或成年人骨骼年齡，找出其成長(zhǎng)發(fā)育和身體狀況規(guī)律，對(duì)于醫(yī)學(xué)研究意義更大。通過本次設(shè)計(jì)，自己從中取得了一些成績(jī)，理論水平也得到了一定的提高，同時(shí)也暴露了一些問題。首先，對(duì)一個(gè)課題必須要閱讀大量的文獻(xiàn)和書籍來

31、獲得一定的感性認(rèn)識(shí)，然后才能有自己的想法，這是一條必經(jīng)之路。其次，理論基礎(chǔ)知識(shí)很重要，論文設(shè)計(jì)了很多的算法，會(huì)用到很多基礎(chǔ)知識(shí)，如果用的時(shí)候再去學(xué)會(huì)浪費(fèi)時(shí)間。最后，要有信心，遇到困難要向別人請(qǐng)教，這樣可以大大加快研究進(jìn)程。以上是我做論文的一些心得體會(huì)，這些對(duì)我以后的工作會(huì)有很大的幫助。附錄A 程序代碼%直接法設(shè)計(jì)低通濾波器M=20;alpha=(M-1)/2; l=0:M-1;wl=(2*pi/M)*l;Hrs=1,1,1,zeros(1,15),1,1;Hdr=1,1,0,0;wdl=0,0.25,0.25,1;k1=0:floor(M-1)/2);k2=floor(M-1)/2)+1:M-

32、1;angH=-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2);H=Hrs.*exp(j*angH);h=real(ifft(H,M);db,mag,pha,grd,w=freqz_m(h,1);Hr,ww,a,L=Hr_Type2(h);subplot(2,2,1);plot(wl(1:11)/pi,Hrs(1:11),'o',wdl,Hdr);axis(0,1,-0.1,1.1);title('Frequency Samples:M=20')xlabel('frequency in pi units');yla

33、bel('Hr(k)')subplot(2,2,2);stem(l,h);axis(-1,M,-0.1,0.3)title('Impulse Response');xlabel('n');ylabel('h(n)');subplot(2,2,3);plot(ww/pi,Hr,wl(1:11)/pi,Hrs(1:11),'o');axis(0,1,-0.2,1.2);title('Amplitude Response')xlabel('frequency in pi units');y

34、label('Hr(w)')subplot(2,2,4);plot(w/pi,db);axis(0,1,-60,10);gridtitle('MagnitudeResponse');xlabel('frequencyinpiunits');ylabel('Decibels');最優(yōu)設(shè)計(jì)法T1=0.5clear all;M=40;alpha=(M-1)/2; l=0:M-1;wl=(2*pi/M)*l;Hrs=ones(1,5),0.5,zeros(1,29),0.5,ones(1,4);Hdr=1,1,0,0;wdl=0,0.2,

35、0.3,1;k1=0:floor(M-1)/2);k2=floor(M-1)/2)+1:M-1;angH=-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2);H=Hrs.*exp(j*angH);h=real(ifft(H,M);db,mag,pha,grd,w=freqz_m(h,1);Hr,ww,a,L=Hr_Type2(h);subplot(2,2,1);plot(wl(1:21)/pi,Hrs(1:21),'o',wdl,Hdr);axis(0,1,-0.1,1.1);title('Frequency Samples:M=40

36、9;)xlabel('frequency in pi units');ylabel('Hr(k)')subplot(2,2,2);stem(l,h);axis(-1,M,-0.1,0.3)title('Impulse Response');xlabel('n');ylabel('h(n)');subplot(2,2,3);plot(ww/pi,Hr,wl(1:21)/pi,Hrs(1:21),'o');axis(0,1,-0.2,1.2);title('Amplitude Response&

37、#39;)xlabel('frequency in pi units');ylabel('Hr(w)')subplot(2,2,4);plot(w/pi,db);axis(0,1,-120,10);gridtitle('MagnitudeResponse');xlabel('frequencyinpiunits');ylabel('Decibels');%最優(yōu)設(shè)計(jì)法T1=0.39clear all;M=40;alpha=(M-1)/2; l=0:M-1;wl=(2*pi/M)*l;Hrs=ones(1,5),0.3

38、9,zeros(1,29),0.39,ones(1,4);Hdr=1,1,0,0;wdl=0,0.2,0.3,1;k1=0:floor(M-1)/2);k2=floor(M-1)/2)+1:M-1;angH=-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2);H=Hrs.*exp(j*angH);h=real(ifft(H,M);db,mag,pha,grd,w=freqz_m(h,1);Hr,ww,a,L=Hr_Type2(h);subplot(2,2,1);plot(wl(1:21)/pi,Hrs(1:21),'o',wdl,Hdr);axis

39、(0,1,-0.1,1.1);title('Frequency Samples:M=40')xlabel('frequency in pi units');ylabel('Hr(k)')subplot(2,2,2);stem(l,h);axis(-1,M,-0.1,0.3)title('Impulse Response');xlabel('n');ylabel('h(n)');subplot(2,2,3);plot(ww/pi,Hr,wl(1:21)/pi,Hrs(1:21),'o'

40、);axis(0,1,-0.2,1.2);title('Amplitude Response')xlabel('frequency in pi units');ylabel('Hr(w)')subplot(2,2,4);plot(w/pi,db);axis(0,1,-120,10);gridtitle('MagnitudeResponse');xlabel('frequency in pi units');ylabel('Decibels');增加兩個(gè)過渡帶clear all;M=60;alpha=(M-1)/2; l=0:M-1;wl=(2*pi/M)*l;