實(shí)驗(yàn)四窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器.doc

《實(shí)驗(yàn)四窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器.doc》由會(huì)員分享，可在線閱讀，更多相關(guān)《實(shí)驗(yàn)四窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器.doc（16頁(yè)珍藏版）》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。

實(shí)驗(yàn)四 窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器 一、實(shí)驗(yàn)?zāi)康? 1、掌握窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器的原理及具體方法。 2、掌握頻率取樣法設(shè)計(jì)FIR數(shù)字濾波器的原理和基本方法。 3、學(xué)習(xí)利用窗函數(shù)法和頻率取樣法設(shè)計(jì)低通、帶通、高通、帶阻數(shù)字濾波器。 二、實(shí)驗(yàn)環(huán)境 計(jì)算機(jī)、MATLAB軟件 三、實(shí)驗(yàn)基礎(chǔ)理論 窗函數(shù)設(shè)計(jì)FIR濾波器 1.基本原理 窗函數(shù)設(shè)計(jì)法的基本思想為，首先選擇一個(gè)適當(dāng)?shù)睦硐氲臑V波器，然后用窗函數(shù)截取它的單位脈沖響應(yīng)，得到線性相位和因果的FIR濾波器。這種方法的重點(diǎn)是選擇一個(gè)合適的窗函數(shù)和理想濾波器，使設(shè)計(jì)的濾波器的單位脈沖響應(yīng)逼近理想濾波器的單位脈沖響應(yīng)。 2.設(shè)計(jì)步驟 （1）給定理想濾波器的頻率響應(yīng)，在通帶上具有單位增益和線性相位，在阻帶上具有零響應(yīng)。一個(gè)帶寬為的低通濾波器由下式給定： 其中為采樣延遲，其作用是為了得到一個(gè)因果系統(tǒng)。 （2）確定這個(gè)濾波器的單位脈沖響應(yīng) 為了得到一個(gè)長(zhǎng)度為N的因果的線性相位FIR濾波器，我們令 （3）用窗函數(shù)截取得到所設(shè)計(jì)FIR數(shù)字濾波器： 3.窗函數(shù)的選擇 常用的窗函數(shù)有矩形（Rectangular）窗，漢寧（Hanning）窗，海明（Hamming）窗、布萊克曼（Blackman）窗、凱瑟（Kaiser）窗等 表4-1 MATLAB中產(chǎn)生窗函數(shù)的命令 MATLAB函數(shù) 窗函數(shù) MATLAB函數(shù) 窗函數(shù) Boxcar 矩形窗函數(shù) Blackman 布萊克曼窗 Hanning 漢寧窗函數(shù) Kaiser 凱瑟窗函數(shù) Hamming 海明窗 表4-2 常用窗函數(shù)的特性 窗函數(shù) 窗函數(shù)頻率特性 加窗后濾波器指標(biāo) 旁瓣峰值dB 主瓣寬度 過(guò)渡帶寬 最小阻帶衰減dB 矩形窗 -13 4π/N 1.8π/N -21 漢寧窗 -31 8π/N 6.2π/N -44 海明窗 -41 8π/N 6.6π/N -53 布萊克曼窗 -57 12π/N 11π/N -74 凱瑟窗是一種廣泛在實(shí)際中廣泛應(yīng)用的窗函數(shù)，它由下式給定： 其中是修正的零階貝塞爾函數(shù)，參數(shù)控制最小阻帶衰減，這種窗函數(shù)對(duì)于相同的N可以提供不同的過(guò)渡帶寬。由于貝塞爾函數(shù)比較復(fù)雜，這種窗函數(shù)的設(shè)計(jì)方程很難推導(dǎo)，然而幸運(yùn)的是，有一些經(jīng)驗(yàn)設(shè)計(jì)方程可以直接使用。 已知給定的指標(biāo) ，濾波器長(zhǎng)度N和凱瑟窗參數(shù)可以按如下凱瑟窗方程給出 過(guò)渡帶帶寬： 頻率取樣設(shè)計(jì)FIR濾波器 1.基本原理 頻率取樣法從頻域出發(fā)，把理想的濾波器等間隔采樣得到,將作為實(shí)際設(shè)計(jì)濾波器的： 得到以后可以由來(lái)確定唯一確定濾波器的單位脈沖響應(yīng),可以由求得： 其中為內(nèi)插函數(shù)： 有求得的頻率響應(yīng)將逼近。 如果我們?cè)O(shè)計(jì)的是線性相位FIR濾波器，則的幅度和相位滿足線性相位濾波器的約束條件。 我們將表示為如下形式 當(dāng)為實(shí)數(shù)，則 由此得到 即為中心偶對(duì)稱。在利用線性相位條件可知，對(duì)于1型和2型線性相位濾波器： 對(duì)于3型和4型線性相位濾波器 2.設(shè)計(jì)步驟 （1）由給定的理想濾波器給出和。 （2）由求得 （3）根據(jù)求得或 四、實(shí)驗(yàn)內(nèi)容 1、設(shè)計(jì)一個(gè)數(shù)字低通FIR濾波器，其技術(shù)指標(biāo)如下： 分別采用矩形窗、漢寧窗、海明窗、布萊克曼窗、凱瑟窗設(shè)計(jì)該濾波器。結(jié)合實(shí)驗(yàn)結(jié)果，分別討論采用上述方法設(shè)計(jì)的數(shù)字濾波器是否都能滿足給定指標(biāo)要求。 （1） 矩形窗 程序代碼： wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp; N=ceil(1.8*pi/tr_width) n=0:N-1; wc=(wst+wp)/2; alpha=(N-1)/2; hd=(wc/pi)*sinc((wc/pi)*(n-alpha)); w_boxcar=boxcar(N); h=hd.*w_boxcar; subplot(221); stem(n,hd,filled); axis tight;xlabel(n);ylabel(hd(n)); [Hr,w1]=zerophase(h); subplot(222); plot(w1/pi,Hr); axis;xlabel(\omega/\pi);ylabel(H(\omega)); subplot(223); stem(n,h,filled); axis tight;xlabel(n);ylabel(h(n)); [H,w]=freqz(h,1); subplot(224); plot(w/pi,20*log10(abs(H)/max(H))); axis tight;xlabel(\omega/\pi);ylabel(dB); grid on; MATLAB圖形： （2）漢寧窗 程序代碼： wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp; N=ceil(6.2*pi/tr_width) n=0:N-1; wc=(wst+wp)/2; alpha=(N-1)/2; hd=(wc/pi)*sinc((wc/pi)*(n-alpha)); w_boxcar=hanning(N); h=hd.*w_boxcar; subplot(221); stem(n,hd,filled); axis tight;xlabel(n);ylabel(hd(n)); [Hr,w1]=zerophase(h); subplot(222); plot(w1/pi,Hr); axis;xlabel(\omega/\pi);ylabel(H(\omega)); subplot(223); stem(n,h,filled); axis tight;xlabel(n);ylabel(h(n)); [H,w]=freqz(h,1); subplot(224); plot(w/pi,20*log10(abs(H)/max(H))); axis tight;xlabel(\omega/\pi);ylabel(dB); grid on; MATLAB圖形： （3）海明窗 程序代碼： wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp; N=ceil(6.6*pi/tr_width) n=0:N-1; wc=(wst+wp)/2;alpha=(N-1)/2; hd=(wc/pi)*sinc((wc/pi)*(n-alpha)); w_boxcar=hamming(N); h=hd.*w_boxcar; subplot(221); stem(n,hd,filled); axis tight;xlabel(n);ylabel(hd(n)); [Hr,w1]=zerophase(h); subplot(222); plot(w1/pi,Hr); axis;xlabel(\omega/\pi);ylabel(H(\omega)); subplot(223); stem(n,h,filled); axis tight; xlabel(n); ylabel(h(n)); [H,w]=freqz(h,1); subplot(224); plot(w/pi,20*log10(abs(H)/max(H))); axis tight;xlabel(\omega/\pi);ylabel(dB); grid on; MATLAB圖形： （4）布萊克曼窗 程序代碼： wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp; N=ceil(11*pi/tr_width) n=0:N-1; wc=(wst+wp)/2;alpha=(N-1)/2; hd=(wc/pi)*sinc((wc/pi)*(n-alpha)); w_boxcar=blackman(N); h=hd.*w_boxcar; subplot(221); stem(n,hd,filled); axis tight;xlabel(n);ylabel(hd(n)); [Hr,w1]=zerophase(h); subplot(222); plot(w1/pi,Hr); axis;xlabel(\omega/\pi);ylabel(H(\omega)); subplot(223); stem(n,h,filled); axis tight;xlabel(n);ylabel(h(n)); [H,w]=freqz(h,1); subplot(224); plot(w/pi,20*log10(abs(H)/max(H))); axis tight;xlabel(\omega/\pi);ylabel(dB); grid on; MATLAB圖形為： （5）凱瑟窗 程序代碼： wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp;As=50; N=ceil((As-7.95)/(2.285*tr_width))+1; beta=0.1102*(As-8.7); n=0:N-1; wc=(wst+wp)/2;alpha=(N-1)/2; hd=(wc/pi)*sinc((wc/pi)*(n-alpha)); w_boxcar=kaiser(N,beta); h=hd.*w_boxcar; subplot(221); stem(n,hd,filled); axis tight;xlabel(n);ylabel(hd(n)); [Hr,w1]=zerophase(h); subplot(222); plot(w1/pi,Hr); axis;xlabel(\omega/\pi);ylabel(H(\omega)); subplot(223); stem(n,h,filled); axis tight;xlabel(n);ylabel(h(n)); [H,w]=freqz(h,1); subplot(224); plot(w/pi,20*log10(abs(H)/max(H))); axis tight;xlabel(\omega/\pi);ylabel(dB); grid on; MATLAB圖形： 2、設(shè)計(jì)一個(gè)數(shù)字帶通FIR濾波器，其技術(shù)指標(biāo)如下： 下阻帶邊緣： 下通帶邊緣： 上通帶邊緣： 上阻帶邊緣： 程序代碼： wp1=0.2*pi;Rp1=1; wst1=0.35*pi;A1=60; width1=wst1-wp1; N1=ceil(11*pi/width1)+1; n1=0:(N1-1); wc1=(wp1+wst1)/2; alpha=(N1-1)/2; wp2=0.65*pi;Rp2=1;wst2=0.8*pi;A2=60; width2=wst2-wp2; N2=ceil(11*pi/width2)+1; n2=0:(N2-1); wc2=(wp2+wst2)/2; alpha=(N2-1)/2; hd=(wc2/pi)*sinc((wc2/pi)*(n2-alpha))-(wc1/pi)*sinc((wc1/pi)*(n1-alpha)); w_w=blackman(N1); h=hd.*w_w; subplot(221); stem(n1,h,filled); subplot(222); [H,w]=freqz(h,1); plot(w/pi,20*log10(abs(H)/max(abs(H)))); subplot(223); [Hr,w1]=zerophase(h); plot(w1/pi,Hr); subplot(224); stem(n1,hd,filled); [Hr,wl]=zerophase(h); grid on; MATLAB圖形： 3.采用頻率取樣法設(shè)計(jì)FIR數(shù)字低通濾波器，滿足以下指標(biāo) （1）取N=20，過(guò)渡帶沒(méi)有樣本。 （2）取N=40，過(guò)渡帶有一個(gè)樣本，T=0.39。 （3）取N=60，過(guò)渡帶有兩個(gè)樣本，T1=0.5925，T2=0.1099。 （4）分別討論采用上述方法設(shè)計(jì)的數(shù)字濾波器是否都能滿足給定的指標(biāo)要求。 （1）程序代碼： N=20; alpha=(N-1)/2; L=0:N-1;wL=(2*pi/N)*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((N-1)/2); k2=floor((N-1)/2)+1:N-1; angH=[-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2)]; H=Hrs.*exp(j*angH); h=ifft(H,N); w=[0:500]*pi/500; [Hr,wr]=zerophase(h); subplot(221); plot(wdL,Hdr,wL(1:N/2+1)/pi,Hrs(N/2+1)); axis([0,1,-0.1,1.1]); xlabel(\omega (\pi)); ylabel(Hr(k)); subplot(222); stem(L,h,filled); axis([0,N-1,-0.1,0.3]); xlabel(n); ylabel(h(n)); subplot(223); plot(wr/pi,Hr,wL(1:N/2+1)/pi,Hrs(1:N/2+1)); axis([0,1,-0.2,1.2]); grid on; xlabel(\omega (\pi)); ylabel(Hr(\omega)); subplot(224); plot(wr/pi,20*log10((abs(Hr)/max(abs(Hr))))); axis([0,1,-50,5]); grid on;xlabel(\omega (\pi));ylabel(dB) MATLAB圖形如下： （2）程序代碼： N=40; alpha=(N-1)/2; L=0:N-1; wL=(2*pi/N)*L; Hrs=[1,1,1,1,1,0.39,zeros(1,29),0.39,1,1,1,1]; Hdr=[1,1,0.39,0,0]; wdL=[0,0.2,0.25,0.3,1]; k1=0:floor((N-1)/2); k2=floor((N-1)/2)+1:N-1; angH=[-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2)]; H=Hrs.*exp(j*angH); h=ifft(H,N); w=[0:500]*pi/500; [Hr,wr]=zerophase(h); subplot(221); plot(wdL,Hdr,wL(1:N/2+1)/pi,Hrs(N/2+1)); axis([0,1,-0.1,1.1]); xlabel(\omega (\pi)); ylabel(Hr(k)); subplot(222); stem(L,h,filled); axis([0,N-1,-0.1,0.3]); xlabel(n); ylabel(h(n)); subplot(223); plot(wr/pi,Hr,wL(1:N/2+1)/pi,Hrs(1:N/2+1)); axis([0,1,-0.2,1.2]); grid on; xlabel(\omega (\pi)); ylabel(Hr(\omega)); subplot(224); plot(wr/pi,20*log10((abs(Hr)/max(abs(Hr))))); axis([0,1,-50,5]); grid on; xlabel(\omega (\pi)); ylabel(dB) grid on; MATLAB圖形如下： 3、程序代碼： N=60; alpha=(N-1)/2; L=0:N-1; wL=(2*pi/N)*L; Hrs=[1,1,1,1,1,1,1,0.5925,0.1099,zeros(1,43),0.1099,0.5925,1,1,1,1,1,1]; Hdr=[1,1,0.5925,0.1099,0,0]; wdL=[0,0.2,7/30,8/30,0.3,1]; k1=0:floor((N-1)/2); k2=floor((N-1)/2)+1:N-1; angH=[-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2)]; H=Hrs.*exp(j*angH); h=ifft(H,N);w=[0:500]*pi/500; [Hr,wr]=zerophase(h); subplot(221); plot(wdL,Hdr,wL(1:N/2+1)/pi,Hrs(N/2+1)); axis([0,1,-0.1,1.1]); xlabel(\omega (\pi)); ylabel(Hr(k)); subplot(222); stem(L,h,filled); axis([0,N-1,-0.1,0.3]); xlabel(n); ylabel(h(n)); subplot(223); plot(wr/pi,Hr,wL(1:N/2+1)/pi,Hrs(1:N/2+1)); axis([0,1,-0.2,1.2]); grid on; xlabel(\omega (\pi)); ylabel(Hr(\omega)); subplot(224); plot(wr/pi,20*log10((abs(Hr)/max(abs(Hr))))); axis([0,1,-50,5]);grid on;xlabel(\omega (\pi));ylabel(dB) grid on; MATLAB圖形: 各階數(shù)的通阻帶指標(biāo)如下表： 階數(shù) 通帶波動(dòng)0.25dB 阻帶衰減50dB 是否滿足指標(biāo) N=20 不滿足 不滿足 否 N=40 不滿足 滿足 否 N=60 滿足 滿足 是 4.采用頻率取樣技術(shù)設(shè)計(jì)下面的高通濾波器 對(duì)于高通濾波器，N必須為奇數(shù)（或1型濾波器）。選擇N=33，過(guò)渡帶有兩個(gè)樣本，過(guò)渡帶的最優(yōu)值為T1=0.1095，T2=0.598. 程序代碼： N=33;alpha=(N-1)/2;L=0:N-1;wL=(2*pi/N)*L; Hrs=[0,0,0,0,0,0,0,0,0,0,0,0.1095,0.598,1,1,1,1,1,1,1,0.598,0.1095,0,0,0,0,0,0,0,0,0,0,0]; Hdr=[0,0,0.1095,0.598,1,1]; wdL=[0,0.6,2/3,24/33,0.8,1]; k1=0:floor((N-1)/2); k2=floor((N-1)/2)+1:N-1; angH=[-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2)]; H=Hrs.*exp(j*angH); h=ifft(H,N); w=[0:500]*pi/500; Hr=h*cos((alpha-L)*w); subplot(221); plot(wdL,Hdr,wL(1:17)/pi,Hrs(1:17)); axis([0,1,-0.1,1.1]); xlabel(\omega (\pi)); ylabel(Hr(k)); subplot(222); stem(L,h,filled); axis([0,N-1,-0.1,0.3]); xlabel(n); ylabel(h(n)); subplot(223); plot(w/pi,Hr,wL(1:17)/pi,Hrs(1:17)); axis([0,1,-0.2,1.2]); xlabel(\omega (\pi)); ylabel(Hr(\omega)); subplot(224); plot(w/pi,20*log10((abs(Hr)/max(abs(Hr))))); axis([0,1,-50,5]);xlabel(\omega (\pi));ylabel(dB) grid on; MATLAB圖形為： 五、實(shí)驗(yàn)心得與體會(huì) 通過(guò)這次實(shí)驗(yàn)，學(xué)習(xí)了運(yùn)用窗函數(shù)法和頻率取樣法設(shè)計(jì)具有線性相位的數(shù)字低通、高通、帶通、帶阻濾波器。同時(shí)也從圖像上加深了對(duì)線性相位的理解和認(rèn)識(shí)。