MATLAB數(shù)字信號處理課件

第5章使用MATLAB實現(xiàn)數(shù)字信號處理本章主要內(nèi)容如下：5．1 數(shù)字信號處理基本內(nèi)容及相應(yīng)的MATLAB工具5．2 信號通過系統(tǒng)的時域分析5．3 信號通過系統(tǒng)的頻域和Z域分析5．4 濾波器設(shè)計5．5 頻譜分析第5章使用MATLAB實現(xiàn)數(shù)字信號處理本章主要內(nèi)容如下：15.1數(shù)字信號處理基本內(nèi)容

及相應(yīng)的MATLAB工具數(shù)字信號處理的基本內(nèi)容通常分為兩部分：離散時間信號與系統(tǒng)分析主要涉及離散時間信號與系統(tǒng)的時域、頻域表示，以及信號通過系統(tǒng)的時域、頻域分析及其變換域分析。

MATLAB函數(shù)庫中提供了filter,conv,convmtx,fft,ifft,freqz,impz,zplane等與之相應(yīng)的函數(shù)。數(shù)字濾波器設(shè)計和譜分析數(shù)字濾波器設(shè)計包括了無限沖激響應(yīng)(IIR)和有限沖激響應(yīng)（FIR）濾波器設(shè)計，譜分析又可進(jìn)一步分為線性譜分析和非線性譜分析。MATLAB為此提供了多種成熟算法的相應(yīng)函數(shù)以及極為豐富的設(shè)計工具。5.1數(shù)字信號處理基本內(nèi)容

及相應(yīng)的MATLAB工具數(shù)字信25.2 時域分析卷積，濾波，單位沖激響應(yīng)5.2 時域分析卷積，濾波，單位沖激響應(yīng)35.2.1卷積MATLAB提供

conv函數(shù)實現(xiàn)標(biāo)準(zhǔn)的一維信號卷積:

例如，若系統(tǒng)h(n)為 >>h=[111]輸入序列x(n)為 >>x=[111]則x(n)經(jīng)過系統(tǒng)h(n)后的MATLAB實現(xiàn)為： >>conv(h,x)或

conv([111],[111])執(zhí)行后即得到y(tǒng)(n)為 ans= 12321注意：使用conv函數(shù)時，h(n)和x(n)都必須是有限長的,，否則不能使用conv函數(shù)。5.2.1卷積MATLAB提供conv函數(shù)實現(xiàn)標(biāo)準(zhǔn)的一維4例5-1時域離散序列的卷積計算

與圖形顯示

例5-1(教材p63)

：已知離散信號x(n)和h(n)，求y(n)=x(n)*h(n)，并用圖形表示。例5-1時域離散序列的卷積計算

與圖形顯示例5-1(教5例5-1的MATLAB程序Nh=20;Nx=10;m=5; %設(shè)定Nx，Nh和位移值mn=0:Nh-1;h1=(0.9).^n; %產(chǎn)生h1(n)h2=h1;nx=0:Nx-1;x1=ones(1,Nx); %產(chǎn)生x1(n)x2=zeros(1,Nx+m);fork=m+1:m+Nx %產(chǎn)生x2(n)=x1(n-m)x2(k)=x1(k-m);end %產(chǎn)生x2(n)y1=conv(x1,h1); %計算y1(n)=x1(n)*h1(n)y2=conv(x2,h2); %計算y2(n)=x2(n)*h2(n)subplot(3,2,1)stem(nx,x1,'.')axis([03001.2]),title(‘x1(n)’) %繪圖……（以下省略）例5-1的MATLAB程序Nh=20;Nx=10;m=5; 65.2.2濾波數(shù)字濾波器的系統(tǒng)函數(shù)H(z)用如下式表示:在MATLAB中，用向量b，a來表示濾波器的系數(shù)b(i)和

a(i)。

5.2.2濾波數(shù)字濾波器的系統(tǒng)函數(shù)H(z)用如下式表示:7濾波器分類當(dāng)n

=

0，m≠0時，稱為AR濾波器，即自回歸（AutoRecurrence）濾波器，具無限沖激響應(yīng)（IIR），也即其單位采樣響應(yīng)h(n)具無限長度；若m

=

0，a(1)≠0，稱為MA濾波器，即滑動平均（MovingAverage）濾波器，其單位采樣響應(yīng)h(n)是有限長度，故稱有限沖激響應(yīng)（FIR）濾波器；如果n、m都大于零，稱為ARMA濾波器，而其沖激響應(yīng)也為IIR。濾波器分類當(dāng)n

=

0，m≠0時，稱為AR濾波器，即自回歸（8filter函數(shù)

MATLAB提供了filter函數(shù)來對離散信號進(jìn)行濾波，表達(dá)信號通過系統(tǒng)后的結(jié)果。與conv不同的是，filter函數(shù)可適用于無限沖激響應(yīng)系統(tǒng)的情況，但信號仍須是有限長的。例如，一個單極點的低通濾波器系數(shù)如下：

>>b=1; %分子系數(shù)向量b(i) >>a=[1-0.9]; %分母系數(shù)向量a(i)

如果用filter函數(shù)實現(xiàn)對信號x濾波，只要調(diào)用： >>y=filter(b,a,x);

就可給出輸入x經(jīng)過濾波以后的輸出y。filter函數(shù)MATLAB提供了filter函數(shù)來對離95.2.3單位沖激響應(yīng)數(shù)字濾波器的單位沖激響應(yīng)定義為輸入為單位樣本序列時數(shù)字濾波器的響應(yīng)，即：h(n)=T[δ(n)] 其中：5.2.3單位沖激響應(yīng)數(shù)字濾波器的單位沖激響應(yīng)定義為輸入10單位沖激響應(yīng)的MATLAB實現(xiàn)MATLAB近似實現(xiàn)單位采樣信號的方法為：

imp=[1;zeros(p,1)]; %zeros(p,1)產(chǎn)生p個零元素組成的列向量，p是正整數(shù)。使用imp后，濾波器的沖激響應(yīng)可近似得到為：

h=filter(b,a,imp);impz函數(shù)可以直接求出數(shù)字濾波器的單位沖激響應(yīng)，即：impz(b,a) 該命令將同時繪出濾波器的單位沖激響應(yīng)，教材p66圖5-2。單位沖激響應(yīng)的MATLAB實現(xiàn)MATLAB近似實現(xiàn)單位采樣信115．3頻域和Z域分析

頻率響應(yīng)，零極點分析

5．3頻域和Z域分析頻率響應(yīng)，零極點分析125.3.1頻率響應(yīng)MATLAB數(shù)字信號處理工具箱有很多函數(shù)提供對模擬和數(shù)字濾波器的頻率響應(yīng)分析。其中，freqz函數(shù)和freqs函數(shù)分別返回數(shù)字和模擬濾波器的頻率響應(yīng)。工具箱中通常使用的單位頻率是Nyquist頻率，即采樣頻率的1/2。注意：就數(shù)字濾波器函數(shù)來說，其頻域指標(biāo)中的所有頻率都以Nyquist頻率進(jìn)行歸一化。因此Nyquist頻率也稱歸一化頻率。5.3.1頻率響應(yīng)MATLAB數(shù)字信號處理工具箱有很多函13關(guān)于Nyquist頻率的說明例如：系統(tǒng)采樣頻率為1000Hz，則若數(shù)字濾波器的截止頻率等于300Hz，經(jīng)Nyquist頻率歸一化后，其歸一化頻率就是300/500=0.6。若將歸一化頻率轉(zhuǎn)換成數(shù)字信號處理教科書中所使用的數(shù)字頻率(rad)，需乘以π；反之，若乘以采樣頻率fs的一半，則將歸一化頻率轉(zhuǎn)換回了模擬域頻率(Hz)。歸一化頻率f應(yīng)滿足0<f<1。對于頻率響應(yīng)而言，歸一化頻率和模擬頻率都可使用。關(guān)于Nyquist頻率的說明例如：系統(tǒng)采樣頻率為1000Hz14freqz命令

freqz

函數(shù)提供了基于FFT算法所得到的數(shù)字濾波器H(z)的頻率響應(yīng)最常用的形式：[h,w]=freqz(b,a,p)

其意義是，對于數(shù)字濾波器： freqz函數(shù)接受濾波器系數(shù)向量b和a，以及一個用于指定所需計算的頻率響應(yīng)樣點數(shù)的整數(shù)p，返回一個與由p個頻率樣本點構(gòu)成的實頻率向量相對應(yīng)的復(fù)頻率響應(yīng)向量h。

freqz命令

freqz函數(shù)提供了基于FFT算法所得到的15freqz的命令形式h=freqz(b,a,w)采用上面的形式時，需先對頻率樣本點向量w作出定義。通常的做法是使用linspace函數(shù)。[h,w]=freqz(b,a)w和p沒有定義，默認(rèn)w由（0~π）上均分的512點構(gòu)成，頻率單位為rad/sample。[h,w]=freqz(b,a,p，’whole’)使用帶參數(shù)’whole’選項的命令形式[h,w]=freqz(b,a,p,fs)用參數(shù)選項指定采樣頻率fsfreqz的命令形式h=freqz(b,a,w)16linspace的命令形式linspace函數(shù)的命令形式有兩種：linspace（x1,x2） %產(chǎn)生[x1,x2]范圍內(nèi)經(jīng)線性分割得到的100點的向量linspace（x1,x2，N） %產(chǎn)生[x1,x2]范圍內(nèi)經(jīng)線性分割得到的N點的向量

例如：w=linspace（0，pi）%

w由（0，pi）上100個等分點組成 w=linspace（0，1000）

%w由（0，1000）上100個等分點組成linspace的命令形式linspace函數(shù)的命令形式有兩17freqz函數(shù)無返回輸出參數(shù)格式的調(diào)用freqz函數(shù)還能以無返回輸出參數(shù)的格式調(diào)用： freqz（b，a，256）或freqz（b，a，256，2000）無返回輸出參數(shù)調(diào)用freqz函數(shù)的好處是，MATLAB能夠自動繪出頻率在（0~π）范圍內(nèi)的幅頻特性和相頻特性圖。注意，無返回輸出參數(shù)調(diào)用freqz函數(shù)繪出的相頻特性不能正確給定在處的值，這是因為在MATLAB中，屬于下半個單位圓。

freqz函數(shù)無返回輸出參數(shù)格式的調(diào)用freqz函數(shù)還能以無18頻率響應(yīng)的實例例：先構(gòu)成一個截止頻率為400Hz的9階巴特沃思（Butterworth）低通數(shù)字濾波器，求出其系數(shù)b，a，再求出其256點頻率響應(yīng)。指定的采樣頻率fs=2000Hz。實現(xiàn)1：先調(diào)用butter函數(shù)，再調(diào)用freqz函數(shù)；實現(xiàn)2：無返回輸出參數(shù)，調(diào)用freqz函數(shù)；實現(xiàn)3：為得到完整的（0~π）范圍內(nèi)的幅相特性圖，可用帶返回輸出參數(shù)的格式調(diào)用freqz，在得到各頻率點上的頻率響應(yīng)的數(shù)值后，使用MATLAB所提供的abs函數(shù)和angle函數(shù)分別提取幅頻特性與相頻特性。注意：為了得到連續(xù)的相頻曲線，需要使用unwrap函數(shù)，該函數(shù)通過在發(fā)生跳變后的各處自動加上，從而除去了相位的跳變，這稱為相位的“解卷繞”。

頻率響應(yīng)的實例例：先構(gòu)成一個截止頻率為400Hz的9階巴特沃19實現(xiàn)1

[b,a]=butter(9,400/1000); %括號中第一個輸入是階數(shù)，第二個是歸一化截止頻率

[h,f]=freqz(b,a,256,2000);或

h=freqz（b，a，256，2000）；實現(xiàn)1 [b,a]=butter(9,400/1000)20實現(xiàn)2實現(xiàn)1中第二條命令的形式可改寫為freqz（b，a，256）

或freqz（b，a，256，2000）可自動繪出頻率在（0~π）范圍內(nèi)的幅頻特性和相頻特性圖。注意無返回輸出參數(shù)調(diào)用freqz函數(shù)繪出的相頻特性不能正確給定在ω=π處的值，因為ω=π屬于下半個單位圓。實現(xiàn)2實現(xiàn)1中第二條命令的形式可改寫為21實現(xiàn)3[b,a]=butter(9,400/1000);w=linspace(0,1000);h=freqz(b,a,w,2000);mag=abs(h);pha=angle(h);%提取濾波器頻率響應(yīng)的幅度mag和相位phasemilogy(w,mag);title('Magnitude');figure;plot(w,pha);title('Phase');出現(xiàn)了幅度為2π的跳變——卷繞實現(xiàn)3[b,a]=butter(9,400/1000);22解卷繞：unwrap函數(shù)為了得到連續(xù)的相頻曲線，需要使用unwrap函數(shù)，該函數(shù)通過在發(fā)生2π跳變后的各處自動加上2π，從而除去了相位的跳變，這稱為相位的“解卷繞”。解卷繞解卷繞：unwrap函數(shù)為了得到連續(xù)的相頻曲線，需要使用u23freqs函數(shù)freqs函數(shù)用于求取由b和a系數(shù)向量定義的模擬濾波器 H（s）的頻率響應(yīng)。使用方法類似freqz函數(shù)。與第二章（p32例2-1）采用數(shù)組相除方法求取頻率響應(yīng)相比，使用freqs函數(shù)要方便很多。freqs函數(shù)freqs函數(shù)用于求取由b和a系數(shù)向量定義的模245.3.2零極點分析zplane函數(shù)用于畫出線性系統(tǒng)在Z平面上的零極點。有兩種使用方法：

1、在已知零極點時，例如某濾波器的零點為-1/2，一對共軛極點為和時，只要輸入命令zer=-0.5;pol=0.9*exp(j*2*pi*[-0.30.3]');zplane(zer,pol)

即可畫出零極點。（見p70圖5-6）5.3.2零極點分析zplane函數(shù)用于畫出線性系統(tǒng)25用zplane函數(shù)繪制零極點圖另一種情況：已知系統(tǒng)的系統(tǒng)函數(shù)系數(shù)向量b和a，則可通過調(diào)用zplane(b,a)繪出零極點。這種情形下，zplane函數(shù)先求得系統(tǒng)函數(shù)的零點和極點，然后繪出零極點圖。用zplane函數(shù)繪制零極點圖另一種情況：已知系統(tǒng)的系統(tǒng)函數(shù)265.4 濾波器設(shè)計濾波器技術(shù)指標(biāo)，IIR濾波器設(shè)計，F(xiàn)IR濾波器設(shè)計5.4 濾波器設(shè)計濾波器技術(shù)指標(biāo)，275.4.1濾波器的技術(shù)指標(biāo)設(shè)計濾波器之前，要定義一組濾波器的技術(shù)指標(biāo)。濾波器的技術(shù)指標(biāo)通常采用幅度指標(biāo)給出，以保證物理上的可實現(xiàn)性。幅度指標(biāo)以容差方式給出，如p71圖5-7所示。具體指標(biāo)：通帶截止頻率Wp；阻帶截止頻率Ws；通帶波紋δp；阻帶衰減δs。其中通帶波紋和阻帶衰減也可以用分貝數(shù)給出5.4.1濾波器的技術(shù)指標(biāo)設(shè)計濾波器之前，要定義一組濾波28濾波器以容差方式給出的幅度指標(biāo)

|H(jw)|10通帶過渡帶阻帶pwswsdpd-1w濾波器以容差方式給出的幅度指標(biāo)|H(jw)|10通帶過渡295.4.2IIR濾波器設(shè)計MATLAB工具箱提供了幾種模擬濾波器原型的產(chǎn)生函數(shù),如巴特沃思(Butterworth),切比雪夫(Chebyshev)濾波器等。還提供了從模擬低通濾波器原型轉(zhuǎn)換為高通、帶通和帶阻的轉(zhuǎn)換函數(shù)，模擬濾波器轉(zhuǎn)換為數(shù)字濾波器的雙線性變換法和沖激響應(yīng)不變法，模擬IIR濾波器的階數(shù)選擇函數(shù)以及數(shù)字濾波器直接設(shè)計函數(shù)等等，使用起來非常方便。本節(jié)以巴特沃思濾波器為例

介紹MATLAB工具箱提供的濾波器設(shè)計函數(shù)使用方法。5.4.2IIR濾波器設(shè)計MATLAB工具箱提供了幾301．巴特沃思濾波器設(shè)計函數(shù)

[b,a]

=

butter(n,Wn,options)缺省情況下，返回一個用有理分式表示的低通數(shù)字濾波器系統(tǒng)函數(shù)的分子分母系數(shù)向量b和a。

設(shè)計的技術(shù)指標(biāo)只需要指定一個歸一化截止頻率Wn（這里Wn的單位為：Nyquist

頻率

=

1Hz）。參數(shù)選項options：選項空缺返回低通濾波器；‘high’則表示返回的是高通數(shù)字濾波器；至于帶通和帶阻濾波器設(shè)計，則輸入的Wn必須是二個元素的向量；選項options空缺返回帶通濾波器，帶阻濾波器需加參數(shù)options為'stop'。

1．巴特沃思濾波器設(shè)計函數(shù)[b,a]

=

butter(n31Butter函數(shù)舉例[B,A]=butter(N,Wn)返回系統(tǒng)函數(shù)的長度為N+1的分子系數(shù)B和分母系數(shù)A。式中Wn為歸一化截止頻率，0<Wn<1.0。[B,A]=butter(N,Wn,'high')設(shè)計一個高通數(shù)字濾波器。如果Wn=[W1W2]，且options空缺，則返回一個2N階的帶通濾波器，其通帶為：W1<W<W2。[B,A]=butter(N,Wn,'stop')，且Wn=[W1W2]設(shè)計2N階的帶阻濾波器。butter(N,Wn,'s')，butter(N,Wn,'high','s')或

butter(N,Wn,'stop','s')調(diào)用butter()函數(shù)附加一個options為‘s’時，表示設(shè)計的是一個模擬濾波器，而指定的截止頻率Wn單位此時必須為rad/s。Butter函數(shù)舉例[B,A]=butter(N,Wn322.巴特沃思濾波器階次選擇函數(shù)buttord函數(shù)用于給出滿足給定頻域指標(biāo)要求的最小階次濾波器的設(shè)計。其調(diào)用格式為：

[n,Wn]

=

buttord(Wp,Ws,Rp,Rs) %用于數(shù)字濾波器[n,Wn]

=

buttord(Wp,Ws,Rp,Rs,'s')%用于模擬濾波器該函數(shù)返回符合技術(shù)指標(biāo)的最小階數(shù)N以及Butterworth濾波器3dB截止頻率Wn。設(shè)計指標(biāo)同前：通帶起伏不超過Rp，阻帶衰減不小于Rs，Wp和Ws分別為歸一化通帶和阻帶截止頻率。

2.巴特沃思濾波器階次選擇函數(shù)buttord函數(shù)用于給出滿足33Buttord函數(shù)應(yīng)用：例5-1

例：設(shè)計一個數(shù)字帶通濾波器，通帶為1000到2000Hz，止帶的起始位置離開通帶兩邊500Hz，也即分別為500Hz和2500Hz，設(shè)采樣率為10KHz，通帶起伏為1dB,止帶最小衰減為60dB。解MATLAB實現(xiàn)如下：

[n,Wn]=buttord([10002000]/5000,[5002500]/5000,1,60)得到： n= 12 Wn= 0.19510.4080因此 [b,a]=butter(n,Wn);給出了滿足要求的最小階次巴特沃思帶通濾波器。Buttord函數(shù)應(yīng)用：例5-1例：設(shè)計一個數(shù)字帶通濾波器343．經(jīng)典的IIR濾波器設(shè)計方法一種常用的IIR濾波器設(shè)計方法是利用已經(jīng)很成熟的模擬濾波器設(shè)計結(jié)果。設(shè)計步驟如下：（1）按照一定規(guī)則將給出的數(shù)字濾波器的技術(shù)指標(biāo)轉(zhuǎn)換為模擬低通濾波器的技術(shù)指標(biāo)。 （2）根據(jù)轉(zhuǎn)換后的技術(shù)指標(biāo)設(shè)計模擬低通濾波器的H(s)。 （3）再按一定規(guī)則將H(s)轉(zhuǎn)換成H(z)。步驟(3)中最常用的方法:沖激響應(yīng)不變法（也克稱為脈沖響應(yīng)不變法）和雙線性變換法。

3．經(jīng)典的IIR濾波器設(shè)計方法一種常用的IIR濾波器設(shè)計方法35沖激響應(yīng)不變法的原理

沖激響應(yīng)不變法是通過對模擬濾波器的單位沖激響應(yīng)ha(nTs),進(jìn)行采樣，將其作為數(shù)字濾波器的單位沖激響應(yīng)h(n)。即：據(jù)此設(shè)計得到的數(shù)字濾波器的頻率響應(yīng)與模擬濾波器的頻率響應(yīng)之間的關(guān)系為：因此存在著頻率混迭。從而限制了其適用范圍。注意：S域頻率軸與Z域單位圓之間不存在一一對應(yīng)關(guān)系。沖激響應(yīng)不變法的原理沖激響應(yīng)不變法是通過對模擬濾波器的單位36沖激響應(yīng)不變設(shè)計的impinvar()函數(shù)impinvar()函數(shù)的調(diào)用格式是：

[Bz,Az]=impinvar(Bs,As,fs)這一函數(shù)把具有[Bs,As]模擬濾波器的轉(zhuǎn)換成了采樣頻率為fs的數(shù)字濾波器的[Bz,Az]。如果沒有給定采樣頻率fs，函數(shù)默認(rèn)為1Hz。

沖激響應(yīng)不變設(shè)計的impinvar()函數(shù)impinvar(37雙線性變換法原理雙線性變換法所定義的s平面與z平面的映射關(guān)系為：給定模擬濾波器的系統(tǒng)函數(shù)，數(shù)字濾波器為：因此S域內(nèi)的有理分式函數(shù)在Z域內(nèi)仍然保持了有理分式函數(shù)形式，從而保證了數(shù)字濾波器的物理可實現(xiàn)性。雙線性變換法原理雙線性變換法所定義的s平面與z平面的映射關(guān)系38雙線性變換法原理（續(xù)）雙線性變換將S平面上的頻率軸一對一映射到了Z平面的單位圓上，消除了頻率混疊問題：但Z平面單位圓和S域頻率軸之間存在嚴(yán)重的非線性關(guān)系；由此也對其適用范圍產(chǎn)生了限制；并且，在濾波器設(shè)計中，在頻域指標(biāo)轉(zhuǎn)換時需要進(jìn)行Pre-warping。

雙線性變換法原理（續(xù)）雙線性變換將S平面上的頻率軸一對一映射39雙線性變換設(shè)計的bilinear函數(shù)bilinear()函數(shù)的調(diào)用格式有三種：

[BzAz]=bilinear(Bs,As,Fs)%把具有[Bs,As]模擬濾波器的轉(zhuǎn)換成數(shù)字濾波器的[Bz,Az]。Fs是采樣頻率。

[ZdPdKd]=bilinear(Z,P,K,Fs)%把模擬濾波器的零點、極點，轉(zhuǎn)換成數(shù)字濾波器的零點、極點模型，F(xiàn)s是采樣頻率。[AdBdCdDd]=bilinear(A,B,C,D,Fs)%把模擬濾波器的狀態(tài)方程模型轉(zhuǎn)換成數(shù)字濾波器的狀態(tài)方程模型。Fs是采樣頻率。注意：這里的采樣頻率均為實際頻率。雙線性變換設(shè)計的bilinear函數(shù)bilinear()函數(shù)40例5-3IIR數(shù)字濾波器的設(shè)計示例

例5-3中直接給出了模擬低通濾波器技術(shù)指標(biāo)（）。如果給定的是數(shù)字低通濾波器技術(shù)指標(biāo)，則應(yīng)先將數(shù)字濾波器的技術(shù)指標(biāo)折合成相應(yīng)的模擬低通濾波器指標(biāo)，然后采用示例所用的設(shè)計方法進(jìn)行濾波器設(shè)計。

例5-3IIR數(shù)字濾波器的設(shè)計示例例5-3中直接給出了41例5-3設(shè)計要求與指標(biāo)（1）設(shè)計一個巴特沃思模擬低通濾波器，指標(biāo)如下： 通帶內(nèi) ，波紋小于10dB； 阻帶內(nèi)，衰減不小于20dB。（2）用雙線性變換法，將上述模擬低通濾波器變換為數(shù)字低通濾波器，采樣間隔分別取T=1s，1.5s，比較設(shè)計結(jié)果。（3）用沖激響應(yīng)不變法，將上述的模擬低通濾波器變換為數(shù)字低通濾波器，采樣間隔分別取T=1s，1.5s，比較設(shè)計結(jié)果。

例5-3設(shè)計要求與指標(biāo)（1）設(shè)計一個巴特沃思模擬低通濾波42例5-3的解由于給定的設(shè)計指標(biāo)是模擬濾波器的頻域指標(biāo)，所以，首先選用buttord()函數(shù)進(jìn)行最小階數(shù)的濾波器設(shè)計。 [N,Wn]

=

buttord(Wp,Ws,Rp,Rs,'s') 注意這里Wp，Ws單位為rad/s。此題中，Wp=0.2*pi,Ws=0.3*pi,Rp=10,Rs=20，所以可用 [N,Wn]=buttord(0.2*pi,0.3*pi,10,20,'s');得到實現(xiàn)該頻域指標(biāo)要求的巴特沃思低通濾波器最低階次N和3dB截止頻率Wn。例5-3的解由于給定的設(shè)計指標(biāo)是模擬濾波器的頻域指標(biāo)，所以，43例5-3的解（續(xù)）再使用butter()函數(shù)，得到模擬低通濾波器的系統(tǒng)函數(shù)的系數(shù)向量Bs和As： [Bs,As]=butter(N,Wn,'s');根據(jù)采樣間隔分別為T=1s，1.5s，可得采樣頻率分別為：Fs1=1Hz，F(xiàn)s2=2/3Hz，調(diào)用bilinear()函數(shù)和impinvar()函數(shù)，可以得到相應(yīng)的兩種數(shù)字濾波器系統(tǒng)函數(shù)的系數(shù)向量Bz和Az。即： [Bz,Az]=bilinear(Bs，As，F(xiàn)s) 或 [Bz,Az]=impinvar(Bs，As，F(xiàn)s)注意，由于低通巴特沃思濾波器并不是帶限的，因此對其用沖激響應(yīng)不變法設(shè)計數(shù)字濾波器時，將會發(fā)生混疊。（參見下張ppt）

例5-3的解（續(xù)）再使用butter()函數(shù)，得到模擬低通44出現(xiàn)混疊無混疊出現(xiàn)混疊無混疊455.4.3FIR濾波器設(shè)計FIR濾波器具有IIR濾波器所沒有的線性相位特性，所以在工程中得到了廣泛的應(yīng)用。

MATLAB中除了cremez函數(shù)以外，所有用于FIR濾波器設(shè)計的函數(shù)都是用于設(shè)計線性相位的FIR濾波器的。5.4.3FIR濾波器設(shè)計FIR濾波器具有IIR濾波器所46FIR濾波器的線性相位條件

FIR數(shù)字濾波器的系統(tǒng)函數(shù)為：線性相位條件：FIR數(shù)字濾波器的h(n)為實數(shù)，且滿足以下任一條件

h(n)=h(N-1-n)（關(guān)于n=(N-1)/2偶對稱）h(n)＝-h(N-1-n)（關(guān)于n=(N-1)/2奇對稱）FIR濾波器的線性相位條件FIR數(shù)字濾波器的系統(tǒng)函數(shù)為：47根據(jù)h(n)的對稱情況，分四種類型h(n)奇對稱，N為奇數(shù)

h(n)奇對稱，N為偶數(shù)

h(n)偶對稱，N為奇數(shù)

h(n)偶對稱，N為偶數(shù)

根據(jù)h(n)的對稱情況，分四種類型h(n)奇對稱，N為奇數(shù)48FIR線性相位濾波器的頻率響應(yīng)具有線性相位的FIR濾波器的頻率響應(yīng)函數(shù)可以用下式表示：式中或，。注意，式中N為脈沖響應(yīng)h(n)的長度，是振幅（Amplitude）響應(yīng)函數(shù)，不是幅度（Magnitude）響應(yīng)函數(shù)，可以有正有負(fù)，與其相關(guān)的相位響應(yīng)也是一個連續(xù)線性函數(shù)。對于前述四種h(n),其幅度函數(shù)有不同的特點。FIR線性相位濾波器的頻率響應(yīng)具有線性相位的FIR濾波器的頻49FIR濾波器設(shè)計：窗函數(shù)法基本原理：將理想低通濾波器的無限長沖激響應(yīng)序列截斷（等效于加矩形窗）后乘以窗函數(shù)，用得到的有限長序列去逼近理想低通濾波器。設(shè)所要求的理想數(shù)字濾波器的頻率響應(yīng)為, 是與其對應(yīng)的單位脈沖響應(yīng)，因此

由于是矩形頻率特性，故一定是無限長的非因果序列。而所要設(shè)計的是有限長的FIR濾波器，因此用有限長的去逼近無限長的最直接的方法就是截斷，即用有限長的窗函數(shù)來截取得到：

這就是窗函數(shù)設(shè)計方法。

FIR濾波器設(shè)計：窗函數(shù)法基本原理：將理想低通濾波器的無限長50窗函數(shù)法設(shè)計步驟（1）設(shè)所要求的FIR濾波器理想頻率響應(yīng)為：（2）由性能指標(biāo)確定窗函數(shù)W(n)和窗口長度N，由過渡帶寬近似于窗口函數(shù)主瓣寬求得窗口長度N；（3）由求得實際濾波器的單位脈沖響應(yīng)h(n)，即求取所設(shè)計的FIR濾波器h(n)的系數(shù)向量；（4）檢驗濾波器的性能。窗函數(shù)法設(shè)計步驟（1）設(shè)所要求的FIR濾波器理想頻率響應(yīng)為：51FIR濾波器設(shè)計：fir1函數(shù)MATLAB中使用窗函數(shù)法設(shè)計線性相位FIR濾波器的常用設(shè)計函數(shù)為fir1()，調(diào)用格式如下：b=fir1(n,Wn)b=fir1(n,Wn,'ftype')b=fir1(n,Wn,window)b=fir1(n,Wn,'ftype',window)FIR濾波器設(shè)計：fir1函數(shù)MATLAB中使用窗函數(shù)法設(shè)計52fir1命令格式的說明fir1命令中n為FIR濾波器的階數(shù)，對于高通、帶阻濾波器，n應(yīng)該取偶數(shù)；Wn為濾波器的截止頻率，0<Wn<1,對于帶通、帶阻濾波器，Wn=[W1W2]，且滿足W1<W2,對于多帶濾波器，Wn=[W1W2W3…]，且滿足0<W1<W2<…<1；b為FIR濾波器的系數(shù)向量，長度為n+1；‘ftype’是濾波器的類型，默認(rèn)時為低通或帶通濾波器，‘high’為高通濾波器，‘stop’為帶阻濾波器；Window為窗函數(shù)，是長度為n+1的列向量。fir1命令格式的說明fir1命令中n為FIR濾波器的階數(shù)，53例5-4：低通濾波器例5-4理想數(shù)字低通濾波器的截止頻率為ω0

（rad），在頻率小于ω0處，其幅度為1，在頻率為ω0至π處，幅度為0，故它的沖激響應(yīng)序列h(n)為：其沖激響應(yīng)為無限長，因此不是因果可實現(xiàn)的。用窗截取其有限長度，再加以因果化，就可以得到一個可實現(xiàn)的線性相位FIR濾波器。例5-4：低通濾波器例5-4理想數(shù)字低通濾波器的截止頻率54例5-4的解例如取長度51，低通截止頻率ω0=

0.4π，則得到:b=0.4*sinc(0.4*(-25:25));這里用了最簡單的矩形窗。下面的命令將顯示這個濾波器的頻率響應(yīng)（見p77圖5-9）：

[H,w]=freqz(b,1,512,2); %b為H(z)分子系數(shù)，1為分母系數(shù)，512個樣本，采樣頻率2Hz

plot(w,abs(H));grid

title('TruncatedSincLowpassFIRFilter');例5-4的解例如取長度51，低通截止頻率ω0=0.4π，則55圖5-9窗函數(shù)法低通濾波器的頻率特性

圖5-9窗函數(shù)法低通濾波器的頻率特性56使用其它窗函數(shù)可以使用其他形式的窗函數(shù)；例如改用長度51的海明窗，即對H(z)的系數(shù)乘以海明窗函數(shù)，則命令為：b=0.4*sinc(0.4*(-25:25));b=b.*hamming(51)';%對H(z)的系數(shù)乘以海明窗函數(shù)[H,w]=freqz(b,1,512,2);plot(w,abs(H)),gridplot(w,20*log10(abs(H))),title('Hamming-WindowedTruncatedSincLPFIRFilter');下一張圖將給出其幅頻特性，注意與矩形窗情況比較。使用其它窗函數(shù)可以使用其他形式的窗函數(shù)；57圖5-10使用海明窗得到的低通濾波器特性圖5-10使用海明窗得到的低通濾波器特性58直接使用FIR濾波器設(shè)計函數(shù)fir1()

b=fir1(51,0.4);%fir1()中窗函數(shù)（window）項空缺代表使用海明窗[H,W]=freqz(b,1,512);plot(W/pi,abs(H));gridtitle(‘Hamming_windowedFIRLPF’)直接使用FIR濾波器設(shè)計函數(shù)fir1()b=fir1(559fir2()函數(shù)函數(shù)fir2()也是FIR濾波器設(shè)計函數(shù)，用于設(shè)計具有任意形狀頻率響應(yīng)的FIR數(shù)字濾波器。例如： n=50;f=[00.40.51];%給定通帶、止帶范圍分別為[00.4]和[0.51]，1為Nyquist頻率m=[1100];%給定通帶，止帶頻率響應(yīng)幅度要求b=fir2(n,f,m) 將返回行向量b為包含n+1個系數(shù)的n階FIR濾波器，其頻率響應(yīng)幅度符合給定的f和m要求。類似于fir1()函數(shù)，可以選擇窗，命令形式為

b=fir2(n,f,m,window)fir2()函數(shù)函數(shù)fir2()也是FIR濾波器設(shè)計函60FIR濾波器設(shè)計：頻率采樣設(shè)計法頻率采樣設(shè)計法是從頻域出發(fā)，把給定的理想頻率響應(yīng)等間隔采樣N個點得到。再對作IDFT，得到然后將h（n）作為所設(shè)計的濾波器的單位沖激響應(yīng)。

FIR濾波器設(shè)計：頻率采樣設(shè)計法頻率采樣設(shè)計法是從頻域出發(fā)，61頻率采樣設(shè)計法的有關(guān)說明對進(jìn)行頻率采樣得到的H(k)，可以表示成：為保證設(shè)計出的FIR濾波器具有線性相位，h(n)應(yīng)在[0,N-1]上具有對稱性。由此進(jìn)行的理論分析表明，頻率采樣值的幅度與相位也須滿足一定的約束條件，分為4種類型。因此，頻率采樣法設(shè)計FIR濾波器可歸結(jié)為從滿足約束條件的H(k)中求出h(n)問題。

頻率采樣設(shè)計法的有關(guān)說明對進(jìn)行頻率采62頻率采樣設(shè)計法的MATLAB實現(xiàn)從滿足約束條件的H(k)中求出h(n)，可利用MATLAB中的ifft函數(shù)。注意：從時域角度看,求出的h(n)與欲設(shè)計的hd(n)是有差別的，兩者關(guān)系如下：從頻域角度看，由頻率采樣設(shè)計所得到的濾波器的頻率特性是由H(k)內(nèi)插形成的，因此對的逼近程度與特性的平滑性有關(guān)。

存在時域混疊

頻率采樣設(shè)計法的MATLAB實現(xiàn)從滿足約束條件的H(k)中求63例5-5：頻率采樣法設(shè)計舉例例5-5所要設(shè)計的FIR低通數(shù)字濾波器技術(shù)指標(biāo)為：，，例5-5：頻率采樣法設(shè)計舉例例5-5所要設(shè)計的FIR低通64例5-5的解在給定的圖5-12理想低通濾波器特性上選擇N=21點進(jìn)行采樣得到|H(k)|，如圖5-13。附加相位約束條件對H(k)求反變換IDFT就可得h(n)。例5-5的解在給定的圖5-12理想低通濾波器特性上選擇N=265圖5-14例5-5頻率采樣法設(shè)計低通濾波器(N=21)

最小阻帶衰減僅為14dB圖5-14例5-5頻率采樣法設(shè)計低通濾波器(N=21)66例5-5的改進(jìn)設(shè)計分析：從圖5-14幅度響應(yīng)曲線可見，N=21時，最小的阻帶衰減為14dB，這在絕大多數(shù)情況下是不能接受的。但僅僅增大采樣點數(shù)，比如N=61時，最小的阻帶衰減仍為18dB，仍然不能滿足通常的止帶指標(biāo)要求。為此，除了增大N以外，還應(yīng)采用在過渡帶插入過渡采樣點的辦法來除去原始的特性中的突變。

例5-5的改進(jìn)設(shè)計分析：從圖5-14幅度響應(yīng)曲線可見，N=267圖5-16例5-5的改進(jìn)設(shè)計

(N=61，增加兩個過渡采樣點)增加兩個過渡采樣點后，在N=61時，最小的阻帶衰減為45dB，已滿足通常的止帶衰減要求圖5-16例5-5的改進(jìn)設(shè)計

(N=61，增加兩個過渡采685.5頻譜分析頻譜分析的主要作用是對信號的頻率構(gòu)成進(jìn)行分析，了解其在頻帶內(nèi)的分布情況；本課程內(nèi)僅對確定性信號進(jìn)行，不涉及隨機信號的頻譜分析；因此主要工具為FFT，不涉及非線性譜分析或即參量譜分析的內(nèi)容。

5.5頻譜分析頻譜分析的主要作用是對信號的頻率構(gòu)成進(jìn)行69離散傅里葉變換（DFT）在離散信號與系統(tǒng)分析中，DFT有著頻譜分析儀的特性。對連續(xù)信號與系統(tǒng)通過時域采樣，也可以使用DFT進(jìn)行譜分析。因此，DFT在頻譜分析中起了核心的作用。另一方面，從技術(shù)上而言，DFT運算存在著高效的快速算法FFT，因此信號的頻譜分析中將主要使用MATLAB的fft函數(shù)以及ifft函數(shù)。

離散傅里葉變換（DFT）在離散信號與系統(tǒng)分析中，DFT有著頻70使用DFT對連續(xù)信號進(jìn)行譜分析工程實際應(yīng)用中，常會遇到連續(xù)時間信號，其頻譜也是連續(xù)函數(shù)。為了利用DFT對信號進(jìn)行頻譜分析，須先對進(jìn)行時域采樣，得到，再對進(jìn)行DFT，得到的則是的傅里葉變換在頻域區(qū)間上的N點等間隔采樣。對帶限信號來說，X(k)與在范圍內(nèi)上的采樣值除相差一個因子外完全相同。以上說明構(gòu)成了用DFT對連續(xù)信號進(jìn)行譜分析的理論基礎(chǔ)。使用DFT對連續(xù)信號進(jìn)行譜分析工程實際應(yīng)用中，常會遇到連續(xù)時71用DFT進(jìn)行連續(xù)信號頻譜分析時注意點

（1）若信號持續(xù)時間有限長，則其頻譜無限寬。（2）若信號頻譜有限，則其信號持續(xù)時間必定無限長。（3）實際應(yīng)用中，所能觀察到的信號長度總是有限的。（4）因此，用DFT對連續(xù)時間信號進(jìn)行頻譜分析時，采樣后產(chǎn)生的頻譜混疊失真是不可避的。（5）得到的譜分析結(jié)果總是近似的，近似的程度取決于信號帶寬、采樣率以及信號的截取長度。

用DFT進(jìn)行連續(xù)信號頻譜分析時注意點（1）若信號持續(xù)時間有72例5-6：連續(xù)信號的頻譜分析

已知模擬信號對x(t)進(jìn)行采樣，得，分別以T=0.5s，0.1s進(jìn)行采樣，樣本數(shù)分別取N=20，100，以使被分析信號等長，用DFT計算的并與比較。例5-6：連續(xù)信號的頻譜分析已知模擬信號73例5-6的分析原信號x(t)的傅立葉變換為其幅度譜為，幅頻特性如圖5-17所示。由圖可知，|X(jΩ)|的-3dB點為Ω=1，而Ω=2π時已衰減到-16dB。因此，如取采樣頻率Ωs=4π，造成的頻率混疊效應(yīng)會比較小，此時相當(dāng)于fs=2,T=0.5s。在T=0.1s時，混迭將更小。例5-6的分析原信號x(t)的傅立葉變換為74例5-6解:x(t)采樣后所得信號的N點DFT按DFT定義求解題給連續(xù)信號x(t)經(jīng)采樣后x(n)的N點DFT實際情況下，信號表達(dá)式事先難以確知，因此解析式通常是無法得到的。

例5-6解:x(t)采樣后所得信號的N點DFT按DFT定義75兩種采樣率所得結(jié)果的比較兩種采樣率所得結(jié)果的比較76例5-6的結(jié)果分析實際應(yīng)用情況下，通常是截取連續(xù)信號的一段樣本進(jìn)行DFT分析。顯然，由于xa(t)本身非帶限，且又僅截取了[0,T0]的一段，因此，采樣后一定有頻率混疊（aliasing）存在，再加上截取后的樣本值x(n)并不等于前面式中所示的xN(n)，所以，x(n)在數(shù)字頻率[0,π]內(nèi)的DFT將不可能與相應(yīng)模擬頻率內(nèi)的Xa(jΩ)完全相同。造成誤差的兩個來源中，時域中x(n)與xN(n)的誤差在實際應(yīng)用中無法消除，但由于采樣率不夠而引入的頻率混疊可以通過減小T而得到減輕。比較例5-6在兩種采樣率下的結(jié)果，可以明顯看到提高采樣率的作用，但其代價為計算量的增加。例5-6的結(jié)果分析實際應(yīng)用情況下，通常是截取連續(xù)信號的一段樣77其他說明實際工作中，被分析數(shù)據(jù)是對信號進(jìn)行截取有限長度得到的，這相當(dāng)于被分析的是原信號乘上了一個窗函數(shù)。這種時域上的截斷，在頻域上就是所得頻譜是原信號的頻譜和窗函數(shù)頻譜的卷積，若信號頻率不位于DFT的譜線位置，將會產(chǎn)生頻譜的泄漏。在對實際信號進(jìn)行DFT分析時，泄漏效應(yīng)是難以避免的，只能通過選擇好的窗函數(shù)來抑制其影響。用DFT分析信號時，要求被分析的樣本序列在所截取的范圍內(nèi)可視為是穩(wěn)定的。利用DFT進(jìn)行頻譜分析時，DFT的計算結(jié)果是樣本序列以窗口時間為周期，進(jìn)行延拓后所得周期波形的傅里葉級數(shù)的系數(shù)，因此與信號的真實頻譜仍有差別。

其他說明實際工作中，被分析數(shù)據(jù)是對信號進(jìn)行截取有限長度得到的78fft函數(shù)與ifft函數(shù)MATLAB提供fft()函數(shù)來計算DFT，fft()函數(shù)是用機器語言而不是以MATLAB指令寫成的，因此它的運行速度很快。fft()函數(shù)的命令形式：

Y=fft(x)

-計算信號x的快速離散傅里葉變換Y。當(dāng)x為矩陣（多通道信號）時，計算x中每一列信號的離散傅立葉變換。

Y=fft(x,n) -計算信號x的n點FFT，當(dāng)x的長度大于n時，截斷x，否則補零。MATLAB還提供ifft()函數(shù)用于FFT的反變換計算，命令格式類似fft()函數(shù)。fft函數(shù)與ifft函數(shù)MATLAB提供fft()函數(shù)來計79例5-7

當(dāng)信號受到噪聲污染后，很難從時域中看出它包含的頻率成分，這時通?？梢杂肍FT來分析信號頻率成分。例5-7中給出了一個仿真產(chǎn)生的受到均勻分布隨機噪聲干擾的信號，例中用fft()函數(shù)進(jìn)行頻譜分析來找出其主要頻率成分，并通過filter()函數(shù)濾去噪聲的干擾。所得結(jié)果分別示于圖5-20與圖5-21中。

例5-7當(dāng)信號受到噪聲污染后，很難從時域中看出它包含的頻率805.6 練習(xí)課后完成p91~p92練習(xí)3-5。5.6 練習(xí)課后完成p91~p92練習(xí)3-5。81第5章使用MATLAB實現(xiàn)數(shù)字信號處理本章主要內(nèi)容如下：5．1 數(shù)字信號處理基本內(nèi)容及相應(yīng)的MATLAB工具5．2 信號通過系統(tǒng)的時域分析5．3 信號通過系統(tǒng)的頻域和Z域分析5．4 濾波器設(shè)計5．5 頻譜分析第5章使用MATLAB實現(xiàn)數(shù)字信號處理本章主要內(nèi)容如下：825.1數(shù)字信號處理基本內(nèi)容

及相應(yīng)的MATLAB工具數(shù)字信號處理的基本內(nèi)容通常分為兩部分：離散時間信號與系統(tǒng)分析主要涉及離散時間信號與系統(tǒng)的時域、頻域表示，以及信號通過系統(tǒng)的時域、頻域分析及其變換域分析。

MATLAB函數(shù)庫中提供了filter,conv,convmtx,fft,ifft,freqz,impz,zplane等與之相應(yīng)的函數(shù)。數(shù)字濾波器設(shè)計和譜分析數(shù)字濾波器設(shè)計包括了無限沖激響應(yīng)(IIR)和有限沖激響應(yīng)（FIR）濾波器設(shè)計，譜分析又可進(jìn)一步分為線性譜分析和非線性譜分析。MATLAB為此提供了多種成熟算法的相應(yīng)函數(shù)以及極為豐富的設(shè)計工具。5.1數(shù)字信號處理基本內(nèi)容

及相應(yīng)的MATLAB工具數(shù)字信835.2 時域分析卷積，濾波，單位沖激響應(yīng)5.2 時域分析卷積，濾波，單位沖激響應(yīng)845.2.1卷積MATLAB提供

conv函數(shù)實現(xiàn)標(biāo)準(zhǔn)的一維信號卷積:

例如，若系統(tǒng)h(n)為 >>h=[111]輸入序列x(n)為 >>x=[111]則x(n)經(jīng)過系統(tǒng)h(n)后的MATLAB實現(xiàn)為： >>conv(h,x)或

conv([111],[111])執(zhí)行后即得到y(tǒng)(n)為 ans= 12321注意：使用conv函數(shù)時，h(n)和x(n)都必須是有限長的,，否則不能使用conv函數(shù)。5.2.1卷積MATLAB提供conv函數(shù)實現(xiàn)標(biāo)準(zhǔn)的一維85例5-1時域離散序列的卷積計算

與圖形顯示

例5-1(教材p63)

：已知離散信號x(n)和h(n)，求y(n)=x(n)*h(n)，并用圖形表示。例5-1時域離散序列的卷積計算

與圖形顯示例5-1(教86例5-1的MATLAB程序Nh=20;Nx=10;m=5; %設(shè)定Nx，Nh和位移值mn=0:Nh-1;h1=(0.9).^n; %產(chǎn)生h1(n)h2=h1;nx=0:Nx-1;x1=ones(1,Nx); %產(chǎn)生x1(n)x2=zeros(1,Nx+m);fork=m+1:m+Nx %產(chǎn)生x2(n)=x1(n-m)x2(k)=x1(k-m);end %產(chǎn)生x2(n)y1=conv(x1,h1); %計算y1(n)=x1(n)*h1(n)y2=conv(x2,h2); %計算y2(n)=x2(n)*h2(n)subplot(3,2,1)stem(nx,x1,'.')axis([03001.2]),title(‘x1(n)’) %繪圖……（以下省略）例5-1的MATLAB程序Nh=20;Nx=10;m=5; 875.2.2濾波數(shù)字濾波器的系統(tǒng)函數(shù)H(z)用如下式表示:在MATLAB中，用向量b，a來表示濾波器的系數(shù)b(i)和

a(i)。

5.2.2濾波數(shù)字濾波器的系統(tǒng)函數(shù)H(z)用如下式表示:88濾波器分類當(dāng)n

=

0，m≠0時，稱為AR濾波器，即自回歸（AutoRecurrence）濾波器，具無限沖激響應(yīng)（IIR），也即其單位采樣響應(yīng)h(n)具無限長度；若m

=

0，a(1)≠0，稱為MA濾波器，即滑動平均（MovingAverage）濾波器，其單位采樣響應(yīng)h(n)是有限長度，故稱有限沖激響應(yīng)（FIR）濾波器；如果n、m都大于零，稱為ARMA濾波器，而其沖激響應(yīng)也為IIR。濾波器分類當(dāng)n

=

0，m≠0時，稱為AR濾波器，即自回歸（89filter函數(shù)

MATLAB提供了filter函數(shù)來對離散信號進(jìn)行濾波，表達(dá)信號通過系統(tǒng)后的結(jié)果。與conv不同的是，filter函數(shù)可適用于無限沖激響應(yīng)系統(tǒng)的情況，但信號仍須是有限長的。例如，一個單極點的低通濾波器系數(shù)如下：

>>b=1; %分子系數(shù)向量b(i) >>a=[1-0.9]; %分母系數(shù)向量a(i)

如果用filter函數(shù)實現(xiàn)對信號x濾波，只要調(diào)用： >>y=filter(b,a,x);

就可給出輸入x經(jīng)過濾波以后的輸出y。filter函數(shù)MATLAB提供了filter函數(shù)來對離905.2.3單位沖激響應(yīng)數(shù)字濾波器的單位沖激響應(yīng)定義為輸入為單位樣本序列時數(shù)字濾波器的響應(yīng)，即：h(n)=T[δ(n)] 其中：5.2.3單位沖激響應(yīng)數(shù)字濾波器的單位沖激響應(yīng)定義為輸入91單位沖激響應(yīng)的MATLAB實現(xiàn)MATLAB近似實現(xiàn)單位采樣信號的方法為：

imp=[1;zeros(p,1)]; %zeros(p,1)產(chǎn)生p個零元素組成的列向量，p是正整數(shù)。使用imp后，濾波器的沖激響應(yīng)可近似得到為：

h=filter(b,a,imp);impz函數(shù)可以直接求出數(shù)字濾波器的單位沖激響應(yīng)，即：impz(b,a) 該命令將同時繪出濾波器的單位沖激響應(yīng)，教材p66圖5-2。單位沖激響應(yīng)的MATLAB實現(xiàn)MATLAB近似實現(xiàn)單位采樣信925．3頻域和Z域分析

頻率響應(yīng)，零極點分析

5．3頻域和Z域分析頻率響應(yīng)，零極點分析935.3.1頻率響應(yīng)MATLAB數(shù)字信號處理工具箱有很多函數(shù)提供對模擬和數(shù)字濾波器的頻率響應(yīng)分析。其中，freqz函數(shù)和freqs函數(shù)分別返回數(shù)字和模擬濾波器的頻率響應(yīng)。工具箱中通常使用的單位頻率是Nyquist頻率，即采樣頻率的1/2。注意：就數(shù)字濾波器函數(shù)來說，其頻域指標(biāo)中的所有頻率都以Nyquist頻率進(jìn)行歸一化。因此Nyquist頻率也稱歸一化頻率。5.3.1頻率響應(yīng)MATLAB數(shù)字信號處理工具箱有很多函94關(guān)于Nyquist頻率的說明例如：系統(tǒng)采樣頻率為1000Hz，則若數(shù)字濾波器的截止頻率等于300Hz，經(jīng)Nyquist頻率歸一化后，其歸一化頻率就是300/500=0.6。若將歸一化頻率轉(zhuǎn)換成數(shù)字信號處理教科書中所使用的數(shù)字頻率(rad)，需乘以π；反之，若乘以采樣頻率fs的一半，則將歸一化頻率轉(zhuǎn)換回了模擬域頻率(Hz)。歸一化頻率f應(yīng)滿足0<f<1。對于頻率響應(yīng)而言，歸一化頻率和模擬頻率都可使用。關(guān)于Nyquist頻率的說明例如：系統(tǒng)采樣頻率為1000Hz95freqz命令

freqz

函數(shù)提供了基于FFT算法所得到的數(shù)字濾波器H(z)的頻率響應(yīng)最常用的形式：[h,w]=freqz(b,a,p)

其意義是，對于數(shù)字濾波器： freqz函數(shù)接受濾波器系數(shù)向量b和a，以及一個用于指定所需計算的頻率響應(yīng)樣點數(shù)的整數(shù)p，返回一個與由p個頻率樣本點構(gòu)成的實頻率向量相對應(yīng)的復(fù)頻率響應(yīng)向量h。

freqz命令

freqz函數(shù)提供了基于FFT算法所得到的96freqz的命令形式h=freqz(b,a,w)采用上面的形式時，需先對頻率樣本點向量w作出定義。通常的做法是使用linspace函數(shù)。[h,w]=freqz(b,a)w和p沒有定義，默認(rèn)w由（0~π）上均分的512點構(gòu)成，頻率單位為rad/sample。[h,w]=freqz(b,a,p，’whole’)使用帶參數(shù)’whole’選項的命令形式[h,w]=freqz(b,a,p,fs)用參數(shù)選項指定采樣頻率fsfreqz的命令形式h=freqz(b,a,w)97linspace的命令形式linspace函數(shù)的命令形式有兩種：linspace（x1,x2） %產(chǎn)生[x1,x2]范圍內(nèi)經(jīng)線性分割得到的100點的向量linspace（x1,x2，N） %產(chǎn)生[x1,x2]范圍內(nèi)經(jīng)線性分割得到的N點的向量

例如：w=linspace（0，pi）%

w由（0，pi）上100個等分點組成 w=linspace（0，1000）

%w由（0，1000）上100個等分點組成linspace的命令形式linspace函數(shù)的命令形式有兩98freqz函數(shù)無返回輸出參數(shù)格式的調(diào)用freqz函數(shù)還能以無返回輸出參數(shù)的格式調(diào)用： freqz（b，a，256）或freqz（b，a，256，2000）無返回輸出參數(shù)調(diào)用freqz函數(shù)的好處是，MATLAB能夠自動繪出頻率在（0~π）范圍內(nèi)的幅頻特性和相頻特性圖。注意，無返回輸出參數(shù)調(diào)用freqz函數(shù)繪出的相頻特性不能正確給定在處的值，這是因為在MATLAB中，屬于下半個單位圓。

freqz函數(shù)無返回輸出參數(shù)格式的調(diào)用freqz函數(shù)還能以無99頻率響應(yīng)的實例例：先構(gòu)成一個截止頻率為400Hz的9階巴特沃思（Butterworth）低通數(shù)字濾波器，求出其系數(shù)b，a，再求出其256點頻率響應(yīng)。指定的采樣頻率fs=2000Hz。實現(xiàn)1：先調(diào)用butter函數(shù)，再調(diào)用freqz函數(shù)；實現(xiàn)2：無返回輸出參數(shù)，調(diào)用freqz函數(shù)；實現(xiàn)3：為得到完整的（0~π）范圍內(nèi)的幅相特性圖，可用帶返回輸出參數(shù)的格式調(diào)用freqz，在得到各頻率點上的頻率響應(yīng)的數(shù)值后，使用MATLAB所提供的abs函數(shù)和angle函數(shù)分別提取幅頻特性與相頻特性。注意：為了得到連續(xù)的相頻曲線，需要使用unwrap函數(shù)，該函數(shù)通過在發(fā)生跳變后的各處自動加上，從而除去了相位的跳變，這稱為相位的“解卷繞”。

頻率響應(yīng)的實例例：先構(gòu)成一個截止頻率為400Hz的9階巴特沃100實現(xiàn)1

[b,a]=butter(9,400/1000); %括號中第一個輸入是階數(shù)，第二個是歸一化截止頻率

[h,f]=freqz(b,a,256,2000);或

h=freqz（b，a，256，2000）；實現(xiàn)1 [b,a]=butter(9,400/1000)101實現(xiàn)2實現(xiàn)1中第二條命令的形式可改寫為freqz（b，a，256）

或freqz（b，a，256，2000）可自動繪出頻率在（0~π）范圍內(nèi)的幅頻特性和相頻特性圖。注意無返回輸出參數(shù)調(diào)用freqz函數(shù)繪出的相頻特性不能正確給定在ω=π處的值，因為ω=π屬于下半個單位圓。實現(xiàn)2實現(xiàn)1中第二條命令的形式可改寫為102實現(xiàn)3[b,a]=butter(9,400/1000);w=linspace(0,1000);h=freqz(b,a,w,2000);mag=abs(h);pha=angle(h);%提取濾波器頻率響應(yīng)的幅度mag和相位phasemilogy(w,mag);title('Magnitude');figure;plot(w,pha);title('Phase');出現(xiàn)了幅度為2π的跳變——卷繞實現(xiàn)3[b,a]=butter(9,400/1000);103解卷繞：unwrap函數(shù)為了得到連續(xù)的相頻曲線，需要使用unwrap函數(shù)，該函數(shù)通過在發(fā)生2π跳變后的各處自動加上2π，從而除去了相位的跳變，這稱為相位的“解卷繞”。解卷繞解卷繞：unwrap函數(shù)為了得到連續(xù)的相頻曲線，需要使用u104freqs函數(shù)freqs函數(shù)用于求取由b和a系數(shù)向量定義的模擬濾波器 H（s）的頻率響應(yīng)。使用方法類似freqz函數(shù)。與第二章（p32例2-1）采用數(shù)組相除方法求取頻率響應(yīng)相比，使用freqs函數(shù)要方便很多。freqs函數(shù)freqs函數(shù)用于求取由b和a系數(shù)向量定義的模1055.3.2零極點分析zplane函數(shù)用于畫出線性系統(tǒng)在Z平面上的零極點。有兩種使用方法：

1、在已知零極點時，例如某濾波器的零點為-1/2，一對共軛極點為和時，只要輸入命令zer=-0.5;pol=0.9*exp(j*2*pi*[-0.30.3]');zplane(zer,pol)

即可畫出零極點。（見p70圖5-6）5.3.2零極點分析zplane函數(shù)用于畫出線性系統(tǒng)106用zplane函數(shù)繪制零極點圖另一種情況：已知系統(tǒng)的系統(tǒng)函數(shù)系數(shù)向量b和a，則可通過調(diào)用zplane(b,a)繪出零極點。這種情形下，zplane函數(shù)先求得系統(tǒng)函數(shù)的零點和極點，然后繪出零極點圖。用zplane函數(shù)繪制零極點圖另一種情況：已知系統(tǒng)的系統(tǒng)函數(shù)1075.4 濾波器設(shè)計濾波器技術(shù)指標(biāo)，IIR濾波器設(shè)計，F(xiàn)IR濾波器設(shè)計5.4 濾波器設(shè)計濾波器技術(shù)指標(biāo)，1085.4.1濾波器的技術(shù)指標(biāo)設(shè)計濾波器之前，要定義一組濾波器的技術(shù)指標(biāo)。濾波器的技術(shù)指標(biāo)通常采用幅度指標(biāo)給出，以保證物理上的可實現(xiàn)性。幅度指標(biāo)以容差方式給出，如p71圖5-7所示。具體指標(biāo)：通帶截止頻率Wp；阻帶截止頻率Ws；通帶波紋δp；阻帶衰減δs。其中通帶波紋和阻帶衰減也可以用分貝數(shù)給出5.4.1濾波器的技術(shù)指標(biāo)設(shè)計濾波器之前，要定義一組濾波109濾波器以容差方式給出的幅度指標(biāo)

|H(jw)|10通帶過渡帶阻帶pwswsdpd-1w濾波器以容差方式給出的幅度指標(biāo)|H(jw)|10通帶過渡1105.4.2IIR濾波器設(shè)計MATLAB工具箱提供了幾種模擬濾波器原型的產(chǎn)生函數(shù),如巴特沃思(Butterworth),切比雪夫(Chebyshev)濾波器等。還提供了從模擬低通濾波器原型轉(zhuǎn)換為高通、帶通和帶阻的轉(zhuǎn)換函數(shù)，模擬濾波器轉(zhuǎn)換為數(shù)字濾波器的雙線性變換法和沖激響應(yīng)不變法，模擬IIR濾波器的階數(shù)選擇函數(shù)以及數(shù)字濾波器直接設(shè)計函數(shù)等等，使用起來非常方便。本節(jié)以巴特沃思濾波器為例

介紹MATLAB工具箱提供的濾波器設(shè)計函數(shù)使用方法。5.4.2IIR濾波器設(shè)計MATLAB工具箱提供了幾1111．巴特沃思濾波器設(shè)計函數(shù)

[b,a]

=

butter(n,Wn,options)缺省情況下，返回一個用有理分式表示的低通數(shù)字濾波器系統(tǒng)函數(shù)的分子分母系數(shù)向量b和a。

設(shè)計的技術(shù)指標(biāo)只需要指定一個歸一化截止頻率Wn（這里Wn的單位為：Nyquist

頻率

=

1Hz）。參數(shù)選項options：選項空缺返回低通濾波器；‘high’則表示返回的是高通數(shù)字濾波器；至于帶通和帶阻濾波器設(shè)計，則輸入的Wn必須是二個元素的向量；選項options空缺返回帶通濾波器，帶阻濾波器需加參數(shù)options為'stop'。

1．巴特沃思濾波器設(shè)計函數(shù)[b,a]

=

butter(n112Butter函數(shù)舉例[B,A]=butter(N,Wn)返回系統(tǒng)函數(shù)的長度為N+1的分子系數(shù)B和分母系數(shù)A。式中Wn為歸一化截止頻率，0<Wn<1.0。[B,A]=butter(N,Wn,'high')設(shè)計一個高通數(shù)字濾波器。如果Wn=[W1W2]，且options空缺，則返回一個2N階的帶通濾波器，其通帶為：W1<W<W2。[B,A]=butter(N,Wn,'stop')，且Wn=[W1W2]設(shè)計2N階的帶阻濾波器。butter(N,Wn,'s')，butter(N,Wn,'high','s')或

butter(N,Wn,'stop','s')調(diào)用butter()函數(shù)附加一個options為‘s’時，表示設(shè)計的是一個模擬濾波器，而指定的截止頻率Wn單位此時必須為rad/s。Butter函數(shù)舉例[B,A]=butter(N,Wn1132.巴特沃思濾波器階次選擇函數(shù)buttord函數(shù)用于給出滿足給定頻域指標(biāo)要求的最小階次濾波器的設(shè)計。其調(diào)用格式為：

[n,Wn]

=

buttord(Wp,Ws,Rp,Rs) %用于數(shù)字濾波器[n,Wn]

=

buttord(Wp,Ws,Rp,Rs,'s')%用于模擬濾波器該函數(shù)返回符合技術(shù)指標(biāo)的最小階數(shù)N以及Butterworth濾波器3dB截止頻率Wn。設(shè)計指標(biāo)同前：通帶起伏不超過Rp，阻帶衰減不小于Rs，Wp和Ws分別為歸一化通帶和阻帶截止頻率。

2.巴特沃思濾波器階次選擇函數(shù)buttord函數(shù)用于給出滿足114Buttord函數(shù)應(yīng)用：例5-1

例：設(shè)計一個數(shù)字帶通濾波器，通帶為1000到2000Hz，止帶的起始位置離開通帶兩邊500Hz，也即分別為500Hz和2500Hz，設(shè)采樣率為10KHz，通帶起伏為1dB,止帶最小衰減為60dB。解MATLAB實現(xiàn)如下：

[n,Wn]=buttord([10002000]/5000,[5002500]/5000,1,60)得到： n= 12 Wn= 0.19510.4080因此 [b,a]=butter(n,Wn);給出了滿足要求的最小階次巴特沃思帶通濾波器。Buttord函數(shù)應(yīng)用：例5-1例：設(shè)計一個數(shù)字帶通濾波器1153．經(jīng)典的IIR濾波器設(shè)計方法一種常用的IIR濾波器設(shè)計方法是利用已經(jīng)很成熟的模擬濾波器設(shè)計結(jié)果。設(shè)計步驟如下：（1）按照一定規(guī)則將給出的數(shù)字濾波器的技術(shù)指標(biāo)轉(zhuǎn)換為模擬低通濾波器的技術(shù)指標(biāo)。 （2）根據(jù)轉(zhuǎn)換后的技術(shù)指標(biāo)設(shè)計模擬低通濾波器的H(s)。 （3）再按一定規(guī)則將H(s)轉(zhuǎn)換成H(z)。步驟(3)中最常用的方法:沖激響應(yīng)不變法（也克稱為脈沖響應(yīng)不變法）和雙線性變換法。

3．經(jīng)典的IIR濾波器設(shè)計方法一種常用的IIR濾波器設(shè)計方法116沖激響應(yīng)不變法的原理

沖激響應(yīng)不變法是通過對模擬濾波器的單位沖激響應(yīng)ha(nTs),進(jìn)行采樣，將其作為數(shù)字濾波器的單位沖激響應(yīng)h(n)。即：據(jù)此設(shè)計得到的數(shù)字濾波器的頻率響應(yīng)與模擬濾波器的頻率響應(yīng)之間的關(guān)系為：因此存在著頻率混迭。從而限制了其適用范圍。注意：S域頻率軸與Z域單位圓之間不存在一一對應(yīng)關(guān)系。沖激響應(yīng)不變法的原理沖激響應(yīng)不變法是通過對模擬濾波器的單位117沖激響應(yīng)不變設(shè)計的impinvar()函數(shù)impinvar()函數(shù)的調(diào)用格式是：

[Bz,Az]=impinvar(Bs,As,fs)這一函數(shù)把具有[Bs,As]模擬濾波器的轉(zhuǎn)換成了采樣頻率為fs的數(shù)字濾波器的[Bz,Az]。如果沒有給定采樣頻率fs，函數(shù)默認(rèn)為1Hz。

沖激響應(yīng)不變設(shè)計的impinvar()函數(shù)impinvar(118雙線性變換法原理雙線性變換法所定義的s平面與z平面的映射關(guān)系為：給定模擬濾波器的系統(tǒng)函數(shù)，數(shù)字濾波器為：因此S域內(nèi)的有理分式函數(shù)在Z域內(nèi)仍然保持了有理分式函數(shù)形式，從而保證了數(shù)字濾波器的物理可實現(xiàn)性。