《數(shù)字信號(hào)處理實(shí)驗(yàn)》指導(dǎo)書
數(shù)字信號(hào)處理實(shí)驗(yàn)指導(dǎo)書 指導(dǎo)教師:朱士永數(shù)字信號(hào)處理實(shí)驗(yàn)指導(dǎo)書實(shí)驗(yàn)一、離散信號(hào)的產(chǎn)生運(yùn)算(一)實(shí)驗(yàn)?zāi)康暮鸵螅寒a(chǎn)生離散信號(hào),進(jìn)行信號(hào)運(yùn)算(二)主要儀器設(shè)備: 計(jì)算機(jī);Matlab軟件;word文字處理軟件(三)實(shí)驗(yàn)原理與方法(1)產(chǎn)生實(shí)指數(shù)序列:n=0:50;x=0.9.n;stem(x);(2)正余弦序列:n=0:50;x=sin(0.075*pi*n);% x=cos(0.075*pi*n);stem(x);(3)產(chǎn)生隨機(jī)序列 x=10.*rand(1,30);%產(chǎn)生長(zhǎng)度為30,大小在(0,10)間的一維隨機(jī)向量stem(x);%繪圖(4)序列反褶:y=fliplr(x);%產(chǎn)生序列x的反折序列(5)移位移位函數(shù)文件:sigshift.mfunctiony,ny=sigshift(x,nx,n0) ;ny=nx+n0;%右移n0位;y=x;matlab實(shí)現(xiàn)調(diào)用文件test1.m文件:x=3,11,7,0,-1,4,2;nx=-3:3;y,ny=sigshift(x,nx,2);%右移2位subplot(3,1,1);stem(nx,x);ylabel('x');subplot(3,1,2);stem(ny,y);ylabel('y(n)');(四)實(shí)驗(yàn)步驟與內(nèi)容1.編程產(chǎn)生長(zhǎng)度為10,大小在(0,5)間的一維隨機(jī)向量,在下面表格記錄序列值,根據(jù)表格值繪出序列的圖形。n12345678910x(n)2.編程產(chǎn)生序列:,在下面表格記錄序列值,根據(jù)表格值繪出序列的圖形。n12345678910y(n)3.重新編程繪圖實(shí)現(xiàn)上兩個(gè)序列的相加序列A(n)=x(n)+y(n)、相乘運(yùn)算序列B(n)=x(n).y(n), 記錄在下表格中,根據(jù)表格繪出序列A(n), B(n)的圖形。n12345678910A(n)B(n)4.實(shí)現(xiàn)2中序列的反褶運(yùn)算序列Z(n),并記錄Z(n)在下列表格中,根據(jù)表格值繪出序列Z(n)的圖形。n12345678910Z(n)(五)實(shí)驗(yàn)報(bào)告要求1.寫出實(shí)驗(yàn)實(shí)驗(yàn)報(bào)告,包括:目的和要求、儀器 2.寫出實(shí)驗(yàn)步驟及內(nèi)容,給出程序,填寫實(shí)驗(yàn)數(shù)據(jù)表格,繪出圖形實(shí)驗(yàn)二、卷積與相關(guān)運(yùn)算(一)實(shí)驗(yàn)?zāi)康暮鸵笥胢atlab編寫DSP程序,觀察輸入序列的長(zhǎng)度與輸出序列長(zhǎng)度的關(guān)系(二)主要儀器設(shè)備計(jì)算機(jī);Matlab軟件;word文字處理軟件(三)實(shí)驗(yàn)原理與方法(1)離散序列的卷積:x=3,11,7,0,-1,4,2;nx=-3:3;h=2,3,0,-5,2,1;nh=-1:4;nyb=nx(1)+nh(1);nye=nx(length(x)+nh(length(h);ny=nyb:nyey=conv(x,h)subplot(3,1,1);stem(nx,x);ylabel('x');subplot(3,1,2);stem(nh,h);ylabel('h');subplot(3,1,3);stem(ny,y);ylabel('y(n)=x(n)*h(n)');(2)離散序列的相關(guān):x=3,11,7,0,-1,4,2;nx=-3:3;h=2,3,0,-5,2,1;nh=-1:4;nyb=nx(1)+nh(1);nye=nx(length(x)+nh(length(h);ny=1:2*max(length(x),length(h)-1;y=xcorr(x,h);%相關(guān)subplot(3,1,1);stem(nx,x);ylabel('x(nx)');subplot(3,1,2);stem(nh,h);ylabel('h(nh)');subplot(3,1,3);stem(ny,y);ylabel('y(ny)');(四)實(shí)驗(yàn)步驟與內(nèi)容1.編程實(shí)現(xiàn)序列:和的卷積,記錄在下表格:n012345678910111213y(n)根據(jù)表格值繪出卷積序列y(n)的圖形。2.編程實(shí)現(xiàn)1中兩個(gè)序列的相關(guān)運(yùn)算序列,記錄在下表格:n12345678910z(n)n111213141516171819z(n)繪出卷積序列的圖形。3.編程實(shí)現(xiàn)序列:和的卷積,并繪圖nz(n)nz(n)4.編程實(shí)現(xiàn)3中兩個(gè)序列的相關(guān)運(yùn)算序列并繪圖nz(n)nz(n)(五)實(shí)驗(yàn)報(bào)告要求1.寫出實(shí)驗(yàn)實(shí)驗(yàn)報(bào)告,包括:目的和要求、儀器 2.寫出實(shí)驗(yàn)步驟及內(nèi)容,給出程序,填寫實(shí)驗(yàn)數(shù)據(jù)表格,繪出圖形實(shí)驗(yàn)三、離散系統(tǒng)分析(一)實(shí)驗(yàn)?zāi)康暮鸵笄蟪鲭x散系統(tǒng)的零極點(diǎn),分析離散系統(tǒng)的幅頻和相頻特性(二)主要儀器設(shè)備計(jì)算機(jī);Matlab軟件;word文字處理軟件(三)實(shí)驗(yàn)原理與方法1 繪制離散系統(tǒng)的頻率響應(yīng):模擬:或 按照矢量w指定的頻率值沿虛坐標(biāo)計(jì)算模擬濾波器的復(fù)頻率響應(yīng)h,當(dāng)w省略時(shí)自動(dòng)取w=200,不代變量輸出的直接繪出曲線。數(shù)字:或按照n個(gè)頻率點(diǎn)計(jì)算頻率響應(yīng)h,n個(gè)頻率點(diǎn)均勻分布在上半單位圓(0),并記錄在w中,相應(yīng)的頻率響應(yīng)記錄在h中。若不帶返回參數(shù),則繪制出幅頻特性和相頻特性曲線。2繪制離散系統(tǒng)零極點(diǎn)圖調(diào)用格式:zplane(z,p)zplane(b,a)其中b、a分別為系統(tǒng)函數(shù)分子分母的系數(shù)向量,z為系統(tǒng)零點(diǎn)向量,p為系統(tǒng)極點(diǎn)向量。例如:數(shù)字濾波器的零極點(diǎn)圖、幅頻特性和相頻特性曲線程序:a=1 0.5 0.1;b=0.2 0.3 0;zplane(b,a); %繪制零極點(diǎn)圖freqz(b,a,128);%繪制128個(gè)頻率點(diǎn)的幅頻特性和相頻特性曲線3數(shù)字系統(tǒng)的單位脈沖響應(yīng)1): t記錄取樣點(diǎn)矢量,取樣點(diǎn)n由系統(tǒng)自動(dòng)選取。2):用戶指定取樣點(diǎn)n,當(dāng)n為標(biāo)量時(shí),t=0:n-1不帶輸出變量的impz將在當(dāng)前圖形窗口利用stem(t,h)繪出脈沖響應(yīng)。4系統(tǒng)函數(shù)轉(zhuǎn)換(1)系統(tǒng):函數(shù)tf2zp可以將系統(tǒng)函數(shù)表示成零極點(diǎn)形式:調(diào)用格式:z,p,k=tf2zp(b,a);%求出系統(tǒng)的零極點(diǎn)其中b、a分別為系統(tǒng)函數(shù)分子分母的系數(shù)向量,z為系統(tǒng)零點(diǎn)向量,p為系統(tǒng)極點(diǎn)向量,k為增益系數(shù)。(2)函數(shù)zp2tf可以將零極點(diǎn)形式轉(zhuǎn)化為分子分母多項(xiàng)式形式調(diào)用格式:b,a=zp2tf(z,p,k)(四)實(shí)驗(yàn)步驟與內(nèi)容1.編程求出系統(tǒng)的零極點(diǎn)記錄在下表,并在Z平面繪制出零極點(diǎn)圖類型零點(diǎn)極點(diǎn)2.編程繪出上述1中系統(tǒng)的幅頻特性和相頻特性曲線3.編程繪出系統(tǒng)的零極點(diǎn)記錄在下表,并在Z平面繪制出零極點(diǎn)圖類型零點(diǎn)極點(diǎn)4.編程繪出上述3中系統(tǒng)的幅頻特性和相頻特性曲線5.繪出系統(tǒng)1和系統(tǒng)2的單位脈沖響應(yīng)(五)實(shí)驗(yàn)報(bào)告要求1.寫出實(shí)驗(yàn)實(shí)驗(yàn)報(bào)告,包括:目的和要求、儀器 2.寫出實(shí)驗(yàn)步驟及內(nèi)容,給出程序,填寫實(shí)驗(yàn)數(shù)據(jù)表格,繪出圖形實(shí)驗(yàn)四、信號(hào)譜分析(一)實(shí)驗(yàn)?zāi)康暮鸵笳莆誇FT算法,利用MATLAB語(yǔ)言對(duì)信號(hào)進(jìn)行頻譜分析(二)主要儀器設(shè)備計(jì)算機(jī);Matlab軟件;word文字處理軟件(三)實(shí)驗(yàn)原理與方法一維快速傅立葉變換函數(shù)fft,調(diào)用格式:(1)X=fft(x) %返回x的相同長(zhǎng)度的fft(2)X=fft(x,n) %計(jì)算n點(diǎn)的fft調(diào)用函數(shù)完成:,其中,則信號(hào)的幅度譜為:或相位譜為:例如:x=1,2,4,1,0,3,5X=fft(x,128)p=sqrt(X.*conj(X); %或ph1=abs(X);ph=angle(X);subplot(2,1,1);plot(p);subplot(2,1,2);plot(ph);(四)實(shí)驗(yàn)步驟與內(nèi)容 (都取128點(diǎn)的FFT)1.編程繪圖實(shí)現(xiàn)數(shù)字信號(hào):x=10,20,30,90,10,60,60,20,5的頻譜(幅頻和相頻)特性2.編程繪圖實(shí)現(xiàn)數(shù)字信號(hào):x=1,2,3,1,2,1的頻譜(幅頻和相頻)特性3.在采樣時(shí)間點(diǎn)t=0:0.001:0.6,對(duì)信號(hào)進(jìn)行離散處理,編程繪圖實(shí)現(xiàn)離散信號(hào)頻譜(幅頻和相頻)特性(五)實(shí)驗(yàn)報(bào)告要求1.寫出實(shí)驗(yàn)實(shí)驗(yàn)報(bào)告,包括:目的和要求、儀器 2.寫出實(shí)驗(yàn)步驟及內(nèi)容,給出程序,繪出圖形實(shí)驗(yàn)五、IIR數(shù)字濾波器的設(shè)計(jì)(一)實(shí)驗(yàn)?zāi)康暮鸵笳莆针p線性變換法及脈沖相應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的具體設(shè)計(jì)方法及其原理。用雙線性變換法及脈沖響應(yīng)不變法編程設(shè)計(jì)低通、高通和帶通IIR數(shù)字濾波器。(二)主要儀器設(shè)備:計(jì)算機(jī);Matlab軟件;word文字處理軟件(三)實(shí)驗(yàn)原理與方法1.低通模擬濾波器原型函數(shù)z,p,k=besselap(n) %貝塞爾 濾波器z,p,k=buttap(n) %Butterworth巴特沃斯濾波器z,p,k=cheb1ap(n,Rp) %雪比契夫1型z,p,k=cheb2ap(n,Rs) %雪比契夫2型z,p,k=ellipap(n,Rp,Rs) %橢圓濾波器2.模擬頻率變換函數(shù)(1)原型低通到低通模擬濾波器變換函數(shù)lp2lp:bt,at=lp2lp(b,a,Wo);其中W0為截止頻率(2)原型低通到帶通模擬濾波器變換函數(shù)lp2bp:bt,at=lp2bp(b,a,Wo,Bw);其中W0為中心頻率,Bw為帶寬(3)原型低通到高通模擬濾波器變換函數(shù)lp2hp:bt,at=lp2hp(b,a,Wo);其中W0為截止頻率(4)原型低通到帶阻模擬濾波器變換函數(shù)lp2bsbt,at=lp2bs(b,a,Wo,Bw) ;其中W0為中心頻率,Bw為帶寬3.雙線性變換法函數(shù)zd,pd,kd=bilinear(z,p,k,Fs) % Fs為取樣頻率zd,pd,kd=bilinear(z,p,k,Fs,Fp)numd,dend=bilinear(num,den,Fs) % num,den為分子分母多項(xiàng)式系數(shù)向量numd,dend=bilinear(num,den,Fs,Fp)4.沖擊響應(yīng)不變法函數(shù)bz,az=impinvar(b,a,Fs)bz,az=impinvar(b,a)例如:取采樣頻率f=1KHz,用雙線性變換法設(shè)計(jì)五階Butterworth低通數(shù)字濾波器,繪出模擬濾波器與數(shù)字濾波器的幅頻與相頻特性,MATLAB程序如下: z,p,k=buttap(5) ;% 設(shè)計(jì)五階Butterworth低通模擬濾波器原型zd,pd,kd=bilinear(z,p,k,1000);%雙線性變換得到低通數(shù)字濾波器b,a=zp2tf(zd,pd,kd);%濾波器類型轉(zhuǎn)換w=128;freqs(b,a,w)figure;freqz(b,a,w)(四)實(shí)驗(yàn)步驟與內(nèi)容1.取采樣頻率f=100Hz,用雙線性變換法設(shè)計(jì)五階Butterworth低通數(shù)字濾波器,繪出模擬濾波器與數(shù)字濾波器的幅頻與相頻特性2.設(shè)高通截止頻率為w0=10000Hz, 取采樣頻率f=20000,用雙線性變換法設(shè)計(jì)六階高通Butterworth數(shù)字濾波器,繪出數(shù)字濾波器的頻譜特性3.設(shè)帶通濾波器的濾波器中心頻率為W0=2KHz,帶寬為BW=100Hz, 取采樣頻率f=10kHZ,用脈沖相應(yīng)不變法設(shè)計(jì),設(shè)計(jì)五階帶通Butterworth數(shù)字濾波器,繪出數(shù)字濾波器的頻譜特性4.直接設(shè)計(jì)五階butterworth帶通濾波器,繪出頻譜圖。(高端與低端截止頻率分別為0.2和0.9)(五)實(shí)驗(yàn)報(bào)告要求1.寫出實(shí)驗(yàn)實(shí)驗(yàn)報(bào)告,包括:目的和要求、儀器 2.寫出實(shí)驗(yàn)步驟及內(nèi)容,給出程序,繪出圖形實(shí)驗(yàn)六、FIR數(shù)字濾波器的設(shè)計(jì)(一)實(shí)驗(yàn)?zāi)康暮鸵笳莆沼么昂瘮?shù)法,頻率采樣法及優(yōu)化設(shè)計(jì)法設(shè)計(jì)FIR濾波器的原理及方法,熟悉響應(yīng)的計(jì)算機(jī)編程,熟悉線性相位FIR濾波器的幅頻特性和相頻特性(二)主要儀器設(shè)備計(jì)算機(jī);Matlab軟件;word文字處理軟件(三)實(shí)驗(yàn)原理與方法1、常用的窗函數(shù)(1)矩形序列: w=boxcar(n);% w序列長(zhǎng)度為n(2)三角窗函數(shù): w=triangle(n)(3)Bartlett窗: w=bartlett(n)(4)漢明窗: w=hamming(n)(5)漢寧窗: w=hanning(n)(6)布萊克曼窗: w=Blackman(n)(7)Chebyshev窗: w=chebwin(n,r)(8)凱澤窗: w=Kaiser(n,beta)2.基于窗函數(shù)的FIR濾波器設(shè)計(jì)函數(shù)fir1b=fir1n,wn;%b為n階FIR濾波器的單位取樣響應(yīng)(長(zhǎng)度為n+1),wn為FIR低通濾波器的截止頻率b=fir1n,wn,ftype;b=fir1n,wn,window;%缺省window時(shí)按照hamming窗截取b=fir1n,wn,ftype,window;3.基于頻率采樣法和窗函數(shù)設(shè)計(jì)具有任意頻率響應(yīng)的FIR濾波器fir2b=fir2n,f,m;%n為階次,f為取樣頻率點(diǎn)向量且0f1遞增,%m為與f對(duì)應(yīng)的取樣相對(duì)幅度向量,m與f長(zhǎng)度相同b=fir2n,f,m,window;%指定窗口函數(shù)截取時(shí)域序列,缺省為hamming窗4.最小二乘法線性相位FIR濾波器設(shè)計(jì)函數(shù)firlsb=firls(n,f,m)b=firls(n,f,m,w)b=firls(n,f,m,ftype)b=firls(n,f,m,w,ftype)例如:分別采樣窗函數(shù)法和頻率取樣法設(shè)計(jì)11階FIR帶通數(shù)字濾波器,截止頻率為0.4和0.9,繪出濾波器單位取樣響應(yīng)和頻率特性。MATLAB設(shè)計(jì)如下:w=0.4,0.9;h1=fir1(11,w,'high');stem(h1);figure;freqz(h1,1);f=0,0.4,0.9,1;m=0,1,1,0;h2=fir2(11,f,m);figure;stem(h2);figure;freqz(h2,1);(四)實(shí)驗(yàn)步驟與內(nèi)容1.分別用矩形窗、漢寧窗、漢明窗、Blackman、Kaise(=8.5)窗設(shè)計(jì)10階(N=11)FIR數(shù)字濾波器,使得其頻譜特性滿足:繪出濾波器單位取樣響應(yīng)和頻率特性 2.采用頻率采樣法和Blackman窗函數(shù)設(shè)計(jì)一個(gè)FIR低通數(shù)字濾波器,取樣響應(yīng)長(zhǎng)度為32,截止頻率取。3采用頻率采樣法和Blackman窗函數(shù)設(shè)計(jì)一個(gè)10階FIR高通數(shù)字濾波器,截止頻率取。(五)實(shí)驗(yàn)報(bào)告要求1.寫出實(shí)驗(yàn)實(shí)驗(yàn)報(bào)告,包括:目的和要求、儀器 2.寫出實(shí)驗(yàn)步驟及內(nèi)容,給出程序,繪出圖形部分實(shí)驗(yàn)解答:實(shí)驗(yàn)五、IIR數(shù)字濾波器的設(shè)計(jì)(三)實(shí)驗(yàn)原理與方法1.低通模擬濾波器原型函數(shù)z,p,k=besselap(n) z,p,k=buttap(n)z,p,k=cheb1ap(n,Rp)z,p,k=cheb2ap(n,Rs)z,p,k=ellipap(n,Rp,Rs)2.模擬頻率變換函數(shù)(1)原型低通到低通模擬濾波器變換函數(shù)lp2lpbt,at=lp2lp(b,a,Wo)(2)原型低通到帶通模擬濾波器變換函數(shù)lp2bpbt,at=lp2bp(b,a,Wo,Bw)(3)原型低通到高通模擬濾波器變換函數(shù)lp2hpbt,at=lp2hp(b,a,Wo)(4)原型低通到帶阻模擬濾波器變換函數(shù)lp2bsbt,at=lp2bs(b,a,Wo,Bw)3.雙線性變換法函數(shù)zd,pd,kd=bilinear(z,p,k,Fs)zd,pd,kd=bilinear(z,p,k,Fs,Fp)numd,dend=bilinear(num,den,Fs)numd,dend=bilinear(num,den,Fs,Fp)4.沖擊響應(yīng)不變法函數(shù)bz,az=impinvar(b,a,Fs)bz,az=impinvar(b,a)(四)實(shí)驗(yàn)步驟與內(nèi)容例.取采樣頻率f=1KHz,用雙線性變換法設(shè)計(jì)五階Butterworth低通數(shù)字濾波器,繪出模擬濾波器與數(shù)字濾波器的幅頻與相頻特性 z,p,k=besselap(5) ;zd,pd,kd=bilinear(z,p,k,1000);b,a=zp2tf(zd,pd,kd);w=128;freqs(b,a,w)figure;freqz(b,a,w)1.取采樣頻率f=100Hz,用雙線性變換法設(shè)計(jì)五階Butterworth低通數(shù)字濾波器,繪出模擬濾波器與數(shù)字濾波器的幅頻與相頻特性z,p,k=besselap(5) ;zd,pd,kd=bilinear(z,p,k,100);b,a=zp2tf(zd,pd,kd);w=128;freqs(b,a,w)figure;freqz(b,a,w)2.設(shè)高通截止頻率為w0=10kHz, 取采樣頻率f=100,用雙線性變換法設(shè)計(jì)六階高通Butterworth數(shù)字濾波器,繪出數(shù)字濾波器的頻譜特性z,p,k=besselap(6) ;%原型低通b,a=zp2tf(z,p,k);%零極點(diǎn)型轉(zhuǎn)化為分式型w=128;w0=10000;bt,at=lp2hp(b,a,w0)%低通轉(zhuǎn)化為高通bz,az=bilinear(bt,at,1000);%雙線性法freqs(bt,at,w); %模擬頻譜figure;freqz(bz,az,w);%數(shù)字頻譜3.設(shè)帶通濾波器的濾波器中心頻率為W0=500KHz,帶寬為BW=10KHz, 取采樣頻率f=100,用脈沖相應(yīng)不變法設(shè)計(jì),設(shè)計(jì)五階帶通Butterworth數(shù)字濾波器,繪出數(shù)字濾波器的頻譜特性4.直接設(shè)計(jì)五階butterworth帶通濾波器,繪出頻譜圖。(高端與低端截止頻率分別為0.2和0.9)第7章 FIR濾波器設(shè)計(jì)第六章我們介紹了無(wú)限沖激響應(yīng)(IIR)濾波器的設(shè)計(jì)方法。其中最常用的由模擬濾波器轉(zhuǎn)換為數(shù)字濾波器的方法為雙線性變換法,因?yàn)檫@種方法無(wú)混疊效應(yīng),效果較好。但通過(guò)前面的例子我們看到,IIR數(shù)字濾波器相位特性不好(非線性,如圖 6-11、圖6-13、圖6-15 ),也不易控制。然而在現(xiàn)代信號(hào)處理中,例如圖像處理、數(shù)據(jù)傳輸、雷達(dá)接收以及一些要求較高的系統(tǒng)中對(duì)相位特性要求較為嚴(yán)格,這種濾波器就無(wú)能為力了。改善相位特性的方法是采用有限沖激響應(yīng)濾波器。本章首先對(duì)FIR濾波器原理及其使用函數(shù)作基本介紹,然后重點(diǎn)介紹窗函數(shù)法設(shè)計(jì)FIR濾波器,并對(duì)最優(yōu)濾波器設(shè)計(jì)函數(shù)進(jìn)行介紹。7.1 FIR濾波器原理概述及濾波函數(shù)7.1.1 FIR濾波器原理及設(shè)計(jì)方法分類根據(jù)第 6 章對(duì)數(shù)字濾波器的介紹,我們知道FIR濾波器的傳遞函數(shù)為: (7-1)可得FIR濾波器的系統(tǒng)差分方程為: 因此,F(xiàn)IR濾波器又稱為卷積濾波器。根據(jù)第 4 章中所描述的系統(tǒng)頻率響應(yīng),F(xiàn)IR濾波器的頻率響應(yīng)表達(dá)式為: (7-2)信號(hào)通過(guò)FIR濾波器不失真條件與(6-6)式所描述的相同,即濾波器在通帶內(nèi)具有恒定的幅頻特性和線性相位特性。理論上可以證明(這里從略):當(dāng)FIR濾波器的系數(shù)滿足下列中心對(duì)稱條件: (7-3)時(shí),濾波器設(shè)計(jì)在逼近平直幅頻特性的同時(shí),還能獲得嚴(yán)格的線性相位特性。線性相位FIR濾波器的相位滯后和群延遲在整個(gè)頻帶上是相等且不變的。對(duì)于一個(gè) N 階的線性相位FIR濾波器,群延遲為常數(shù),即濾波后的信號(hào)簡(jiǎn)單地延遲常數(shù)個(gè)時(shí)間步長(zhǎng)。這一特性使通帶頻率內(nèi)信號(hào)通過(guò)濾波器后仍保持原有波形形狀而無(wú)相位失真。本章主要介紹的FIR數(shù)字濾波器設(shè)計(jì)方法及 MATLAB 信號(hào)處理工具箱提供的FIR數(shù)字濾波器設(shè)計(jì)函數(shù),見表7-1。由于篇幅所限,本章我們主要介紹窗函數(shù)法和最優(yōu)化設(shè)計(jì)方法。表7-1 FIR濾波器設(shè)計(jì)的主要方法函數(shù)設(shè)計(jì)方法說(shuō)明工具函數(shù)窗函數(shù)法理想濾波器加窗處理fir1(單頻帶) , fir2(多頻帶) , kaiserord最優(yōu)化設(shè)計(jì)平方誤差最小化逼近理想幅頻響應(yīng)或Park-McClellan 算法產(chǎn)生等波紋濾波器firls , remez,remezord約束最小二乘逼近在滿足最大誤差限制條件下使整個(gè)頻帶平方誤差最小化fircls,fircls1升余弦函數(shù)具有光滑、正弦過(guò)渡帶的低通濾波器設(shè)計(jì)Fircos 7.1.2 FIR數(shù)字濾波器濾波函數(shù)相對(duì)于IIR 濾波器的濾波函數(shù),F(xiàn)IR數(shù)字濾波器濾波函數(shù)除了dimpulse和dstep僅適用于IIR濾波器外,其他各種函數(shù)可直接應(yīng)用于FIR濾波器,只是輸入的分母多項(xiàng)式向量a=1。另外,MATLAB還提供了一個(gè)函數(shù)fftfilt,該函數(shù)利用效率高的基于FFT算法實(shí)現(xiàn)對(duì)數(shù)據(jù)的濾波,該函數(shù)只適用于FIR濾波器,調(diào)用形式為:y=fftfilt(b,x,n)式中,b為FIR濾波器的系數(shù)向量;x為輸入數(shù)據(jù);n為FFT長(zhǎng)度,缺省時(shí),函數(shù)選用最佳的FFT長(zhǎng)度,y為濾波器的輸出。該函數(shù)執(zhí)行下面的操作:n=length(x);y=ifft(fft(x).*fft(b,n)./fft(a,n);應(yīng)注意,y=fftfilt(b,x)等價(jià)于y=filter(b,a,x)。7.2 FIR濾波器的窗函數(shù)設(shè)計(jì)7.2.1 窗函數(shù)的基本原理FIR濾波器設(shè)計(jì)的主要任務(wù)是根據(jù)給定的性能指標(biāo)確定濾波器的系數(shù)b,即系統(tǒng)單位脈沖序列h(n),它是一個(gè)有限長(zhǎng)序列。FIR濾波器的理想頻率響應(yīng),可寫成復(fù)數(shù)形式的Fourier級(jí)數(shù)形式: (7-4)式中,hd(n)是對(duì)應(yīng)的單位脈沖響應(yīng)序列。這說(shuō)明濾波器的頻率響應(yīng)和單位脈沖響應(yīng)互為Fourier變換對(duì)。因此其單位脈沖響應(yīng)可由下式求得, (7-5)求得序列后,通過(guò)z變換,可得到 (7-6)注意,這里為無(wú)限長(zhǎng)序列,因此是物理上不可實(shí)現(xiàn)的。如何變成物理上可實(shí)現(xiàn)呢?一個(gè)自然的想法是只取其中的某些項(xiàng),即只截取中的一部分,比如n=0,N-1,N為正整數(shù)。這種處理相當(dāng)于將,n=-與函數(shù)w(n)相乘,w(n)具有下列形式:w(n)相當(dāng)于一個(gè)矩形,我們稱之為矩形窗。即我們可采用矩形窗函數(shù)w(n)將無(wú)限脈沖響應(yīng)截取一段h(n)來(lái)近似為,這種截取在數(shù)學(xué)上表示為: h(n)= w(n) (7-7)這里應(yīng)該強(qiáng)調(diào)的是,加窗函數(shù)不是可有可無(wú)的,而是將設(shè)計(jì)變?yōu)槲锢砜蓪?shí)現(xiàn)所必須的。截取之后的濾波器傳遞函數(shù)變?yōu)椋?(7-8)式中,N為窗口寬度,H(z)是物理可實(shí)現(xiàn)系統(tǒng)。為了獲得線性相位,F(xiàn)IR濾波器h(n)必須滿足中心對(duì)稱條件(即7-3式),序列h(n)的延遲為。這種方法的基本原理是用一定寬度的矩形窗函數(shù)截取無(wú)限脈沖響應(yīng)序列獲得有限長(zhǎng)的脈沖響應(yīng)序列,從而得到FIR濾波器的脈沖響應(yīng),故稱為FIR濾波器的窗函數(shù)設(shè)計(jì)法。經(jīng)過(guò)加矩形窗后所得的濾波器實(shí)際頻率響應(yīng)能否很好地逼近理想頻率響應(yīng)呢?圖 7-1 示意給出了理想濾波器加矩形窗后的情況。理想低通濾波器的頻率響應(yīng)如圖中左上角圖,矩形窗的頻率響應(yīng)為左下角圖。時(shí)間域內(nèi)的乘積(7-7)式要求實(shí)際頻率響應(yīng)為這兩個(gè)頻率響應(yīng)函數(shù)在頻域內(nèi)的卷積(卷積定理),即得到圖形為圖7-1(右圖)。圖 7-1 FIR濾波器理想與實(shí)際頻率響應(yīng) 由圖可看出,加矩形窗后使實(shí)際頻率響應(yīng)偏離理想頻率響應(yīng),主要影響有三個(gè)方面:(1)理想幅頻特性陡直邊緣處形成過(guò)渡帶,過(guò)渡帶寬取決于矩形窗函數(shù)頻率響應(yīng)的主瓣寬度。(2)過(guò)渡帶兩側(cè)形成肩峰和波紋,這是矩形窗函數(shù)頻率響應(yīng)的旁瓣引起的,旁瓣相對(duì)值越大,旁瓣越多,波紋越多。(3)隨窗函數(shù)寬度N的增大,矩形窗函數(shù)頻率響應(yīng)的主瓣寬度減小,但不改變旁瓣的相對(duì)值。為了改善FIR濾波器性能,要求窗函數(shù)的主瓣寬度盡可能窄,以獲得較窄的過(guò)渡帶;旁瓣相對(duì)值盡可能小,數(shù)量盡可能少,以獲得通帶波紋小,阻帶衰減大,在通帶和阻帶內(nèi)均平穩(wěn)的特點(diǎn),這樣可使濾波器實(shí)際頻率響應(yīng)更好地逼近理想頻率響應(yīng)。這里我們明確兩個(gè)概念:截?cái)嗪皖l譜泄漏。信號(hào)是無(wú)限長(zhǎng)的,而在進(jìn)行信號(hào)處理時(shí)只能采取有限長(zhǎng)信號(hào),所以需要將信號(hào)“截?cái)唷?。在信?hào)處理中, “截?cái)唷北豢闯墒怯靡粋€(gè)有限長(zhǎng)的“窗口”看無(wú)限長(zhǎng)的信號(hào),或者從分析的角度是無(wú)限長(zhǎng)的信號(hào)x(t)乘以有限長(zhǎng)的窗函數(shù)w(t)。由傅立葉變換性質(zhì)可知,時(shí)間域內(nèi)的乘積對(duì)應(yīng)于頻率域的卷積,即 (7-9)這里,x(t)是頻寬有限信號(hào),而w(t)是頻寬無(wú)限信號(hào),表示互為Fourier變換對(duì)。截?cái)嗪蟮男盘?hào)也必須是頻寬無(wú)限信號(hào),這樣就是有限頻帶的信號(hào)分散到無(wú)限頻帶中去,這樣就產(chǎn)生了所謂頻譜泄漏。從能量的角度來(lái)看,頻譜泄漏也是能量的泄漏,因?yàn)榧哟昂笫乖瓉?lái)信號(hào)集中的窄頻帶內(nèi)的能量分散到無(wú)限的頻帶寬度范圍內(nèi)。頻譜泄漏是不可避免的,但要盡量減小。上邊只考慮了矩形窗,如果我們使窗的主瓣寬度盡可能地窄,旁瓣盡可能地小,可以獲得性能更好的濾波器,能否改變窗的形狀而達(dá)到這個(gè)目的呢?回答是肯定的。其實(shí)數(shù)字信號(hào)處理的前驅(qū)者們?cè)O(shè)計(jì)了不同于矩形窗的很多窗函數(shù),這些窗函數(shù)在主瓣和旁瓣特性方面各有特點(diǎn),可滿足不同的要求。為此,用窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器時(shí),要根據(jù)給定的濾波器性能指標(biāo)選擇窗口寬度N和窗函數(shù)w(n)。下面我們介紹窗函數(shù)。7.2.2 MATLAB信號(hào)處理中提供的窗函數(shù)(1)矩形窗:前面分析中所用的矩形窗可用下面函數(shù)來(lái)實(shí)現(xiàn)w=boxcar (N),N 為窗的長(zhǎng)度(以下函數(shù)與此同),w為返回的窗函數(shù)序列。(2)漢寧窗:w=hanning(N)漢寧窗的表達(dá)式為: (7-10)(3)哈明窗:w=hamming(N)哈明窗的表達(dá)式為: (7-11)(4)Bartlett窗:w=bartlett(N)Bartlett 窗的表達(dá)式為:當(dāng) N 為奇數(shù)時(shí), (7-12)當(dāng) N 為偶數(shù)時(shí), (7-13)(5) Blackman 窗:w= blackman(N)Blackman 窗的表達(dá)式為: (7-14)Blackman 窗比其他相同尺寸窗 (哈明窗,漢寧窗) 具有主瓣較寬和旁瓣泄漏較小的特點(diǎn)。(6)三角窗:w=triang(N)三角窗的表達(dá)式為:當(dāng) N 為奇數(shù)時(shí), (7-15)當(dāng) N 為偶數(shù)時(shí), (7-16)三角窗和Bartlett窗十分類似。三角窗的兩端值不為零,而Bartlett窗則為零,這一點(diǎn)可從例7-1中看出。(7)Kaiser窗:w=kaiser(n,beta)其中,beta是Kaiser窗參數(shù),影響窗旁瓣幅值的衰減率。Kaiser窗表達(dá)式: (7-17)式中, I0.是修正過(guò)的零階 Bessel 函數(shù)。Kaiser窗用于濾波器設(shè)計(jì)時(shí),若旁瓣幅值為,則 ( 7-18 )(8) Chebyshev窗:w=chebwin(n,r)式中, r 是窗口的旁瓣幅值在主瓣以下的分貝數(shù)。Chebyshev窗具有主瓣寬度最小,而旁瓣等高,高度可調(diào)整的特點(diǎn)。下面我們?cè)贛ATLAB觀看各種窗函數(shù)的形狀和頻率域圖象來(lái)驗(yàn)證上述所講特點(diǎn)。【例7-1】 用MATLAB編程繪制各種窗函數(shù)的形狀。窗函數(shù)的長(zhǎng)度為21。%Samp7_1clfNwin=21;n=0:Nwin-1; %數(shù)據(jù)總數(shù)和序列序號(hào)figure(1)for ii=1:4 switch ii case 1 w=boxcar(Nwin); %矩形窗 stext='矩形窗' case 2 w=hanning(Nwin); %漢寧窗 stext='漢寧窗' case 3 w=hamming(Nwin); %哈明窗 stext='哈明窗' case 4 w=bartlett(Nwin); %Bartlett窗 stext='Bartelett窗' end posplot='2,2,' int2str(ii); %指定繪制窗函數(shù)的圖形位置 subplot(posplot); stem(n,w); %繪出窗函數(shù) hold on plot (n ,w,'r'); %繪制包絡(luò)線 xlabel('n'); ylabel('w(n)'); title(stext); hold off; grid onendfigure(2)for ii=1:4 switch ii case 1 w=blackman(Nwin); %Blackman 窗 stext='Blackman窗' case 2 w=triang(Nwin); %三角窗 stext='三角窗' case 3 w=kaiser(Nwin,4); %Kaiser窗 stext='Kaiser窗(Beta=4)' case 4 w=chebwin(Nwin,40); %Chebyshev 窗 stext='Chebyshev窗(r=40)' end posplot='2,2,' int2str(ii); subplot(posplot); stem(n,w); %繪出窗函數(shù) hold on plot (n,w,'r'); %繪出包絡(luò)線 xlabel('n');ylabel('w(n)');title(stext); hold off;grid on;end程序運(yùn)行結(jié)果見圖 7-2 。可以看到各種窗函數(shù)的形狀。圖 7-2 各種窗函數(shù)的時(shí)間域形狀【例 7-2】 用 MATLAB 編程,采用512個(gè)頻率點(diǎn)繪制各種窗函數(shù)的幅頻特性。%Samp7_2clf;Nf=512; %窗函數(shù)復(fù)數(shù)頻率特性的數(shù)據(jù)點(diǎn)數(shù)Nwin=20; %窗函數(shù)數(shù)據(jù)長(zhǎng)度f(wàn)igure(1)for ii=1:4 switch ii case 1 w=boxcar(Nwin); %矩形窗 stext='矩形窗' case 2 w=hanning(Nwin); %漢寧窗 stext='漢寧窗' case 3 w=hamming (Nwin); %哈明窗 stext='哈明窗' case 4 w=bartlett(Nwin); % Bartlett窗 stext='Bartelett窗' end y,f=freqz(w,1,Nf); %求解窗函數(shù)的幅頻特性,窗函數(shù)相當(dāng)于一個(gè)數(shù)字濾波器 mag=abs(y);%求得窗函數(shù)幅頻特性 posplot='2,2,' int2str(ii); subplot(posplot); plot(f/pi,20* log10(mag/max(mag); %繪制窗函數(shù)的幅頻特性 xlabel('歸一化頻率');ylabel('振幅/dB'); title(stext);grid on;endfigure(2)for ii=1:4 switch ii case 1 w=blackman(Nwin); %Blackman 窗 stext='Blackman窗' case 2 w=triang(Nwin); %三角窗 stext='三角窗' case 3 w=kaiser(Nwin,4); %Kaiser窗 stext='Kaiser窗(Beta=4)' case 4 w=chebwin(Nwin,40); %Chebyshev 窗 stext='Chebyshev窗(r=40)' end y,f=freqz(w,1,Nf); %以 Nf點(diǎn)數(shù)求解窗函數(shù)的幅頻響應(yīng) mag=abs(y); %求得窗函數(shù)幅頻響應(yīng) posplot='2,2,' int2str(ii); subplot(posplot);plot(f/pi,20* log10(mag/max(mag); %繪制幅頻響應(yīng) xlabel('歸一化頻率');ylabel('振幅/dB'); title(stext);grid on;end程序運(yùn)行結(jié)果見圖 7-3 ??梢钥吹礁鞣N窗函數(shù)的幅頻形狀。對(duì)照該圖可知這些窗函數(shù)具有上面所分析的窗函數(shù)的特征。圖 7-3 各種窗函數(shù)的幅頻形狀 由圖 7-3 可見,各種窗函數(shù)都具有明顯的主瓣(Mainlobe)和旁瓣(Sidelobe)。主瓣頻寬和旁瓣的幅值衰減特性決定了窗函數(shù)的應(yīng)用場(chǎng)合。矩形窗具有最窄的主瓣,但也有最大的旁瓣峰值(第一旁瓣衰減為 13 dB);Blackman窗具有最大的旁瓣衰減,但也具有最寬的主瓣寬度。不同窗函數(shù)在這兩方面的特點(diǎn)是不同的,因此應(yīng)根據(jù)具體的問(wèn)題進(jìn)行選擇。通常來(lái)講,哈明窗和漢寧窗的主瓣具有較小的旁瓣和較大的衰減速度,是較為常用的窗函數(shù)。表7-2總結(jié)了各種窗函數(shù)主瓣和旁瓣的特征(理論分析可參考其他的數(shù)字信號(hào)處理教材),大家可對(duì)照窗函數(shù)的幅頻形狀(圖7-3)認(rèn)真理解體會(huì)。表7-2 各種窗函數(shù)的特點(diǎn)窗函數(shù)主瓣寬第一旁瓣相對(duì)主瓣衰減(分貝)矩形窗-13漢寧窗-31哈明窗-41Bartlett窗-25Blackman 窗-57三角窗-25Kaiser窗可調(diào)整可調(diào)整Chebyshev 窗可調(diào)整可調(diào)整 主旁瓣頻率寬度還與窗函數(shù)長(zhǎng)度N有關(guān)。增加窗函數(shù)長(zhǎng)度N將減小窗函數(shù)的主瓣寬度,但不能減小旁瓣幅值衰減的相對(duì)值(分貝數(shù)),這個(gè)值是由窗函數(shù)決定的。這個(gè)特點(diǎn)可由下面的例子清楚地看出?!纠?-3】繪制矩形窗函數(shù)的幅頻響應(yīng),窗長(zhǎng)度分別為:(1)N=10;(2)N=20; (3)N=50;(4)N=100.%Samp7_3clf;Nf=512;for ii=1:4 switch ii case 1 Nwin=10; case 2 Nwin=20; case 3 Nwin=50; case 4 Nwin=100; end w=boxcar(Nwin); %矩形窗 y,f=freqz(w,1,Nf); %用不同的窗長(zhǎng)度求得復(fù)數(shù)頻率特性 mag=abs(y); %求得幅頻特性 posplot='2,2,' int2str(ii); %指定繪圖位置 subplot(posplot); plot (f/pi,20*log10(mag/max(mag); %繪出幅頻形狀 xlabel('歸一化頻率');ylabel('振幅/dB'); stext='N=' int2str(Nwin); %給出標(biāo)題,指出所用的數(shù)據(jù)個(gè)數(shù) title(stext);grid on;end 圖 7-4 數(shù)據(jù)長(zhǎng)度不同的矩形窗的幅頻形狀 程序運(yùn)行結(jié)果見圖7-4。顯然,隨著N的增大,主瓣和旁瓣都變窄,但第一旁瓣相對(duì)主瓣的幅值下降分貝數(shù)相同,第二旁瓣相對(duì)第一旁瓣幅值下降的分貝數(shù)也相同。這樣,隨著N的增大,旁瓣也得到抑制,有力地減少了頻譜泄漏,但不能完全消除。減少主瓣寬度和抑制旁瓣是一對(duì)矛盾,不可兼得,只能根據(jù)不同用途折衷處理。7.2.3 運(yùn)用窗函數(shù)設(shè)計(jì)數(shù)字濾波器用于信號(hào)分析中的窗函數(shù)可根據(jù)用戶的不同要求選擇。用于濾波器的窗函數(shù),一般要求窗函數(shù)的主瓣寬度窄,以獲得較好的過(guò)渡帶;旁瓣相對(duì)值盡可能少,增加通帶的平穩(wěn)度和增大阻帶的衰減?;诖昂瘮?shù)的FIR數(shù)字濾波器設(shè)計(jì)的算法十分簡(jiǎn)單,其主要步驟為:(1)對(duì)濾波器理想頻域幅值響應(yīng)進(jìn)行傅立葉逆變換獲得理想濾波器的單位脈沖響應(yīng)hd(n)。一般假定理想低通濾波器的截止頻率為,其幅頻特性滿足 (7-19)則根據(jù)傅立葉逆變換,單位脈沖響應(yīng)為: (7-20)其中,為信號(hào)延遲。(2)由性能指標(biāo)(阻帶衰減的分貝數(shù))根據(jù)表7-2第3列的值確定滿足阻帶衰減的窗函數(shù)類型w(n)。濾波器的階數(shù)越高,濾波器的幅頻特性越好,但數(shù)據(jù)處理的費(fèi)用也越高,因此像IIR濾波器一樣,F(xiàn)IR濾波器也要確定滿足性能指標(biāo)的濾波器最小階數(shù)。由前面的討論(圖7-1)可知,濾波器的主瓣寬度相當(dāng)于過(guò)渡帶寬,因此,使過(guò)渡帶寬近似于窗函數(shù)主瓣寬(表7-2中的第二列)可求得滿足性能指標(biāo)的窗口長(zhǎng)度N,此時(shí),信號(hào)延遲為(N-1)/2。(3)求實(shí)際濾波器的單位脈沖響應(yīng)h(n):根據(jù)h(n)=hd(n)*w(n)。(4)檢驗(yàn)濾波器的性能。可設(shè)定一些信號(hào)采用 7.1.2 節(jié)指出的函數(shù)或6.3.2節(jié)所給的函數(shù)進(jìn)行濾波。下面采用實(shí)例說(shuō)明如何根據(jù)上面步驟設(shè)計(jì)FIR濾波器。【例 7-4】 用窗函數(shù)設(shè)計(jì)一個(gè)線性相位FIR低通濾波器,并滿足性能指標(biāo):通帶邊界的歸一化頻率wp=0.5,阻帶邊界的歸一化頻率ws=0.66,阻帶衰減不小于30dB,通帶波紋不大于3dB。假設(shè)一個(gè)信號(hào),其中f1=5Hz,f2=20Hz。信號(hào)的采樣頻率為50Hz。試將原信號(hào)與通過(guò)濾波器的信號(hào)進(jìn)行比較。由題意,阻帶衰減不小于30dB,根據(jù)表7-2,選取漢寧窗,因?yàn)闈h寧窗的第一旁瓣相對(duì)主瓣衰減為31dB,滿足濾波要求。在窗函數(shù)設(shè)計(jì)法中,要求設(shè)計(jì)的頻率歸一化到0區(qū)間內(nèi),Nyquist頻率對(duì)應(yīng)于,因此通帶和阻帶邊界頻率為0.5和0.66。程序如下%Samp7_4wp=0.5*pi;ws=0.66*pi; %濾波器邊界頻率wdelta=ws-wp; %過(guò)渡帶寬N=ceil(8*pi/wdelta) %根據(jù)過(guò)渡帶寬等于表 7-2中漢寧窗函數(shù)主瓣寬求得濾波器所用窗函數(shù)的最小長(zhǎng)度Nw=N;wc=(wp+ws)/2; %截止頻率在通帶和阻帶邊界頻率的中點(diǎn)n=0:N-1;alpha=(N-1)/2; %求濾波器的相位延遲m=n-alpha+eps; %eps為MATLAB系統(tǒng)的精度hd=sin(wc*m)./(pi*m); %由(7-20)式求理想濾波器脈沖響應(yīng)win=hanning(Nw); %采用漢寧窗h=hd.*win' %在時(shí)間域乘積對(duì)應(yīng)于頻率域的卷積b=h; figure(1)H,f=freqz(b,1,512,50); %采用 50 Hz 的采樣頻率繪出該濾波器的幅頻和相頻響應(yīng)subplot(2,1,1),plot(f,20*log10(abs(H)xlabel('頻率/Hz');ylabel('振幅/dB');grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel('頻率/Hz');ylabel('相位/o');grid on;%impz(b,1); %可采用此函數(shù)給出濾波器的脈沖響應(yīng)%zplane(b,1); %可采用此語(yǔ)句給出濾波器的零極點(diǎn)圖%grpdelay(b,1); %可采用此函數(shù)給出濾波器的群延遲f1=3;f2=20; %檢測(cè)輸入信號(hào)含有兩種頻率成分dt=0.02; t=0:dt:3; %采樣間隔和檢測(cè)信號(hào)的時(shí)間序列x=sin(2*pi*f1* t)+cos(2* pi*f2* t); %檢測(cè)信號(hào)%y=filter(b,1,x); %可采用此函數(shù)給出濾波器的輸出y=fftfilt(b,x); %給出濾波器的輸出figure(2)subplot(2,1,1), plot(t,x),title('輸入信號(hào)') %繪輸入信號(hào)subplot(2,1,2),plot(t,y) % 繪輸出信號(hào)hold on; plot(1 1*(N-1)/2*dt,ylim, 'r') %繪出延遲到的時(shí)刻xlabel('時(shí)間/s'),title('輸出信號(hào)')程序運(yùn)行結(jié)果見圖7-5和圖7-6。該例設(shè)計(jì)通帶邊界wp=0.5,阻帶邊界頻率ws=0.66,對(duì)應(yīng)于50Hz的采樣頻率通帶邊界頻率為fp=Fs/2*Fnormal=50/2*0.5=12.5Hz, fs=50/2*0.66=16.5Hz, 其中Fs為采樣頻率,F(xiàn)normal為歸一化頻率。由圖7-5上圖可以看到,在小于12.5Hz的頻段上,幾乎看不到下降,即滿足通帶波紋不大于3dB的要求。在大于16.5Hz的頻段上,阻帶衰減大于30dB,滿足設(shè)計(jì)要求。由圖7-5下圖可見,在通帶范圍內(nèi),相位頻率為一條直線,表明該濾