數(shù)字信號(hào)處理數(shù)字濾波器的設(shè)計(jì)

1、 數(shù)字信號(hào)處理1 fir（有限脈沖響應(yīng)）數(shù)字濾波器設(shè)計(jì)1.1設(shè)計(jì)題目：用漢寧窗設(shè)計(jì)線性相位fir低通數(shù)字濾波器。1.2設(shè)計(jì)原理：窗函數(shù)設(shè)計(jì)的基本思想是用firdf（有限長(zhǎng)單位脈沖響應(yīng)數(shù)字濾波器）逼近希望濾波特性。設(shè)希望逼近的濾波器的頻率響應(yīng)函數(shù)為，其脈沖響應(yīng)用表示。為了計(jì)算簡(jiǎn)便方便，通常選擇為具有片段常數(shù)特性的理想濾波器。因此是無(wú)限長(zhǎng)非因果序列，不能直接作為firdf的單位脈沖響應(yīng)。窗函數(shù)設(shè)計(jì)法就是截取為有限長(zhǎng)的一段因果序列，并且合適的窗函數(shù)進(jìn)行加權(quán)作為firdf的單位脈沖響應(yīng)。后面看到，截取的長(zhǎng)度和加權(quán)窗函數(shù)的類(lèi)型都直接影響逼近精度（濾波器指標(biāo)）。具體的設(shè)計(jì)步驟如下： 1） 構(gòu)造希望逼近的頻

2、率響應(yīng)函數(shù)。以低通線性相位firdf設(shè)計(jì)為例，一般選擇為線性相位理想低通濾波器，即， ， 2 ）求出。對(duì) 進(jìn)行ift得到3 ）加窗得到firdf的單位脈沖響應(yīng) 式中，成為窗函數(shù)，其長(zhǎng)度為n。1.3窗函數(shù)：1.3.1 漢寧窗（hanning） w=hanning（n）1.3.2漢寧窗表達(dá)式， k=1,2.n1.3.3漢寧窗w的繪圖程序:n=21;n=0:n-1; %數(shù)據(jù)總數(shù)和序列序號(hào)figure(1) w=hanning(n); %漢寧窗 posplot=2,2, int2str(i); %指定繪制窗函數(shù)的圖形位置 stem(n,w); %繪出窗函數(shù) hold on xlabel(n); yla

3、bel(w(n); title(漢寧窗); hold off; grid on；1.3.4漢寧窗的仿真波形： 1.4設(shè)計(jì)指標(biāo)：通帶截止頻 率=0.4 rad，阻帶截止頻率=0.6 rad，阻帶最小衰減=40db，1.5設(shè)計(jì)過(guò)程：1） 選擇窗函數(shù)，計(jì)算窗函數(shù)長(zhǎng)度n。由阻帶衰減40db知道，hanning和hamming窗 函數(shù)都滿足要求，這兩種窗函數(shù)過(guò)渡帶寬度相同，選擇hanning窗。hanning窗函數(shù)的過(guò)渡帶寬度為8/n,所以應(yīng)當(dāng)滿足： 8/nb；解得n8/b=402）確定希望逼近的理想低通濾波器截止頻率（該點(diǎn)為過(guò)渡中心）處濾波器頻響幅度為0.5（等價(jià)于衰減6db），所以 =0.53）計(jì)算

4、4)加窗得到設(shè)計(jì)結(jié)果 1.6 加窗后h（n）及損耗的繪圖源程序：wp=0.4*pi;ws=0.6*pi;b=ws-wp;n=ceil(8*pi/b)wc=(wp+ws)/2/pi;hn=fir1(n-1,wc) m=1024;hk=fft(hn,m);n=0:n-1;subplot(1,2,1);stem(n,hn,.);xlabel(n);ylabel(h(n);grid onk=1:m/2+1;w=2*(0:m/2)/m;subplot(1,2,2);plot(w,20*log10(abs(hk(k);axis(0,1,-100,5);xlabel(w/pi);ylabel(20lg|hg

5、(w)|);1.7 運(yùn)行結(jié)果：n =40hn = columns 1 through 12 -0.0009 -0.0010 0.0013 0.0018 -0.0025 -0.0034 0.0047 0.0062 -0.0081 -0.0104 0.0133 0.0168 -0.0211 -0.0267 0.0341 0.0443 -0.0598 -0.0869 0.1483 0.4503 0.4503 0.1483 -0.0869 -0.0598 0.0443 0.0341 -0.0267 -0.0211 0.0168 0.0133 -0.0104 -0.0081 0.0062 0.0047

6、-0.0034 -0.0025 0.0018 0.0013 -0.0010 -0.00091.8 加窗后h（n）及損耗的仿真波形由圖顯示：1)由于加窗使得，所以，存在誤差。 2）用漢寧窗進(jìn)行數(shù)字濾波器的設(shè)計(jì)簡(jiǎn)單，實(shí)用，但濾波器的邊界頻率不易精確控制。所設(shè)計(jì)的濾波器在阻帶邊界頻率附近的衰減最小，距阻帶邊界頻率越遠(yuǎn)，衰減越大。因此所設(shè)計(jì)的濾波器性?xún)r(jià)比較低。2.iir（無(wú)限脈沖響應(yīng)）模擬濾波器設(shè)計(jì)2.1 設(shè)計(jì)題目：設(shè)計(jì)一個(gè)切比雪夫i型模擬濾波器2.2設(shè)計(jì)原理：模擬濾波器的設(shè)計(jì)是通過(guò)查表或計(jì)算確定濾波器的系統(tǒng)函數(shù)，再進(jìn)行安裝調(diào)試，一般的設(shè)計(jì)步驟如下：1） 根據(jù)信號(hào)處理要求確定設(shè)計(jì)指標(biāo)2） 選擇濾波器的

7、類(lèi)型3） 計(jì)算濾波器的階數(shù)n4） 通過(guò)查表或計(jì)算確定濾波器的系統(tǒng)函數(shù)5） 綜合實(shí)現(xiàn)并安裝調(diào)試2.3技術(shù)指標(biāo)如下： 通帶邊界頻率wp=0.2*pi，通帶最小衰減rp=1阻帶截止頻率wr=0.4*pi，阻帶最小衰減rs=50db2.4 設(shè)計(jì)源程序：% 主程序：wp=0.2*pi;wr=0.4*pi;rp=1;rs=50; b,a=afd_chb1(wp,wr,rp,rs) c,b,a=sdir2cas(b,a)db,mag,pha,w=freqs_m(b,a,pi);ha,x,t=impulse(b,a);subplot(2,2,1);plot(w/pi,mag);title(幅度響應(yīng)|ha(jo

8、mega)|);subplot(2,2,2);plot(w/pi,db);title(幅度響應(yīng)(db);subplot(2,2,3);plot(w/pi,pha/pi);title(相位響應(yīng));axis(0,1,-1,1);subplot(2,2,4);plot(t,ha);title(單位脈沖響應(yīng)ha(t);axis(0,max(t),min(ha),max(ha);% 切比雪夫i型模擬濾波器的設(shè)計(jì)主程序functionb,a=afd_chb1(wp,wr,rp,rs);ep=sqrt(10(rp/10)-1);a=10(rs/20);omegac=rp;omegar=wr/wp;g=sqr

9、t(a*a-1)/ep;n=ceil(log10(g+sqrt(g*g-1)/log10(omegar+sqrt(omegar*omegar-1);fprintf(n=%2.0fn,n);b,a=u_chblap(n,rp,omegac);% 設(shè)計(jì)歸一化切比雪夫i型模擬濾波器原型子程序functionb,a=u_chblap(n,rp,omegac);z,p,k=cheb1ap(n,rp);a=real(poly(p);ann=a(n+1);p=p*omegac;a=real(poly(p);anu=a(n+1);k=k*anu/ann;b0=k;b=real(poly(z);b=k*b;%系

10、統(tǒng)函數(shù)的幅度響應(yīng)和相位響應(yīng)子程序function db,mag,pha,w=freqs_m(b,a,wmax);w=0:1:500*wmax/500;h=freqs(b,a,w);mag=abs(h);db=20*log10(mag+eps)/max(mag);pha=angle(h);%直接形式轉(zhuǎn)化成級(jí)聯(lián)形式子程序 function c,b,a=sdir2cas(b,a);na=length(a)-1;nb=length(b)-1;b0=b(1);b=b/b0;a0=a(1);a=a/a0;c=b0/a0;p=cplxpair(roots(a);k=floor(na/2);if k*2=na

11、 a=zeros(k,3); for n=1:2:na arow=p(n:1:n+1,:);arow=poly(arow); a(fix(n+1)/2),:)=real(arow); endelseif na=1 a=0 real(poly(p);else a=zeros(k+1,3); for n=1:2:2*k arow=p(n:1:n+1,:);arow=ploy(arow); a(fix(n+1)/2),:)=real(arow); end a(k+1,:)=0 real(poly(p(na);endz=cplxpair(roots(b);k=floor(nb/2);if nb=0 b

12、=0 0 poly(z);elseif k*2=nb b=zeros(k,3); for n=1:2:nb brow=z(n:1:n+1,:);brow=poly(brow); b(fix(n+1)/2),:)=real(brow); endelseif nb=1 b=0 real(poly(z); for n=1:2:2*k brow=z(n:1:n+1,:);brow=poly(brow); b(fix(n+1)/2),:)=real(brow); end b=0 real(poly(z(nb);end2.5 運(yùn)行結(jié)果:n= 6b = 0.0614a = 1.0000 0.9283 1.9

13、308 1.2021 0.9393 0.3071 0.0689c = 0.0614b = 0 0 1a = 1.0000 0.4641 0.1247 1.0000 0.3398 0.5577 1.0000 0.1244 0.99072.6仿真波形：3 iir（無(wú)限脈沖響應(yīng)）數(shù)字濾波器設(shè)計(jì)3.1 設(shè)計(jì)題目：用脈沖響應(yīng)不變法設(shè)計(jì)低通濾波器 3.2設(shè)計(jì)指標(biāo)：通帶截止頻率wp=0.1*pi，通帶衰減ap=1db，阻帶截止頻率wr=0.6*pi，阻帶衰減ar=15。 3.3設(shè)計(jì)原理：先設(shè)計(jì)與該數(shù)字濾波器的指標(biāo)相同的模擬濾波器，再通過(guò)s平面映射到z平面。一般的設(shè)計(jì)步驟如下：1）確定數(shù)字低通濾波器的技術(shù)指標(biāo)

14、：通帶截止頻率wp，通帶衰減ap，阻帶截止頻率wf，阻帶衰減ar=15db2）將數(shù)字低通濾波器的技術(shù)指標(biāo)換成模擬信號(hào)的技術(shù)指標(biāo)。3）按照模擬低通濾波器的指標(biāo)設(shè)計(jì)數(shù)字低通濾波器的指標(biāo)。4）將模擬低通濾波器的，從s平面映射到z平面，得到低通數(shù)字濾波器的系統(tǒng)函數(shù)h(z)。3.4 設(shè)計(jì)源程序：%主程序wp=0.1*pi;wr=0.6*pi;ap=1;ar=15;t=1;omegap=wp/t;omegar=wr/t;cs,ds=afd_chb1(omegap,omegar,ap,ar)c,b,a=sdir2cas(cs,ds);db,mag,pha,omega=freqs_m(cs,ds,pi);su

15、bplot(3,2,3);plot(omega/pi,mag);title(模擬濾波器幅度響應(yīng)|ha(jomega)|);b,a=imp_invr(cs,ds,t);h,n=impz(b,a);c,b,a=dir2par(b,a)db,mag,pha,grd,w=freqz_m(b,a);subplot(3,2,1);plot(w/pi,mag);title(數(shù)字濾波器幅度響應(yīng)|ha(jomega)|);subplot(3,2,4);plot(w/pi,db);title(數(shù)字濾波器幅度響應(yīng)（db）);subplot(3,2,2);plot(w/pi,pha/pi);title(數(shù)字濾波器相位

16、響應(yīng));%脈沖響應(yīng)不變法程序function b,a=imp_invr(c,d,t)r,p,k=residuez(c,d);p=exp(p*t);b,a=residuez(r,p,k);b=real(b).*t;a=real(a);%數(shù)字濾波器響應(yīng)子程序function db,mag,pha,grd,w=freqz_m(b,a);h,w=freqz(b,a,1000,whole);h=(h(1:501);w=(w(1:501);mag=abs(h);db=20*log10(mag+eps)/max(mag);pha=angle(h);grd=grpdelay(b,a,w);%直接型轉(zhuǎn)換成并聯(lián)型

17、子程序function c,b,a=dir2par(b,a);m=length(b);n=length(a);r1,p1,c=residuez(b,a);p=cplxpair(p1,10000000*eps);i=cplxcomp(p1,p);r=r1(i);k=floor(n/2);b=zeros(k,2);a=zeros(k,3);if k*2=n; for i=1:2:(n-2) brow=r(i:1:(i+1),:); arow=p(i:1:(i+1),:); brow,arow=residuez(brow,arow,); b(fix(i+1)/2),:)=real(brow); a(

18、fix(i+1)/2),:)=real(arow); end brow,arow=residuez(r(n-1),p(n-1),); b(k,:)=real(brow) 0;a(k,:)=real(arow) 0;else for i=1:2:(n-1) brow=r(i:1:(i+1),:); arow=p(i:1:(i+1),:); brow,arow=residuez(brow,arow,); b(fix(i+1)/2),:)=real(brow); a(fix(i+1)/2),:)=real(arow); endend%比較兩個(gè)含同樣標(biāo)量元素但（可能）有不同下標(biāo)的復(fù)數(shù)對(duì)及其相應(yīng)留數(shù)向量

19、子程序function i=cplxcomp(p1,p2);i=;for j=1:length(p2) for i=1:length(p1) if(abs(p1(i)-p2(j)0.0001) i=i,i; end endend3.5 運(yùn)行結(jié)果：n= 2cs = 0.0970ds =1.0000 0.3449 0.1088c =b = 0.0970 -0.0923a = 1.0000 -1.6171 0.70833.6 仿真波形：4利用fft對(duì)一信號(hào)進(jìn)行頻譜分析 4.1設(shè)計(jì)指標(biāo)：x(n)=0.1*sin(2*pi*f1*n/fs)+0.5*cos(2*pi*f2*n/fs)+0.5*sin(2

20、*pi*f3*n/fs)中含有三種頻率成分，f1=2hz，f2=2.05hz，f3=1.9hz，采樣頻率為fs=10hz，序列：分別取n1=64,n2=128點(diǎn)有效數(shù)據(jù)作頻譜特性分析，分別在四個(gè)圖形窗口繪出x(n)， x(k)64點(diǎn)dft，x(k)補(bǔ)零到128點(diǎn)dft，x(k)128點(diǎn)dft。比較得出在哪種情況下可以清楚地分辨出信號(hào)的頻譜成分。4.2 設(shè)計(jì)源程序如下： f1=2;f2=2.05;f3=1.9;fs=10;n=64;n=0:n-1; xn=0.1*sin(2*pi*f1*n/fs)+0.5*cos(2*pi*f2*n/fs)+0.5*sin(2*pi*f3*n/fs);k=0:n

21、-1;wn=exp(-j*2*pi/n); %旋轉(zhuǎn)因子nk=n*k; %產(chǎn)生一個(gè)臺(tái)nk值的n乘n維矩陣wnnk=wn.nk; %dft矩陣xk=xn*wnnk; %dft系數(shù)的行subplot(2,2,1); %將多個(gè)圖形畫(huà)到一個(gè)平面plot(n,xn); % 繪制函數(shù)曲線 title(a)輸入序列);subplot(2,2,2);plot(k,abs(xk); %繪制幅頻特性圖title(b)x(k)64點(diǎn)dft);n=128;n=0:n-1;xn=0.1*sin(2*pi*f1*n/fs)+0.5*cos(2*pi*f2*n/fs)+0.5*sin(2*pi*f3*n/fs);k=0:n-1;wn=exp(-j*2*pi/n);nk=n*k;wnnk=wn.nk;xk=xn*wnnk;subplot(2,2,3);plot(1:128,abs(xk);title(c)x(k)補(bǔ)零到128點(diǎn)dft);subplot(2,2,4);plot(k,abs(xk);title(d)x(k)128點(diǎn)dft);4.3運(yùn)行結(jié)果：4.4仿真