%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);

%w=-pi:0.1:pi;
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