sy5.txt

几个Matlab编程示例

%基于Matlab的数字滤波器设计

%1、IIR(Butterworth)
ws=[0.3 0.75];wp=[0.45 0.65];rp=1;rs=40;
[n1,wn1]=buttord(wp,ws,rp,rs);
[num,den]=butter(n1,wn1);
[g,w,pha]=gain(num,den);
subplot(1,2,1);plot(w/pi,g);grid
axis([0 1 -60 5]);
xlabel('\omega/\pi');ylabel('Gain in dB');
title('IIR,gain');
subplot(1,2,2);plot(w/pi,pha);grid
xlabel('\omega/\pi');ylabel('phase');
title('IIR,phase');

%子程序【M文件】
function [g,w,pha]=gain(num,den)
w=0:pi/255:pi;
h=freqz(num,den,w);
g=20*log10(abs(h));
pha=angle(h);



%2、FIR(Blackman)
ws1=0.3*pi;wp1=0.45*pi;
wp2=0.65*pi;ws2=0.75*pi;
As=40;
tr_width=min((wp1-ws1),(ws2-wp2))
M=ceil(11*pi/tr_width)+1
wc1=(ws1+wp1)/2;wc2=(wp2+ws2)/2;
hd=ideal_lp(wc2,M)-ideal_lp(wc1,M);
w_bla=(blackman(M))';
h=hd.*w_bla;
[db,mag,pha,grd,w]=freqz_m(h,[1]);
delta_w=2*pi/1000;
rp=-min(db(wp1/delta_w+1:1:wp2/delta_w))
As=-round(max(db(ws2/delta_w+1:1:501)))
subplot(1,2,2);plot(w/pi,pha);title('FIR,phase')
xlabel('频率（单位：pi）');ylabel('相位');subplot(1,2,1);plot(w/pi,db);
title('幅度响应（单位：dB）');grid;
xlabel('频率（单位：pi）');ylabel('分贝数')
axis([0 1 -150 10]);
set(gca,'XTickMode','manual','XTick',[0,0.3,0.45,0.65,0.75,1])
set(gca,'YTickMode','manual','YTick',[-40,0])

%子程序【M文件】
function hd=ideal_lp(wc,M);
alpha=(M-1)/2;
n=[0:1:(M-1)];
m=n-alpha+eps;
hd=sin(wc*m)./(pi*m);

%子程序【M文件】
function[db,mag,pha,grd,w]=freqz_m(b,a);
[H,w]=freqz(b,a,1000,'whole')
H=(H(1:1:501))';
w=(w(1:1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
grd=grpdelay(b,a,w);