dtransform.m

这是数字信号处理——matlab释义与实现一书的源代码

% 《数字信号处理教程——MATLAB释义与实现》第四章演示程序
%  电子工业出版社出版  陈怀琛编著 2004年9月
%  数字滤波器频带变换演示
% dtranform.m  数字滤波器设计，从低通到其他类型滤波器的频率对应关系

function  [out0,out1]=dtransform(w,string)
%下面的参数为设计滤波器时给定的参数
wu=0.7*pi;
wl=0.5*pi;
wc=0.4*pi; %最终设计的低通滤波器的截止频率
wh=0.7*pi; %最终设计的高通滤波器的下截止频率
wo=0.2*pi; %低通滤波器3dB截止频率

a=cos((wu+wl)/2)/cos((wu-wl)/2);
k=cot((wu-wl)/2)*tan(wo/2);
al=sin((wo-wc)/2)/sin((wo+wc)/2);
ah=-cos((wh+wo)/2)/cos((wh-wo)/2);

zl=exp(j*w);

if (string=='lp2lp')
  %低通到低通频率的变换关系
   for i=1:length(w)
      pl2l(i)=(zl(i)+al)/(1+al*zl(i));
   end
   out1=pl2l;
   out0=pl2l;

elseif (string=='lp2hp')
   %低通到高通频率的变换关系
   for i=1:length(w)
     pl2h(i)=-(zl(i)+ah)/(1+ah*zl(i));
   end
   out1=pl2h;
   out0=pl2h;

elseif (string=='lp2bp')
   %低通到带通频率的变换关系！
   for i=1:length(w)
      b2=zl(i)*(k-1)/(k+1)+1;
      b1=(2*a*k*zl(i)/(k+1)+2*a*k/(k+1));
      b0=(zl(i)+(k-1)/(k+1));
      pole=roots([b2,-b1,b0]);
     if (imag(pole(1))>0)
      pl2bp0(i)=pole(1);
      pl2bp1(i)=pole(2);
     else
     pl2bp0(i)=pole(2);
     pl2bp1(i)=pole(1);
    end
   end
   out0=pl2bp0;
   out1=pl2bp1;
   
else
%低通到带阻频率的变换关系！
  for i=1:length(w)
     c2=zl(i)*(1-k)/(k+1)-1;
     c1=(2*a*k*zl(i)/(k+1)-2*a*k/(k+1));
     c0=(zl(i)-(1-k)/(k+1));
     pole=roots([c2,-c1,c0]);
    if (imag(pole(1))>0)
     pl2bs0(i)=pole(1);
     pl2bs1(i)=pole(2);
   else
     pl2bs0(i)=pole(2);
     pl2bs1(i)=pole(1);
    end
  end
  out0=pl2bs0;
  out1=pl2bs1;
end