基于MATLAB的語音信號濾波處理畢業(yè)設計(論文).doc

學畢業(yè)設計(論文)題目: 基于matlab的語音信號濾波處理作 者 夏杰 屆 別 2015屆 院 別 信息與通信工程學院 專 業(yè) 電子信息工程 指導教師 羅朝明 職 稱 副教授 完成時間 2015年5月17日 摘要隨著計算機和信息科學的飛速發(fā)展，信號處理逐漸發(fā)展成為一門獨立的學科，成為信息科學的重要組成部分，在語音處理、雷達、圖像處理、通信、生物醫(yī)學工程等眾多領域中得到廣泛應應用。本論文研究了濾波器的設計，并通過設計好的濾波器語來分析語音信號。在計算機中錄入一段語音信號，通過matlab對語音信號進行時頻變換和分析；再通過matlab構造一定頻率的噪音信號，再將兩個信號進行合成；然后設計不同類型的數(shù)字濾波器對其噪音進行濾波處理。關鍵詞：matlab；語音信號；濾波器；仿真；abstractwith the rapid development of computers and information science, signal processing is becoming an independent discipline, and playing an important part in information science. it has been widely applied in speech processing, radar, image processing, communication, bio-medical engineering and many other fields. this paper has studied the filter design, and it analyses the voice signal through the designed filter. after recording a voice signal in the computer, then doing the time-frequency transformation and analyses for the voice signal through matlab. we construct a fixed-frequency signal, and synthesis the two signals, and it deals with filter processing on the speech signal with noise, and analyses the signals in the time domain and frequency domain.keyword: matlab；voice signal；filter；simulation；目 錄摘要iabstractii第一章 緒論11.1 語音信號濾波處理的簡介11.2 數(shù)字濾波器的概念11.3 matlab軟件的介紹31.4 論文各章節(jié)的安排4第二章 語音信號的頻譜分析52.1 快速傅里葉算法的原理52.2語音信號的采集及頻譜分析102.3 語音信號的加噪處理及頻譜分析12第三章 數(shù)字濾波器的設計153.1數(shù)字濾波器的設計原理153.2 用窗函數(shù)法設計fir數(shù)字濾波器183.3 用雙線性法設計iir數(shù)字濾波器21第四章 語音信號的濾波244.1 語音信號的濾波測試244.2 加噪語音信號的濾波處理27第五章 總結和展望295.1 課題總結295.2 課題展望29參考文獻30致謝31ii第1章 緒論1.1 語音信號濾波處理的簡介語音信號的基本組成單位是“音素”。音素主要由“濁音”和“清音”這兩大類組成。如果將語音只有背景噪聲稱為“無聲”，那么可以將音素分成“清音”、“濁音”、“無聲”這三類。一個音節(jié)由輔音和元音兩部分組成。因為元音都是濁音，所以元音在音節(jié)中占據(jù)著主要的部分。在普通話中，“輔音元音”構成了每個音節(jié)。語音信號的濾波處理是數(shù)字信號處理領域目前發(fā)展最為迅速的信息科學研究領域的核心技術之一，通過語音傳遞信息是人類最重要、最有效、最常用和最方便的交換信息形式。語音信號處理是現(xiàn)代信息處理的基本內(nèi)容，數(shù)字信號的處理更是重中之重。數(shù)字信號處理的研究內(nèi)容主要是語音信號和圖像信號，而研究語音信號對于現(xiàn)代語音通信有著積極的意義。研究語音信號又分為時域和頻域兩個方面。（1） 語音信號的時域分析處理：對語音信號進行時域分析的本質(zhì)就是分析和提取語音信號的時域參數(shù)。當進行語音分析時，它的時域波形圖是很直觀的。其實語音信號的實質(zhì)就是時域信號，所以使用時域分析是最早也是最廣泛的一種分析方法。時域分析一般用于參數(shù)的分析和應用，如語音的大分類、分割、預處理等。這種分析方法有以下特點：第一，比較直觀的語音信號和明確其物理意義。第二，可以簡單的實現(xiàn)并且運算少。第三，獲得語音信號的重要參數(shù)。最后，它的實現(xiàn)用到了示波器等常規(guī)設備，所以應用起來較簡單。（2）語音信號的頻譜分析處理：語音信號的頻域分析就是對語音信號的頻域持征進行分析。語音信號的頻域分析包括語音信號的功率譜、倒頻譜、頻譜、頻譜包絡分析等，而一般常用的頻域分析方法包括帶傅里葉變換法、通濾波器組法等幾種。本文介紹的是語音信號的傅里葉分析法。1.2 數(shù)字濾波器的概念數(shù)字信號處理是一門發(fā)展迅速、應用廣泛的前沿性學科，其理論和實踐性都很強。為了系統(tǒng)的掌握和理解信號的收集、處理、傳輸、顯示和儲存，本課題介紹了用matlab來處理語音信號。在21世紀這個信息技術高速發(fā)展的時代，信號處理逐漸趨向于數(shù)字化、軟件化這方面發(fā)展。數(shù)字信號處理主要是對語音信號和圖像信號的研究，其中，語音信號的研究可以分為從時域和頻域這兩個方面來進行。時域的分析處理有以下兩種方法：一種是對語音信號進行分析，這一分析方法屬于線性處理的范疇，其原理主要是通過信號的卷積、倍乘、時移、加減、以及求相關函數(shù)等來實現(xiàn)；另外一種是通過處理生成和變換成各種調(diào)制信號，這一分析方法屬于非線性的范疇，其原理是壓縮和擴張信號平均累加器的動態(tài)范圍，用門限法對噪音進行抑制。通過對信號的頻域分析處理，可以對信號的頻率特性在頻譜中加以分析研究，這是對不確定信號分析的主要方法，主要是擴展了信號分析的范圍。通常在世紀應用中，信號的時、頻域分析經(jīng)常同時進行。數(shù)字濾波是數(shù)字信號處理理論的一部分。具體來說，凡是用數(shù)字方式對信號進行濾波、變換、調(diào)制、壓縮、增強、均衡、解調(diào)、估值、識別、產(chǎn)生等加工處理，都可以納入數(shù)字信號處理的領域。 數(shù)字信號處理這門學科主要是對數(shù)字濾波器設計方法進行研究。對于數(shù)字濾波器，早在49年代末期，就有人對其可能性進行過研究討論，在50年代的時候，也有人研究討論過數(shù)字濾波這個具有跨時代意義的問題。但可惜的是直到60年代中期，人們才是開始慢慢形成關于數(shù)字濾波器的一套完整的理論，在這一時期，通過不同形式的數(shù)字濾波結構（有的以最小運算誤差為特點，有的以最快運算速度為特點，而有的則囊括以上兩者）。在數(shù)字濾波器領域中，有限沖激響應（fir）和無限沖激響應（iir）的關系認識轉(zhuǎn)化是一個重大的發(fā)現(xiàn)。在初期，人們一般認為iir濾波器的運算效率比fir濾波器更高，所以明顯青睞于前者，但當有人提出可以用快速傅里葉變換（fft）來實現(xiàn)其卷積運算的理念后，發(fā)現(xiàn)高階fir濾波器統(tǒng)一可以用比較高的運算效率來實現(xiàn)，這就激勵人們大量的研究和探索數(shù)字濾波器的fir高性能濾波器的設計方法和頻域設計方法，從而有了數(shù)字濾波器設計中頻域方法和時域方法的出世。然而，這些數(shù)字濾波器都是一些人們的早期研究。早期的數(shù)字濾波器盡管在聲納、語音、地震和醫(yī)學和的信號處理中曾經(jīng)起過非常大的作用，但是出于當時計算機主機的價格比較昂貴，在一定程度上阻礙了個別專門用途數(shù)字濾波器的發(fā)展。在科學技術蓬勃發(fā)展的70年代，引進了計算機輔助的設計方法，這使數(shù)字濾波器的設計以及實現(xiàn)僅僅是逼近響應模擬濾波器。通常來說，變換模擬濾波器函數(shù)來設計數(shù)字濾波器，幾乎很難逼近任意沖激響應或頻率響應，然而采用計算機輔助設計這一方法則可能實現(xiàn)時域或頻域的最佳逼近，或時域頻域的聯(lián)合最佳逼近，這樣，不但可以使數(shù)字濾波器的分析和設計內(nèi)容變的豐富起來，而且各種新的數(shù)字信號處理系統(tǒng)也都可以用專用數(shù)字硬件加以實現(xiàn)。1.3 matlab軟件的介紹matlab是由美國mathworks公司推出的用于圖形處理和數(shù)值計算的數(shù)學計算環(huán)境。在matlab環(huán)境下，用戶可以集成地進行程序設計、數(shù)值計算、輸入輸出、圖形繪制、文件管理等各項操作。所以相比其他同類軟件，它卓越的數(shù)據(jù)可視化和優(yōu)秀的數(shù)字計算能力使其脫穎而出。最初的matlab系統(tǒng)是由cleve moler公司用fortran語言設計的，現(xiàn)在matlab程序是基于mathworks公司在c語言環(huán)境下開發(fā)的。1984年開始發(fā)行它的第一版（dos版本1.0），而后經(jīng)過10多年的不斷完善和改進，版本不斷升級 ，其功能越來越豐富，工具也越來越豐富，因此它應用的范圍也越來越廣發(fā)。matlab語言是一種對數(shù)據(jù)分析和處理功能十分強大的據(jù)算計應用軟件，它可以將聲音文件變換為離散的數(shù)據(jù)文件，然后利用其強大的矩陣運算能力處理數(shù)據(jù)，如數(shù)字濾波、時域和頻域分析、傅立葉變換、聲音回放以及各種圖的呈現(xiàn)等，信號處理是matlab的最重要的應用領域之一。matlab語言之所以能如此迅速地普及是由于它有著不同于其他語言的特點：1.簡潔緊湊的語言使用非常方便，并且?guī)旌瘮?shù)非常豐富。matlab程序書寫形式自由，利用其豐富的庫函數(shù)巧妙的避開了子程序編程任務的復雜，壓縮了一切不必要的編程工作。2.運算符豐富。由于matlab是用c語言編寫的，所以matlab提供了和c語言幾乎一樣多的運算符，靈活使用matlab的運算符將使程序變得極為簡短。3.matlab既具有結構化的控制語句，又有面向?qū)ο缶幊痰奶匦浴?.語法限制不嚴格，程序設計自由度大。5.程序的可移植性很好，基本上不做什么修改就可以運行在各種型號的計算機和操作系統(tǒng)上。6.matlab的圖形功能強大。在c和fortran語言里，繪圖都不是很容易，但是在matlab里通過數(shù)據(jù)的可視化則顯得非常簡單。此外，matlab編輯圖形界面的能力非常強。7.matlab的一項特色是擁有功能強大的工具箱。matlab包含兩個部分：核心部分和各種可選的工具箱。核心部分有數(shù)百個核心內(nèi)部函數(shù)。本設計只使用了上面的部分功能，就是用matlab對含噪的語音信號同時在時域和頻域進行濾波處理和分析。1.4 論文各章節(jié)的安排第一章 介紹了數(shù)字信號處理的發(fā)展趨勢和所用的的開發(fā)軟件matlb。第二章 介紹了語音信號的采集和頻譜分析的方法和原理。第三章 介紹了數(shù)字濾波器的設計方法。第四章 介紹了用不同的數(shù)字濾波器對語音信號進行濾波處理。第五章 介紹了本人對這次畢業(yè)設計的總結和展望未來數(shù)字濾波器的發(fā)展。第2章 語音信號的頻譜分析2.1 快速傅里葉算法的原理dft是信號分析與處理中的一種重要變換。因為dft直接計算的計算量與變換區(qū)間長度（n的平方）成正比，當n較大時，造成計算量太大，因而在快速傅里葉變換（簡稱fft）出現(xiàn)以前，直接用dft算法對信號進行頻譜分析和它的實時處理是不切實際的。直到1965年發(fā)現(xiàn)了dft這種快速算法，這種情況才發(fā)生了逆轉(zhuǎn)。自1965年圖基（j.w.tuky）和庫力(t.w.coody)在書名為計算機數(shù)學（math.computation,vol,19.1965)雜志上發(fā)表了著名的一篇機器計算傅里葉級數(shù)的一種算法論文后，桑德（g.sand）-圖基等快速傅里葉算法的相繼出現(xiàn)，在人們的不斷進行改進和不斷完善，便很快形成了一套高效運算方法，從而才有了我們現(xiàn)在的快速傅里葉變換，簡稱為fft（fast fourier transform）。dft的運算效率在這種算法的基礎上可以提高12個數(shù)量級，這為數(shù)字信號處理技術應用于各種信號的實時處理創(chuàng)造了良好的條件，極大的推動了數(shù)字信號處理技術在當今社會的發(fā)展。fft算法的概念：(1)“級”概念：將n 點dft先分成兩個n/2點dft，再是四個n/4點dft直至n/2個兩點dft每分一次稱為“一”級運算。因為n=2m所以n點dft可分成m級。(2)“組”概念：每一級都有n/2個蝶形單元，例如：n=16,則每級都有8個蝶形單元。每一級的n/2個蝶形單元可以分成若干組，每一組具有相同的結構，相同的 因子分布，第m級的組數(shù)為： n/2m+1 n=16=24，分4級。(3)因子的分布: 每由后向前（m由m-1-0級）推進一級，則此系數(shù)為后級系數(shù)中偶數(shù)序號的那一半。(4)按時間抽取法：由于每一步分解都是基于在每級按輸入時間序列的次序是屬于偶數(shù)還是奇數(shù)來分解為兩個更短的序列，所以稱為“按時間抽取法”。如圖2.1所示fft算法基本上分為兩大類：時域抽取法fft（decimation-in-time fft,簡稱dit-fft）和頻域抽取法fft（decimation-in-frequency fft,簡稱dif-fft）。下面就介紹dif-fft算法。設序列長度為,且滿足 (為正整數(shù)）按的奇偶性把分解成為兩個點的子序列：2點2點2點2點2點2點2點2點兩個2點dft兩個2點dft兩個2點dft兩個2點dft兩個4點dtf兩個4點dft兩個n/2點的dft2.1按時間抽取法 （2.1）則可將dft化為 （2.2）由于故式2.2可表示成 （2.3）與分別是及的點dft： （2.4） （2.5）由此，我們可以看到，一個點dft已分解成兩個的dft。這兩個點的dft再按照2.4和2.5式組合成一個點dft。 （2.6）這里應該看到,只有個點，即。而x(k)卻有個點，即,故用上式計算得到的只是的前一半結果，要用,來表達全部的值，還必須應用系數(shù)的周期性，即這樣可得到 （2.7）同理可得 （2.8）前兩式說明了后半部分值（）所對應的，分別等于前半部分值（）所對應的，。再考慮的以下性質(zhì)： 這樣，把上述各式帶入，就可以將表達為前后兩部分： k=0,1, （2.9）因此，只要求出0到區(qū)間的所有，值，即可求出0到區(qū)間內(nèi)的所有值，顯然節(jié)省了運算量?？梢杂玫涡盘柫鲌D表示。 一個點dft可以分解為兩個點dft，每一個點dft只需次復數(shù)乘法以及次復數(shù)加法。兩個點dft共需次復數(shù)乘法和次復數(shù)加法。把兩個點dft和成為點dft時，其中有個蝶形運算還需要次復數(shù)乘法及次復數(shù)加法。因而通過第一次分解后，總共還需要次復數(shù)乘法以及次復數(shù)加法。由因此，通過這樣分解后，它的運算工作量至少節(jié)省了一半。既然這樣分解是有效的，由于，因而仍是偶數(shù)，所以可以再一次把每個點子序列再其奇偶部分分解為兩個點的子序列。 （2.10） （2.11）且： 式中： （2.12） （2.13）根據(jù)上面同樣的分析可以得出：利用四個點的dft及兩級蝶形組合運算來計算點dft，相比只用一次分解蝶形組合方式的計算量又減少了一大半。它的算法步驟如下：（1）分組，變量置換：dft變換： 先將按的奇偶分為兩組，作變量置換:n=偶數(shù)時，令n=2r;得到：則其dft可化為兩部分：前半部分： 后半部分：代入dft中 （2）代入dft變換式：（3）求出子序列的dft： （2.14）（4）由上面的公式我們可以知道：一個點的dft被分解為兩個點dft。，這兩個點的dft按照式2.15又合成點dft中的前半部分再應用系數(shù)的周期性，求出用，表達的后半部的值。 （2.15）（5）求出后半部的表示式： （2.16）看出：后半部的值所對應的，則完全重復了前半部分的值所對應的，的值。又（6）最后得出：頻域中的個點頻率成分為：（2.17）結論：只要求出區(qū)間內(nèi)的各個整數(shù)值所對結應的，值，即可以求出在整個區(qū)間內(nèi)全部值，這就是為什么fft能大量節(jié)省計算的關鍵所在。由于，因此仍為偶數(shù)，依照上面的方法可以進一步把每個點子序列，再按輸入的奇偶分解為兩個點的子序列，再按這種方法不斷劃分下去，直到最后剩下的是2點dft，兩點dft實際上只是加減運算。上面頻域中前/后半部分表示式可以用蝶形信號流圖表示。如圖2.2所示 2.2 蝶形信號流圖2.2語音信號的采集及頻譜分析采樣頻率在進行模擬數(shù)字信號的轉(zhuǎn)換過程中，當采樣頻率fs.max大于信號中最好頻率fmax的2倍時，即：fs.max=2fmax，則采樣之后的數(shù)字信號完整的保留了原始信號中的信息，一般實際應用中保證采樣頻率為信號最高頻率的510倍，采樣定理又稱為耐奎斯特定理。1924年耐奎斯特（nyquist）就推導出在理想低通信道的最高大碼元傳輸速率的公式：理想低通信道的最高大碼元傳輸速率=2w*log2n（其中w是理想低通信道的帶寬，n是電平強度） 采樣位數(shù)即采樣值或取樣值，用來衡量聲音波動變化的參數(shù)，是指聲卡在采集和播放聲音文件時所使用數(shù)字聲音信號的二進制位數(shù)。采樣頻率是指錄音設備在一秒鐘內(nèi)對聲音信號的采集次數(shù)，采樣頻率越高，聲音的還原就越真實越自然。采樣位數(shù)和采樣頻率對于音頻接口來說事最為重要的兩個指標，也是選擇音頻接口的兩個重要標準。無論采樣頻率如何，理論上來說采樣的位數(shù)決定了音頻數(shù)據(jù)的最大力度范圍。每增加一個采樣位數(shù)，相當于力度范圍增加了6db。采樣位數(shù)越多則捕捉到的信號越精確。對于采樣率來說你可以想象它類似于一個照相機，44.1khz意味著音頻流進入計算機時計算機每秒會對其拍照達到441000次。顯然采樣頻率越高，計算機攝取的圖片越多，對于原始音頻的還原也越加精確。 該設計采用的是優(yōu)美的鈴聲為分析樣本，利用windows下的錄音機錄制一小段語音，然后在matlab軟件平臺下，利用wavread函數(shù)對語音信號進行采樣，采樣頻率為fs=44100hz，得到的采樣頻率和采樣點數(shù)為：y,fs,bit=wavread(c:代碼wy.wav); 采樣頻率：fs=44100hz bits=16 采樣點數(shù)：n=1329408 對原始語音信號進行采集并播放： x1,fs,bits=wavread(f:yuyinpsb.wav); sound(x1,fs,bits); 信號頻譜分析是將信號源發(fā)出的信號強度按頻率順序展開，使其成為頻率的函數(shù)，并考察變化規(guī)律，稱為頻譜分析。頻譜分析主要分析信號是由哪些頻率的正弦信號疊加得到的，以及這些正弦信號的振幅。信號頻譜分析是將信號源發(fā)出的信號強度按頻率順序展開，使其成為頻率的函數(shù)，并考察變化規(guī)律，稱為頻譜分析。頻譜分析主要分析信號是由哪些頻率的正弦信號疊加得到的，以及這些正弦信號的振幅。首先畫出語音信號的時域波形，然后用函數(shù)fft對語音信號進行傅立葉變換，得到信號的頻譜特性。然后通過sound函數(shù)回放原語音信號，我們聽到一段優(yōu)美的英文歌。其程序如下clearclcx1,fs,bits=wavread(c:代碼wy.wav);size(x1)%求y的行數(shù)和列數(shù)y1=x1(:,1);%對信號進行分列處理(單聲道處理)，前面表示所有行,后面表示取單列n1=length(y1);%取y的長度，也就是總采樣點數(shù)t=1/fs;%采樣間隔ts=t*1000%化成mst1=(0:n1-1)*ts;%設置波形圖橫坐標為時間figure(1);%創(chuàng)建繪圖窗口plot(t1,y1);%畫出音樂信號時域波形圖xlabel(時間（ms）);%橫坐標命名ylabel(化到-1到1之間的幅值（v）);%縱坐標命名 title(原始音樂信號時域圖)%標題n=length(x1) %計算音頻信號的長度figure(2);ftemp=-n/2:n/2-1.*fs/n;y1=fft(x1); %做length(x1)點的fftplot(ftemp,fftshift(abs(y1); %畫出原始語音信號的頻譜圖title(原始語音信號的頻譜);grid on;sound(x1,fs,bits);如圖3.1和3.2所示 圖3.1 原信號時域圖 圖3.2 原信號頻譜圖2.3 語音信號的加噪處理及頻譜分析在matlab中，運用randn函數(shù)產(chǎn)生與原信號等長度的隨機噪音信號。在通過高通濾波器濾除4000hz以前的噪音信號，在將其頻譜與原信號頻譜疊加得到最終的噪音信號。其程序如下：n=length(x1) %計算音頻信號的長度noise=0.1*randn(n,2); %產(chǎn)生等長度的隨機噪聲信號(這里的噪聲的大小取決于隨機函數(shù)的幅度倍數(shù)）x1_z=noise; %將兩個信號疊加成一個新的信號加噪聲處理 n=length(x1); %選取變換的點數(shù) x1_zp=fft(x1_z,n); %對n點進行傅里葉變換到頻域f=fs*(0:n/2-1)/n; % 對應點的頻率figure(2)subplot(2,1,1);plot(x1_z); %加噪語音信號的時域波形圖title(加噪語音信號時域波形);xlabel(時間軸)ylabel(幅值 a)subplot(2,1,2);plot(f,abs(x1_zp(1:n/2); %加噪語音信號的頻譜圖title(加噪語音信號頻譜圖);xlabel(頻率hz);ylabel(頻率幅值);%提取4000hz以后的噪音信號%wp1=2*pi*3800/fs;ws1=2*pi*4000/fs;wc1=(wp1+ws1)/2/pi;bt=ws1-wp1;n1=ceil(1.8*pi/bt);n2=n1+mod(n1+1,2);window=boxcar(n2); %長度為n1的矩形窗windowb1=fir1(n2-1,wc1,high,window);x1_high = filter(b1,1, x1_z); %對信號進行低通濾波figure(5);ftemp=-n/2:n/2-1.*fs/n;plot(ftemp,fftshift(abs(fft(x1_high);title(信號經(jīng)過高通濾波器（頻域）);z1=y1+fft(x1_high);figure(6);plot(ftemp,fftshift(abs(z1);title(疊加信號的頻譜)grid on;如圖3.3和3.4所示 圖3.3 噪音信號的頻譜圖 圖3.4 疊加信號的頻譜圖第3章 數(shù)字濾波器的設計3.1數(shù)字濾波器的設計原理濾波器分為模擬濾波器和數(shù)字濾波器。模擬濾波器由電阻，電容，電感等原器件構成；而數(shù)字濾波器則可由數(shù)字電路實現(xiàn)或軟件實現(xiàn)。按脈沖響應來分類，數(shù)字濾波器分為iir和fir，即無限沖激響應濾波器和有限沖激響應濾波器；其中iir網(wǎng)絡中有反饋回路，fir網(wǎng)絡中沒有反饋回路。從小的方面分：a.按所處理的信號分為兩種：模擬濾波器和數(shù)字濾波器。b.按所通過信號的頻段分為低通、高通、帶通和帶阻濾波器四種。 數(shù)字濾波器按單位取樣響應時域特性可分為:iir(infinite impulse response) 濾波器(即無限長單位沖激響應濾波器)fir(finite impulse response)濾波器(即有限長單位沖激響應濾波器)。低通、高通、帶通和帶阻濾波器都稱作是經(jīng)典濾波器，它們每一種又有模擬濾波器和數(shù)字濾波器兩種形式。 數(shù)字濾波器是數(shù)字信號處理的基礎,可以對信號進行過濾、參數(shù)估計和檢測等處理，在圖像、通信、雷達、語音等許多領域都有著廣泛的應用。特別是在數(shù)據(jù)壓縮和圖像處理等方面取得了非常大的成就。所以數(shù)字濾波器的設計顯得非常重要。matlab的出現(xiàn)為數(shù)字濾波的研究和應用提供了一個很好的平臺。它在以矩陣運算的基礎上，把程序設計、計算和可視化融合到了一個交叉式的環(huán)境中。數(shù)字濾波器用來對信號進行過濾、檢測與參數(shù)估計等處理,在通信、圖像、語音、雷達等許多領域都有著十分廣泛的應用。尤其在圖像處理、數(shù)據(jù)壓縮等方面取得非常大的成就。這使各個領域的研究人員可以直觀方便地進行科學研究與工程應用。當然設計一數(shù)字濾波器首先要清楚其原理。fir線性相位濾波器的原理如下:設是長度為的線性相位濾波器的脈沖響應，它的系統(tǒng)函數(shù)為：它的頻率響應函數(shù)為：它具有一個線性相位約束：其中，為一個常數(shù)。則有以下面的特性：稱之為對稱脈沖響應；或者稱之為反對稱脈沖響應。 于是,根據(jù)值的不同以及m的奇偶性,就分別產(chǎn)生了四種類型的線性fir相位濾波器.它們是:1型線性相位fir濾波器:為奇數(shù), 以中心點對稱;2型線性相位fir濾波器:為偶數(shù), 以點與點之間的中心對稱;3型線性相位fir濾波器:為奇數(shù), 以點為中心反對稱;4型線性相位fir濾波器:為偶數(shù), 以點與點之間的中心反對稱;根據(jù)這四類fir濾波器,可得到其響應的頻率響應的特征:其中為振幅響應,與之區(qū)別的是它與幅值特性不同,前者的取值可正可負,而后者只能為非負值。iir數(shù)字濾波器具有以下特點：(1) iir數(shù)字濾波器的系統(tǒng)函數(shù)可以寫成封閉函數(shù)的形式。 (2) iir數(shù)字濾波器采用遞歸型結構結構上帶有反饋環(huán)路。iir濾波器運算結構一般通常由相加、延時和乘以系數(shù)等基本運算組成，即可以構成級聯(lián)型、直接型、并聯(lián)型、正準型四種結構形式，它們都是具有反饋回路這一結構。由于在運算中會出現(xiàn)舍入處理，會造成誤差不斷累積，偶爾會產(chǎn)生微弱的寄生振蕩。 (3) iir數(shù)字濾波器在設計上可以借鑒成熟的模擬濾波器成果，如橢圓濾波器、巴特沃斯和契比雪夫等，所以在設計模擬器的時候可以查看它們的設計數(shù)據(jù)以及圖表，這樣我們的設計工作量相對而言就少了許多，其次，對計算工具的要求也不高。通過根據(jù)指標寫出模擬濾波器的公式以及通過一定變換，可以將模擬濾波器轉(zhuǎn)換為數(shù)字濾波器的公式。 (4) iir數(shù)字濾波器是非線性的，所以其相位特性不好控制，當對相位要求較高時，則需要添加相位校準網(wǎng)絡。 iir數(shù)字濾波器可用一個n階差分方程表示： （2.2） 或用它的z域系統(tǒng)函數(shù): （2.3） 對照模擬濾波器的傳遞函數(shù): （2.4） 可以看出,設計數(shù)字濾波器的思路與模擬濾波器相仿,設計實質(zhì)都是尋找一組系數(shù)b,a,去逼近所求的頻率響應,以便在性能上滿足預定的技術要求；不同的是數(shù)字濾波器是在z平面尋找合適的h(z), 模擬濾波器則是在s平面用數(shù)學逼近法尋找近似的所需特性h(s)。由于模擬濾波器的設計在理論上已十分成熟,因此數(shù)字濾波器設計的關鍵是將h(s)變換至h(z),即利用復值映射將模擬濾波器離散化。3.2 用窗函數(shù)法設計fir數(shù)字濾波器 窗函數(shù)設計法是一種通過截短和計權的方法使無限長非因果序列成為有限長脈沖應響應序列的設計方法。一般在設計濾波器之前，濾波器的技術指標是要根據(jù)具體的工程應用確定的。通常在運用中，因為數(shù)字濾波器往往被用來實現(xiàn)選頻操作，所以指標形式通常為在頻域中以分貝值給出的相應響應和幅度響應。窗函數(shù)法設計fir濾波器的步驟如下：(1) 根據(jù)阻帶衰減及過度帶寬，選擇窗函數(shù)的類型并估計其窗口長度n（或階數(shù)m=n-1）。窗函數(shù)的類型可以根據(jù)最小阻帶衰減as獨立選擇，原因是因為窗口長度n的大小對最小阻帶衰減as幾乎沒有影響。一旦窗函數(shù)類型確定以后，可以根據(jù)過度帶寬小于給定指標確定所擬用的窗函數(shù)的窗口長度n，假設待設計的濾波器的過渡帶寬為，它與窗口長度n的大小近似成反比，所以只要窗函數(shù)類型確定后，那么它的計算公式也就確定了。但是這些公式都是近似的，得出的窗口長度n還需要在計算中不斷修正以達到最好的要求。原則是：在保證阻帶衰減滿足其規(guī)定要求的情況下，選擇比較小的n，在n和窗函數(shù)類型確定后，就可以調(diào)用matlab軟件中的窗函數(shù)求出窗函數(shù)w（n）。 (2) 根據(jù)待求濾波器的理想頻率響應求出理想單位脈沖響應hd（n），如果給出基于窗函數(shù)法的fir數(shù)字低通濾波器設計。在一般情況下，hd（n）是不能用封閉公式表示的，需要采用數(shù)值方法表示；從w=0到w=2采樣n點，采用離散傅里葉反變換（idft）即可求出。 (3) 計算濾波器的單位脈沖響應h（n），它是理想單位脈沖響應和窗函數(shù)的乘積。 (4) 計算技術指標是否滿足要求，為了計算數(shù)字濾波器在頻域中的特性，可調(diào)用matlab中的freqz函數(shù)，如果不滿足要求的情況下，可根據(jù)具體情況適當調(diào)整窗函數(shù)長度或類型，直到滿足要求為止。 通常上述的幾點難以同時滿足。實際中設計fir數(shù)字濾波器往往要求是線性相位的，因此要求)(nw滿足線性相位的條件，即要求w(n)滿足：w(n)=w(n-1-n) （2） 所以，窗函數(shù)不僅有截短的作用，而且能夠起到平滑的作用，在很多領域得到了應用。以下為利用窗函數(shù)法設計低通、高通和帶通濾波器。低通濾波器設計程序：x1,fs,bits=wavread(c:代碼wy.wav);wp1=2*pi*3800/fs;ws1=2*pi*4000/fs;wc1=(wp1+ws1)/2/pi;bt=ws1-wp1;n1=ceil(1.8*pi/bt);window=boxcar(n1+1); %長度為n1的矩形窗windowb1=fir1(n1,wc1,window);figure(1);freqz(b1,1,512);title(低通濾波器的頻率響應);如3.5圖所示圖3.5 低通濾波器頻率響應高通濾波器設計程序：x1,fs,bits=wavread(c:代碼wy.wav);wp1=2*pi*3800/fs;ws1=2*pi*4000/fs;wc1=(wp1+ws1)/2/pi;bt=ws1-wp1;n1=ceil(1.8*pi/bt);n2=n1+mod(n1+1,2);window=boxcar(n2); %長度為n1的矩形窗windowb1=fir1(n2-1,wc1,high,window);figure(3);freqz(b1,1,512);title(高通濾波器的頻率響應)如圖3.6所示圖3.6 高通濾波器頻率響應帶通濾波器程序：clearclcfs=44100;wp1=2*pi*3205/fs;wp2=2*pi*4410/fs; %帶通濾波器的通帶范圍ws1=2*pi*3000/fs;ws2=2*pi*4600/fs;bt=wp1-ws1;wc=(wp1+ws1)/2/pi,(wp2+ws2)/2/pi;n3=ceil(1.8*pi/bt);window=boxcar(n3); %長度為n1的矩形窗windowb1=fir1(n3-1,wc,window);figure(9);freqz(b1,1,512); %數(shù)字濾波器頻率響應title(帶通濾波器的頻率響應);如3.7圖所示圖3.7 帶通濾波器頻率響應3.3 用雙線性法設計iir數(shù)字濾波器采用非線性頻率壓縮方法，將整個頻率軸上的頻率范圍壓縮到-/t/t之間，再用z=est轉(zhuǎn)換到z平面上。也就是說，第一步先將整個s平面壓縮映射到s1平面的-/t/t一條橫帶里；第二步再通過標準變換關系z=es1t將橫帶轉(zhuǎn)換到z平面上去。這樣就使z平面與s平面對應的建立了單值關系，也就是說，消除了多值變換性就代表著消除了頻譜混疊現(xiàn)象，映射關系如圖3.8所示。圖3.8 雙線性變換的映射關系為了將s平面的整個虛軸j壓縮到s1平面j1軸上的-/t到/t段上，可以通過以下的正切變換實現(xiàn) (3.6)式中,t仍是采樣間隔。當1由-/t經(jīng)過0變化到/t時，由-經(jīng)過0變化到+，也即映射了整個j軸。將式（2-6）寫成 (3.7)將此關系解析延拓到整個s平面和s1平面，令j=s，j1=s1，則得 (3.8)再將s1平面通過以下標準變換關系映射到z平面 (3.9)從而得到s平面和z平面的單值映射關系為： (3.10) (3.11)式（3.10）與式（3.10）是s平面與z平面之間的單值映射關系，這種變換都是兩個線性函數(shù)之比，因此稱為雙線性變換式（3.6）與式（3.10的雙線性變換符合映射變換應滿足的兩點要求。首先,把z=ej，可得 （3.12）即s平面的虛軸映射到z平面的單位圓。其次，將s=+j代入式（2-12），得 （3.13）因此 （3.14）由此看出，當0時，|z|0時，|z|1。也就是說，s平面的左半平面映射到z平面的單位圓內(nèi)，s平面的右半平面映射到z平面的單位圓外，s平面的虛軸映射到z平面的單位圓上。因此，穩(wěn)定的模擬濾波器經(jīng)雙線性變換后所得的數(shù)字濾波器也一定是穩(wěn)定的。本文只對巴特沃斯濾波器進行設計。用雙線性變換法設計巴特沃斯低通濾波器的程序如下:x1,fs,bits=wavread(c:代碼wy.wav);wp=2*3800/fs;ws=2*4000/fs;rp=3;as=10;n,wc=buttord(wp,ws,rp,as,s);b,a=butter(n,wc,s);bz,az=bilinear(b,a,fs);nd,wdc=buttord(wp,ws,rp,as);b,a=butter(nd,wdc);freqz(b,a,1000,fs);圖3.9 巴特沃斯低通濾波器的頻率響應 巴特沃斯低通濾波器零極點圖第4章 語音信號的濾波4.1 語音信號的濾波測試本節(jié)是通過將采集的原始語音信號分別通過已設計好的低通、高通和帶通濾波器，并觀察和分析它們各自的頻譜圖，并利用函數(shù)sound函數(shù)播放濾波后的語音號，仔細聆聽其中聲音的不同。將信號通過4000hz的低通濾波器后，因為只剩下低頻段的信號了，我們可以明顯發(fā)現(xiàn)其聲音與原始信號的聲音相比聲音較低沉。其信號經(jīng)過低通濾波器程序：x1,fs,bits=wavread(c:代碼wy.wav);n=length(x1);wp1=2*pi*3800/fs;ws1=2*pi*4000/fs;wc1=(wp1+ws1)/2/pi;bt=ws1-wp1;n1=ceil(1.8*pi/bt);window=boxcar(n1+1); %長度為n1的矩形窗windowb1=fir1(n1,wc1,window);z1 = filter(b1,1, x1); %對信號進行低通濾波figure(9);ftemp=-n/2:n/2-1.*fs/n;plot(ftemp,fftshift(abs(fft(z1);title(信號經(jīng)過低通濾波器（頻域）);sound(z1,fs);wavwrite(z1,fs,c:代碼語音播放低通.wav);%將低通濾波后的音樂信號保存在語音播放中。如4.1圖所示。將信號通過4000hz的高通濾波器后，因為只剩下高頻段的信號了，我們可以明顯發(fā)現(xiàn)其聲音與原始信號的聲音相比聲音比較刺耳，讓人聽著感覺難受。其信號進過高通濾波器程序：x1,fs,bits=wavread(c:代碼wy.wav);n=length(x1);圖4.1 經(jīng)過低通濾波器的信號頻譜wp1=2*pi*3800/fs;ws1=2*pi*4000/fs;wc1=(wp1+ws1)/2/pi;bt=ws1-wp1;n1=ceil(1.8*pi/bt);n2=n1+mod(n1+1,2);window=boxcar(n2); %長度為n1的矩形窗windowb1=fir1(n2-1,wc1,high,window);z1 = filter(b1,1, x1); %對信號進行高通濾波ftemp=-n/2:n/2-1.*fs/n;plot(ftemp,fftshift(abs(fft(z1);title(信號經(jīng)過高通濾波器（頻域）);sound(z1,fs);wavwrite(z1,fs,c:代碼語音播放高通.wav);%如圖4.2所示 圖4.2 經(jīng)過高通濾波器的信號頻譜將信號通過下面的帶通濾波器后，其播放的聲音讓人聽不出有音樂的感覺了。其信號進過帶通濾波器程序：x1,fs,bits=wavread(c:代碼wy.wav);n=length(x1);wp1=2*pi*3205/fs;wp2=2*pi*4410/fs; %帶通濾波器的通帶范圍ws1=2*pi*3000/fs;ws2=2*pi*4600/fs;bt=wp1-ws1;wc=(wp1+ws1)/2/pi,(wp2+ws2)/2/pi;n3=ceil(1.8*pi/bt);window=boxcar(n3); %長度為n1的矩形窗windowb1=fir1(n3-1,wc,window);z1 = filter(b1,1, x1); %對信號進行高通濾波figure(9);ftemp=-n/2:n/2-1.*fs/n;plot(ftemp,fftshift(abs(fft(z1);title(信號經(jīng)過帶通濾波器（頻域）);sound(z1,fs);wavwrite(z1,fs,c:代碼語音播放帶通.wav);%將帶通濾波后的音樂信號保存在語音播放中如圖4.3所示 圖4.3 經(jīng)過帶通濾波器的信號頻譜將信號通過4000hz的巴特沃斯低通濾波器后，因為只剩下低頻段的信號了，我們可以明顯發(fā)現(xiàn)其聲音與原始信號的聲音相比聲音較低沉。其信號經(jīng)過低通濾波器程序：clcclearx1,fs,bits=wavread(c:代碼wy.wav);n=length(x1) %計算音頻信號的長度wp=2*3800/fs;ws=2*4000/fs;rp=3;as=10;n,wc=buttord(wp,ws,rp,as,s);b,a=butter(n,wc,s);bz,az=bilinear(b,a,fs);nd,wdc=buttord(wp,ws,rp,as);b,a=butter(nd,wdc);freqz(b,a,1000,fs);z2= filter(b,a,x1); ftemp=-n/2:n/2-1.*fs/n;figure(1);plot(ftemp,fftshift(abs(fft(z2);title(信號經(jīng)過低通濾波器（頻域）);grid on;如圖4.4所示：圖4.44.2 加噪語音信號的濾波處理在以上的程序中，我們已對原始信號進行了加噪處理(其加雜的是4000hz以后的高頻成分)，所以我們現(xiàn)在要想把語音信號進行還原，則必須要通過4000hz的低通濾波器對其進行濾波處理，然后我們就可以得到與原信號無限逼近的語音信號，通過以下程序，我們可以看到其頻譜圖與原信號頻譜圖幾乎一樣，聲音頻率也很相似。這樣，我們的加噪語音信號的濾波處理就實現(xiàn)了。%4000hz的低通濾波器對噪音信號進行濾波%wp1=2*pi*3800/fs;ws1=2*pi*4000/fs;wc1=(wp1+ws1)/2/pi;bt=ws1-wp1;n1=ceil(1.8*pi/bt);window=boxcar(n1+1); %長度為n1的矩形窗windowb1=fir1(n1,wc1,window);z2 = filter(b1,1, ifft(z1); %對信號進行低通濾波figure(9);ftemp=-n/2:n/2-1.*fs/n;plot(ftemp,fftshift(abs(fft(z2);title(信號經(jīng)過低通濾波器（頻域）);grid on;so