matlab的數字濾波器的仿真怎麽實現

fs = 35;

Fs = 100;

wp = 2*pi*fp/Fs;

ws = 2*pi*fs/Fs;

wp = tan(wp/2);

ws = tan(ws/2); % 通帶最大衰減為0.5dB，阻帶最小衰減為40dB

[N, wn] = buttord(wp, ws, 0.5, 40, 's'); % 模擬低通濾波器極零點

[z, p, k] = buttap(N); % 由極零點獲得轉移函數參數

[b, a] = zp2tf(z, p, k); % 由原型濾波器獲得實際低通濾波器

[B, A] = lp2lp(b, a, wp);

[bz, az] = bilinear(B, A, .5);

[h, w] = freqz(bz, az, 256, Fs);

figure

plot(w, abs(h))

grid on

圖1 巴特沃斯數字低通濾波器

1-2基於Butterworth模擬濾波器原型，使用雙線性狀換設計數字濾波器：各參數值為：通帶截止頻率Omega=0.2*pi,阻帶截止頻率Omega=0.3*pi,通帶波動值Rp=1dB,阻帶波動值Rs=15dB,設Fs=4000Hz。

代碼：

wp=0.2*pi;ws=0.3*pi;

Fs=4000;T=1/Fs;

OmegaP=(2/T)*tan(wp/2);

OmegaS=(2/T)*tan(ws/2);

rp=1;rs=15;as=15;

ripple=10^(-rp/20);attn=10^(-rs/20);

[n,wn]=buttord(OmegaP,OmegaS,rp,rs,'s');

[z,p,k]=Buttap(n);

[b,a]=zp2tf(z,p,k);

[bt,at]=lp2lp(b,a,wn);

[b,a]=bilinear(bt,at,Fs);

[db,mag,pha,grd,w]=freqz_m(b,a);

%

%下面繪出各條曲線

subplot(2,2,1);plot(w/pi,mag);title('Magnitude Frequency幅頻特性');

xlabel('w(/pi)');ylabel('|H(jw)|');

axis([0,1,0,1.1]);

set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);

set(gca,'YTickMode','manual','YTick',[0 attn ripple 1]);grid

subplot(2,2,2);plot(w/pi,db);title('Magnitude Frequency幅頻特性(db)');

xlabel('w(/pi)');ylabel('dB');

axis([0,1,-30,5]);

set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);

set(gca,'YTickMode','manual','YTick',[-60 -as -rp 0]);grid

subplot(2,2,3);plot(w/pi,pha/pi);title('Phase Frequency相頻特性');

xlabel('w(/pi)');ylabel('pha(/pi)');

axis([0,1,-1,1]);

subplot(2,2,4);plot(w/pi,grd);title('Group Delay群延時');

xlabel('w(/pi)');ylabel('Sample');

axis([0,1,0,15]);

set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);grid

運行結果：

圖2巴特沃思數字低通濾波器幅頻-相頻特性

1-3設計壹巴特沃斯高通數字濾波器，要求通帶截止頻率0.6*pi，通帶衰減不大於1dB，阻帶衰減15DB，抽樣T=1。

Wp=0.6*pi;

Ws=0.4*pi;

Ap=1;

As=15;

[N,wn]=buttord(Wp/pi,Ws/pi,Ap,As); %計算巴特沃斯濾波器階次和截止頻率

%頻率變換法設計巴特沃斯高通濾波器

[db,mag,pha,grd,w]=freqz_m(b,a); %數字濾波器響應

plot(w,mag);

title('數字濾波器幅頻響應|H(ej\Omega)|')

圖3巴特沃斯數字高通濾波器

2-1用窗函數法設計壹個線性相位FIR低通濾波器，並滿足性能指標:通帶邊界頻率

Wp=0.5*pi,阻帶邊界頻率Ws=0.66*pi,阻帶衰減不小於40dB,通帶波紋不大於3dB。選擇漢寧窗。

代碼：

wp =0.5*pi;

ws=0.66*pi;

wdelta =ws-wp;

N= ceil(8*pi/wdelta)

if rem(N,2)==0

N=N+1;

end

);