摘要
經(jīng)典濾波器的濾波思路是從頻率域上將噪聲濾掉,關(guān)鍵是設(shè)計相應(yīng)的濾波器傳遞函數(shù)H(s)、H(z),分別對應(yīng)著模擬濾波器和數(shù)字濾波器的實現(xiàn)。模擬濾波器主要是通過電感(L)、電容(C)、電阻(R)和運放(OPA)等元器件搭建傳遞函數(shù)為H(s)或者近似為H(s)的硬件電路來實現(xiàn),比如RC濾波電路和有源濾波器等。數(shù)字濾波器(DF)從實現(xiàn)的結(jié)構(gòu)上或者是單位脈沖響h(n)上可以分為無限長脈沖響應(yīng)(IIR)和有限長脈沖響應(yīng)(FIR)濾波器。兩者在結(jié)構(gòu)上的區(qū)別是:IIR有反饋回路,即當(dāng)前輸出y(n)中包含以前輸出y(n-k)(k>0);FIR則沒有反饋回路,當(dāng)前輸出y(n)中只包含輸入x(n)和以前的輸入x(n-k)(k>0)。正是有了反饋回路,導(dǎo)致了IIR單位脈沖響應(yīng)h(n)的無限長。
FIR濾波器
設(shè)計流程如下圖所示(以窗函數(shù)法為例,其它方法不涉及)
IIR濾波器
設(shè)計流程如下圖所示(以巴特沃斯低通濾波器采用雙線性變換法為例)
數(shù)字濾波器的實現(xiàn)
由上面分析,數(shù)字濾波器的兩種形式在得到的最后表達(dá)式上都可以歸一化為N階線性常系數(shù)差分方程的一般形式:
y(n)=∑Mk=0akx(n?k)?∑Nr=1bry(n?r)
這樣一個系統(tǒng)的實現(xiàn)可以分別從硬件和軟件上實現(xiàn):硬件上只需要包括移位寄存器、乘法器和加法器就可以;軟件上實現(xiàn)就較為簡單了,直接高級語言描述。
兩種經(jīng)典數(shù)字濾波器的比較
如下圖所示:
基于Matlab的實驗
用窗函數(shù)法設(shè)計FIR數(shù)字濾波器
設(shè)計一線性相位低通FIR數(shù)字濾波器,通帶截止頻率wp=0.2π,通帶最大衰減0.25dB,阻帶起始頻率ws=0.3π,阻帶最小衰減50dB。請給出N和所選窗函數(shù),求出h(n),并繪制相應(yīng)的幅頻特性(dB)曲線,根據(jù)該圖給出實際的阻帶衰減,繪制窗函數(shù)的時域頻域特性曲線。
實驗過程分析:
阻帶最小衰減為50dB,對于各種窗函數(shù)的基本參數(shù),不能用矩形窗和漢寧窗函數(shù),因為它們的阻帶最小衰減不能滿足設(shè)計的要求,因此可以選用漢明和布萊克曼窗進(jìn)行FIR數(shù)字濾波器的設(shè)計。為了對比漢寧窗不能滿足設(shè)計要求,同樣做出漢寧窗設(shè)計的FIR數(shù)字濾波器。
根據(jù)通帶截止頻率和阻帶起始頻率得到過渡帶的帶寬,漢明窗對應(yīng)的數(shù)字濾波器的過渡帶寬為8π÷N,推出用漢明窗設(shè)計時的數(shù)字濾波器的至少N=80,為了取用類型I,即偶對稱、N為奇數(shù),取用N=81;對于用布萊克曼窗對應(yīng)的數(shù)字濾波器的過渡帶寬為12π÷N,推出N=120,取N=121。漢寧窗與漢明窗的N值計算相同。
首先給出漢明窗的FIR實驗結(jié)果:
察漢明窗設(shè)計的FIR濾波器的幅頻特性曲線,在通帶內(nèi),最大衰減小于0.25db,在阻帶內(nèi),最小衰減大于50db,滿足了設(shè)計要求。
然后給出布萊克曼窗的FIR實驗結(jié)果:
觀察布萊克曼窗設(shè)計的FIR濾波器的幅頻特性曲線,在通帶內(nèi),最大衰減小于0.25db,在阻帶內(nèi),最小衰減大于50db,滿足了設(shè)計要求。
最后給出漢寧窗的FIR實驗結(jié)果作為對比:
觀察漢寧窗設(shè)計的FIR濾波器的幅頻特性曲線,在通帶內(nèi),最大衰減小于0.25db,但在阻帶內(nèi)最小衰減小于了50db,所以不滿足設(shè)計要求。
一般地,隨著N值增大,過渡帶愈加陡峭,衰減更快。旁瓣越多對應(yīng)的通帶起伏越大。為了獲得更加平穩(wěn)的通帶和更加陡峭的過渡帶,雖然不能兼得,一般處理時是通過增加主瓣寬度來換取對旁瓣的抑制。
用雙線性變換法設(shè)計IIR數(shù)字濾波器
用雙線性變換法設(shè)計一個IIR巴特沃斯低通數(shù)字濾波器。設(shè)計指標(biāo)參數(shù)為:在通帶內(nèi)頻率低于0.2π時,最大衰減小于1dB;在阻帶內(nèi)[0.3π,π]頻率區(qū)間上,最小衰減大于15dB。用所設(shè)計的濾波器對實際心電圖信號采樣序列進(jìn)行仿真濾波處理,并分別繪出濾波前后的心電圖信號波形圖。觀察總結(jié)濾波作用與效果。
實驗結(jié)果:
實驗結(jié)果分析:
由幅頻特性曲線可以看出設(shè)計的IIR濾波器滿足要求,同時經(jīng)過心電信號的測試,能夠達(dá)到濾除高頻噪聲的目的。
本實驗中,經(jīng)IIR設(shè)計的濾波器傳遞函數(shù)H(z)為:
計算得到一般的線性系統(tǒng)表達(dá)式:
y(n)=0.0007378x(n)+0.004427x(n?1)+0.01107x(n?2)+0.01476x(n?3)+0.01107x(n?4)+0.004427x(n?5)+0.0007378x(n?6)+3.184y(n?1)?4.622y(n?2)+3.779y(n?3)?1.814y(n?4)+0.48y(n?5)?0.05445y(n?6)
一般設(shè)計濾波器就是為了得到上面的表達(dá)式,有了這樣的表達(dá)式之后就可以很方便的通過硬件或者軟件實現(xiàn)濾波器了。
Matlab代碼
FIR濾波器
理想低通濾波器的構(gòu)建:
function hd= ideal_lp( wc,N ) %hd=點0到N-1之間的理想脈沖響應(yīng) % wc=截止頻率(弧度) %N=濾波器的長度 tao=(N-1)/2; n=[0:(N-1)]; %+eps轉(zhuǎn)換成浮點數(shù) m=n-tao+eps; hd=sin(wc*m)./(pi*m); end
布萊克曼窗FIR濾波器:
%阻帶最小衰減為50DB,所以可以選擇漢明窗 wp=0.2*pi;ws=0.3*pi; deltaw=ws-wp;%過渡帶寬 N0=12*pi/deltaw; %N0=ceil(6.6*pi/deltaw); N=N0+mod(N0+1,2);%為實現(xiàn)FIR類型I偶對稱濾波器,確保N為奇數(shù)。 windows=(blackman(N))'; wc=(ws+wp)/2; hd=ideal_lp(wc,N); b=hd.*windows; [H,w]=freqz(b,1); n=0:N-1; dw=2*pi/100; dbH=20*log10((abs(H)+eps)/max(abs(H)));%化為分貝值 %Rp=-(min(db(1:wp/dw+1)))%檢驗通帶波動 %As=-round(max(db(ws/dw+1:501)))%檢驗最小阻帶衰減 %作圖 subplot(221) stem(n,hd); axis([0,N,1.1*min(hd),1.1*max(hd)]);title('理想脈沖響應(yīng)'); xlabel('n'); ylabel('hd(n)'); subplot(222); stem(n,windows); axis([0,N,0,1.1]); title('布萊克曼窗函數(shù)特性'); xlabel('n'); ylabel('wd(n)'); subplot(223),stem(n,b); axis([0,N,1.1*min(b),1.1*max(b)]);title('實際脈沖響應(yīng)'); xlabel('n'); ylabel('h(n)'); subplot(224); plot(w/pi,dbH); title('幅度頻率響應(yīng)'); axis([0,1,-200,10]); xlabel('頻率(單位:pi)'); ylabel('H(e^{jomega})'); set(gca,'XTickMode','manual','XTick',[0,wp/pi,ws/pi,1]); set(gca,'YTickMode','manual','YTick',[-200,-50,-0.25,0]); grid
漢明窗FIR濾波器:
%阻帶最小衰減為50DB,所以可以選擇漢明窗 wp=0.2*pi;ws=0.3*pi; deltaw=ws-wp;%過渡帶寬 N0=8*pi/deltaw; %N0=ceil(6.6*pi/deltaw); N=N0+mod(N0+1,2);%為實現(xiàn)FIR類型I偶對稱濾波器,確保N為奇數(shù)。 windows=(hamming(N))'; %windows=ones(1,N); wc=(ws+wp)/2; hd=ideal_lp(wc,N); b=hd.*windows; [H,w]=freqz(b,1); n=0:N-1; dw=2*pi/100; dbH=20*log10((abs(H)+eps)/max(abs(H)));%化為分貝值 %Rp=-(min(db(1:wp/dw+1)))%檢驗通帶波動 %As=-round(max(db(ws/dw+1:501)))%檢驗最小阻帶衰減 %作圖 figure(1) %subplot(221) stem(n,hd); axis([0,N,1.1*min(hd),1.1*max(hd)]);title('理想脈沖響應(yīng)'); xlabel('n'); ylabel('hd(n)'); figure(2) %subplot(222); stem(n,windows); axis([0,N,0,1.1]); title('漢明窗函數(shù)特性'); xlabel('n'); ylabel('wd(n)'); figure(3) %subplot(223), stem(n,b); axis([0,N,1.1*min(b),1.1*max(b)]);title('實際脈沖響應(yīng)'); xlabel('n'); ylabel('h(n)'); figure(4) %subplot(224); plot(w/pi,dbH); title('幅度頻率響應(yīng)'); axis([0,1,-80,10]); xlabel('頻率(單位:pi)'); ylabel('H(e^{jomega})'); set(gca,'XTickMode','manual','XTick',[0,wp/pi,ws/pi,1]); set(gca,'YTickMode','manual','YTick',[-80,-50,-0.25,0]); grid
漢寧窗FIR濾波器:
%阻帶最小衰減為50DB,所以可以選擇漢明窗 wp=0.2*pi;ws=0.3*pi; deltaw=ws-wp;%過渡帶寬 N0=8*pi/deltaw; %N0=ceil(6.6*pi/deltaw); N=N0+mod(N0+1,2);%為實現(xiàn)FIR類型I偶對稱濾波器,確保N為奇數(shù)。 windows=(hanning(N))'; wc=(ws+wp)/2; %計算得到理想低通濾波器,已經(jīng)移位(N-1)/2 hd=ideal_lp(wc,N); %加窗施加在hn上 b=hd.*windows; % [H,w]=freqz(b,1); n=0:N-1; dw=2*pi/100; dbH=20*log10((abs(H)+eps)/max(abs(H)));%化為分貝值 %Rp=-(min(db(1:wp/dw+1)))%檢驗通帶波動 %As=-round(max(db(ws/dw+1:501)))%檢驗最小阻帶衰減 %作圖 subplot(221) stem(n,hd); axis([0,N,1.1*min(hd),1.1*max(hd)]);title('理想脈沖響應(yīng)'); xlabel('n'); ylabel('hd(n)'); subplot(222); stem(n,windows); axis([0,N,0,1.1]); title('漢寧窗函數(shù)特性'); xlabel('n'); ylabel('wd(n)'); subplot(223),stem(n,b); axis([0,N,1.1*min(b),1.1*max(b)]);title('實際脈沖響應(yīng)'); xlabel('n'); ylabel('h(n)'); subplot(224); plot(w/pi,dbH); title('幅度頻率響應(yīng)'); axis([0,1,-100,10]); xlabel('頻率(單位:pi)'); ylabel('H(e^{jomega})'); set(gca,'XTickMode','manual','XTick',[0,wp/pi,ws/pi,1]); set(gca,'YTickMode','manual','YTick',[-100,-50,-0.25,0]); grid
IIR數(shù)字濾波器
IIR濾波器設(shè)計:
%雙線性變換前預(yù)畸 Fs=500; wp=(100/Fs)*2*pi; ws=(200/Fs)*2*pi; Rp=2; Rs=15; wp2=2*Fs*tan(wp/2); ws2=2*Fs*tan(ws/2); %選擇濾波器的最小階數(shù) [N,wn]=buttord(wp2,ws2,Rp,Rs,'s');%注意此處輸入的是畸變后的指標(biāo),輸出N為符合要求的模擬濾波器的最小階數(shù),wn為3dB帶寬 %創(chuàng)建butterworth模擬濾波器 [Z,P,K]=buttap(N); %把濾波器零極點模型轉(zhuǎn)化為傳遞函數(shù)模型 [Bap,Aap]=zp2tf(Z,P,K); %把模擬濾波器原型轉(zhuǎn)換為截止頻率為wn的模擬低通濾波器 [b,a]=lp2lp(Bap,Aap,wn); %用雙線性法實現(xiàn)模擬濾波器到數(shù)字濾波器的轉(zhuǎn)換 [bz,az]=bilinear(b,a,Fs); %繪制頻率響應(yīng)曲線 [H,W]=freqz(bz,az); subplot(2,1,1); plot(W/pi,abs(H)); grid; xlabel('頻率w/pi'); ylabel('幅度絕對值'); subplot(2,1,2); plot(W/pi,20*log10((abs(H)+eps)/max(abs(H)))); grid; xlabel('頻率w/pi'); ylabel('幅度dB');
心電信號濾波:
%數(shù)字濾波器指標(biāo): wc=0.2*pi;ws=0.3*pi;Rp=1;As=15; ripple=10^(-Rp/20);Attn=10^(-As/20); %轉(zhuǎn)換成為模擬指標(biāo): Fs=1;T=1/Fs; Omgp=(2/T)*tan(wc/2); Omgs=(2/T)*tan(ws/2); %模擬原型濾波器計算: [n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s');%計算階數(shù)和截止頻率 [z0,p0,k0]=buttap(n); ba=k0*real(poly(z0)); aa=real(poly(p0)); [ba1,aa1]=lp2lp(ba,aa,Omgc); %雙線性變換法計算數(shù)字濾波器系數(shù) [bd,ad]=bilinear(ba1,aa1,Fs);%雙線性變換法求數(shù)字濾波器系數(shù)b,a tf(bd,ad,1) [sos,g]=tf2sos(bd,ad)%由直接型轉(zhuǎn)變成級聯(lián)型 %求數(shù)字系統(tǒng)的頻率特性: [H,w]=freqz(bd,ad); dbH=20*log10((abs(H)+eps)/max(abs(H)));%化為分貝值 subplot(2,1,1); plot(w/pi,abs(H)); % 加上對應(yīng)取值軸來驗證是否滿足設(shè)計要求 set(gca,'XTickMode','Manual','XTick',[0,0.2,0.3,pi]) set(gca,'YTickMode','Manual','YTick',[0,Attn,ripple,1]) ylabel('|H|'); xlabel('數(shù)字頻率(*pi)') title('幅度響應(yīng)'); axis([0,1/2,0,1.1]); grid on subplot(212); plot(w/pi,dbH); xlabel('數(shù)字頻率(*pi)'); ylabel('幅度的分貝值DB'); % 加上對應(yīng)取值軸來驗證是否滿足設(shè)計要求 set(gca,'XTickMode','Manual','XTick',[0,0.2,0.3,pi]) set(gca,'YTickMode','Manual','YTick',[-50,-As,-Rp,0]) axis([0,1/2,-50,0.1]); grid on %心電信號: xn=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,... 0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,... 6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0] figure; subplot(211); stem(xn); title('原始心電信號') yn=filter(bd,ad,xn); subplot(212); stem(yn); title('經(jīng)過低通濾波后的心電信號') figure; subplot(211); stem(abs(fft(xn))); title('原始心電信號的fft頻譜') yn=filter(bd,ad,xn); subplot(212); stem(abs(fft(yn)));
title('經(jīng)過低通濾波后的心電信號的fft頻譜');
-
濾波器
+關(guān)注
關(guān)注
161文章
7816瀏覽量
178110 -
fir濾波器
+關(guān)注
關(guān)注
1文章
95瀏覽量
19041 -
IIR濾波器
+關(guān)注
關(guān)注
0文章
31瀏覽量
11512
原文標(biāo)題:經(jīng)典濾波器設(shè)計
文章出處:【微信號:zhuyandz,微信公眾號:FPGA之家】歡迎添加關(guān)注!文章轉(zhuǎn)載請注明出處。
發(fā)布評論請先 登錄
相關(guān)推薦
評論