octave_packages/signal-1.1.3/ncauer.m

   1 ## Copyright (C) 2001 Paulo Neis <p_neis@yahoo.com.br>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## usage: [Zz, Zp, Zg] = ncauer(Rp, Rs, n)
  17 ##
  18 ## Analog prototype for Cauer filter.
  19 ## [z, p, g]=ncauer(Rp, Rs, ws)
  20 ## Rp = Passband ripple
  21 ## Rs = Stopband ripple
  22 ## Ws = Desired order
  23 ##
  24 ## References:
  25 ##
  26 ## - Serra, Celso Penteado, Teoria e Projeto de Filtros, Campinas: CARTGRAF,
  27 ##   1983.
  28 ## - Lamar, Marcus Vinicius, Notas de aula da disciplina TE 456 - Circuitos
  29 ##   Analogicos II, UFPR, 2001/2002.
  30
  31 function [zer, pol, T0]=ncauer(Rp, Rs, n)
  32
  33   ## Cutoff frequency = 1:
  34   wp=1;
  35
  36   ## Stop band edge ws:
  37   ws=__ellip_ws(n, Rp, Rs);
  38
  39   k=wp/ws;
  40   k1=sqrt(1-k^2);
  41   q0=(1/2)*((1-sqrt(k1))/(1+sqrt(k1)));
  42   q= q0 + 2*q0^5 + 15*q0^9 + 150*q0^13; %(....)
  43   D=(10^(0.1*Rs)-1)/(10^(0.1*Rp)-1);
  44
  45   ##Filter order maybe this, but not used now:
  46   ##n=ceil(log10(16*D)/log10(1/q))
  47
  48   l=(1/(2*n))*log((10^(0.05*Rp)+1)/(10^(0.05*Rp)-1));
  49   sig01=0; sig02=0;
  50   for m=0 : 30
  51     sig01=sig01+(-1)^m * q^(m*(m+1)) * sinh((2*m+1)*l);
  52   end
  53   for m=1 : 30
  54     sig02=sig02+(-1)^m * q^(m^2) * cosh(2*m*l);
  55   end
  56   sig0=abs((2*q^(1/4)*sig01)/(1+2*sig02));
  57
  58   w=sqrt((1+k*sig0^2)*(1+sig0^2/k));
  59   #
  60   if rem(n,2)
  61     r=(n-1)/2;
  62   else
  63     r=n/2;
  64   end
  65   #
  66   wi=zeros(1,r);
  67   for ii=1 : r
  68     if rem(n,2)
  69       mu=ii;
  70     else
  71       mu=ii-1/2;
  72     end
  73     soma1=0;
  74     for m=0 : 30
  75       soma1 = soma1 + 2*q^(1/4) * ((-1)^m * q^(m*(m+1)) * sin(((2*m+1)*pi*mu)/n));
  76     end
  77     soma2=0;
  78     for m=1 : 30
  79       soma2 = soma2 + 2*((-1)^m * q^(m^2) * cos((2*m*pi*mu)/n));
  80     end
  81     wi(ii)=(soma1/(1+soma2));
  82   end
  83   #
  84   Vi=sqrt((1-(k.*(wi.^2))).*(1-(wi.^2)/k));
  85   A0i=1./(wi.^2);
  86   sqrA0i=1./(wi);
  87   B0i=((sig0.*Vi).^2 + (w.*wi).^2)./((1+sig0^2.*wi.^2).^2);
  88   B1i=(2 * sig0.*Vi)./(1 + sig0^2 * wi.^2);
  89
  90   ##Gain T0:
  91   if rem(n,2)
  92     T0=sig0*prod(B0i./A0i)*sqrt(ws);
  93   else
  94     T0=10^(-0.05*Rp)*prod(B0i./A0i);
  95   end
  96
  97   ##zeros:
  98   zer=[i*sqrA0i, -i*sqrA0i];
  99
 100   ##poles:
 101   pol=[(-2*sig0*Vi+2*i*wi.*w)./(2*(1+sig0^2*wi.^2)), (-2*sig0*Vi-2*i*wi.*w)./(2*(1+sig0^2*wi.^2))];
 102
 103   ##If n odd, there is a real pole  -sig0:
 104   if rem(n,2)
 105           pol=[pol, -sig0];
 106   end
 107
 108   ##
 109   pol=(sqrt(ws)).*pol;
 110   zer=(sqrt(ws)).*zer;
 111
 112 endfunction
 113
 114 ## usage: ws = __ellip_ws(n, rp, rs)
 115 ## Calculate the stop band edge for the Cauer filter.
 116 function ws=__ellip_ws(n, rp, rs)
 117   kl0 = ((10^(0.1*rp)-1)/(10^(0.1*rs)-1));
 118   k0  = (1-kl0);
 119   int = ellipke([kl0 ; k0]);
 120   ql0 = int(1);
 121   q0  = int(2);
 122   x   = n*ql0/q0;
 123   kl  = fminbnd(@(y) __ellip_ws_min(y,x) ,eps, 1-eps);
 124   ws  = sqrt(1/kl);
 125 endfunction
 126
 127 ## usage: err = __ellip_ws_min(kl, x)
 128 function err=__ellip_ws_min(kl, x)
 129   int=ellipke([kl; 1-kl]);
 130   ql=int(1);
 131   q=int(2);
 132   err=abs((ql/q)-x);
 133 endfunction