幼之梨
2008-07-11, 23:41
Wp=30;Ws=50;Rp=1;As=30;
Ripple=10^(-Rp/20);
Attn=10^(-As/20);
[b,a]=afd_butt(Wp,Ws,Rp,As)
[C,B,A]=sdir2cas(b,a)
[db,mag,pha,w]=freqs_m(b,a,50);
[ha,x,t]=impulse(b,a);
figure(1);clf;
subplot(2,2,1);plot(w,mag);title('Magnitude Response');
xlabel('Analog frequency in rad/s');
ylabel('H');
%
axis([0,50,0,1.1])
set(gca,'Xtickmode','manual','Xtick',[0,30,40,50]);
set(gca,'Ytickmode','manual','Ytick',[0,Attn,Ripple,1]);
grid
subplot(2,2,2);plot(w,db);title('Magnitude in dB');
xlabel('Analog frequency in rad/s');
ylabel('decibels');
%
axis([0,50,-40,5])
set(gca,'Xtickmode','manual','Xtick',[0,30,40,50]);
set(gca,'Ytickmode','manual','Ytick',[-40,-As,-Rp,0]);
grid
subplot(2,2,3);plot(w,pha/pi);title('Phase Response');
xlabel('Analog frequency in rad/s');
ylabel('radians');
%
axis([0,50,-1.1,1.1])
set(gca,'Xtickmode','manual','Xtick',[0,30,40,50]);
set(gca,'Ytickmode','manual','Ytick',[-1,-0.5,0,0.5,1]);
grid
subplot(2,2,4);plot(t,ha);title('Impulse Response');
xlabel('time in seconds');
ylabel('ha(t)');
axis([0,max(t)+0.05,min(ha),max(ha)+0.025]);
set(gca,'Xtickmode','manual','Xtick',[0,0.1,max(t)]);
set(gca,'Ytickmode','manual','Ytick',[0,0.1,max(ha)]);
grid
%巴特沃兹模拟滤波器的设计程序
function[b,a]=afd_butt(Wp,Ws,Rp,As);
if Wp<=0
error('Passband edge must be larger than 0')
end
if Ws<=Wp
error('Stopband edge must be larger than Passed edge')
end
if(Rp<=0)|(As<0)
error('PB ripple and /Or SB attenuation must be larger than 0')
end
N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(Wp/Ws)));
OmegaC=Wp/((10^(Rp/10)-1)^(1/(2*N)));
[b,a]=u_buttap(N,OmegaC);
%设计非归一化巴特沃兹模拟低通滤波器原型子程序
function [b,a]=u_buttap(N,OmegaC);
[z,p,k]=buttap(N);
p=p*OmegaC;
k=k*OmegaC^N;
B=real(poly(z));
b0=k;
b=k*B;
a=real(poly(z));
%计算系统函数的幅度响应和相位响应子程序
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);
%直接形式转换成级联形式子程序
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
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);
end
elseif 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=poly(Arow);
A(fix((n+1)/2),:)=real(Arow);
end
A(K+1,:)=[0 real(poly(p(Na)))];
end
z=cplxpair(roots(b));K=floor(Nb/2);
if Nb==0
B=[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);
end
elseif Nb==1
B=[0 real(poly(z))];
else
B=zeros(K+1,3);
for n=1:2:2*K
Brow=z(n:1:n+1,:);Brow=poly(Brow);
B(fix((n+1)/2),:)=real(Brow);
end
B(K+1,:)=[0 real(poly(z(Nb)))];
end
这是一个巴特沃兹滤波器的设计,执行后出现??? One or more output arguments not assigned during call to 'sdir2cas'.
请教高手这是什么问题呢?怎么解决?谢谢!!!!
Ripple=10^(-Rp/20);
Attn=10^(-As/20);
[b,a]=afd_butt(Wp,Ws,Rp,As)
[C,B,A]=sdir2cas(b,a)
[db,mag,pha,w]=freqs_m(b,a,50);
[ha,x,t]=impulse(b,a);
figure(1);clf;
subplot(2,2,1);plot(w,mag);title('Magnitude Response');
xlabel('Analog frequency in rad/s');
ylabel('H');
%
axis([0,50,0,1.1])
set(gca,'Xtickmode','manual','Xtick',[0,30,40,50]);
set(gca,'Ytickmode','manual','Ytick',[0,Attn,Ripple,1]);
grid
subplot(2,2,2);plot(w,db);title('Magnitude in dB');
xlabel('Analog frequency in rad/s');
ylabel('decibels');
%
axis([0,50,-40,5])
set(gca,'Xtickmode','manual','Xtick',[0,30,40,50]);
set(gca,'Ytickmode','manual','Ytick',[-40,-As,-Rp,0]);
grid
subplot(2,2,3);plot(w,pha/pi);title('Phase Response');
xlabel('Analog frequency in rad/s');
ylabel('radians');
%
axis([0,50,-1.1,1.1])
set(gca,'Xtickmode','manual','Xtick',[0,30,40,50]);
set(gca,'Ytickmode','manual','Ytick',[-1,-0.5,0,0.5,1]);
grid
subplot(2,2,4);plot(t,ha);title('Impulse Response');
xlabel('time in seconds');
ylabel('ha(t)');
axis([0,max(t)+0.05,min(ha),max(ha)+0.025]);
set(gca,'Xtickmode','manual','Xtick',[0,0.1,max(t)]);
set(gca,'Ytickmode','manual','Ytick',[0,0.1,max(ha)]);
grid
%巴特沃兹模拟滤波器的设计程序
function[b,a]=afd_butt(Wp,Ws,Rp,As);
if Wp<=0
error('Passband edge must be larger than 0')
end
if Ws<=Wp
error('Stopband edge must be larger than Passed edge')
end
if(Rp<=0)|(As<0)
error('PB ripple and /Or SB attenuation must be larger than 0')
end
N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(Wp/Ws)));
OmegaC=Wp/((10^(Rp/10)-1)^(1/(2*N)));
[b,a]=u_buttap(N,OmegaC);
%设计非归一化巴特沃兹模拟低通滤波器原型子程序
function [b,a]=u_buttap(N,OmegaC);
[z,p,k]=buttap(N);
p=p*OmegaC;
k=k*OmegaC^N;
B=real(poly(z));
b0=k;
b=k*B;
a=real(poly(z));
%计算系统函数的幅度响应和相位响应子程序
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);
%直接形式转换成级联形式子程序
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
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);
end
elseif 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=poly(Arow);
A(fix((n+1)/2),:)=real(Arow);
end
A(K+1,:)=[0 real(poly(p(Na)))];
end
z=cplxpair(roots(b));K=floor(Nb/2);
if Nb==0
B=[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);
end
elseif Nb==1
B=[0 real(poly(z))];
else
B=zeros(K+1,3);
for n=1:2:2*K
Brow=z(n:1:n+1,:);Brow=poly(Brow);
B(fix((n+1)/2),:)=real(Brow);
end
B(K+1,:)=[0 real(poly(z(Nb)))];
end
这是一个巴特沃兹滤波器的设计,执行后出现??? One or more output arguments not assigned during call to 'sdir2cas'.
请教高手这是什么问题呢?怎么解决?谢谢!!!!