/usr/share/octave/packages/signal-1.3.2/ncauer.m

## Copyright (C) 2001 Paulo Neis <p_neis@yahoo.com.br>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.

## usage: [Zz, Zp, Zg] = ncauer(Rp, Rs, n)
##
## Analog prototype for Cauer filter.
## [z, p, g]=ncauer(Rp, Rs, ws)
## Rp = Passband ripple
## Rs = Stopband ripple
## Ws = Desired order
##
## References:
##
## - Serra, Celso Penteado, Teoria e Projeto de Filtros, Campinas: CARTGRAF,
##   1983.
## - Lamar, Marcus Vinicius, Notas de aula da disciplina TE 456 - Circuitos
##   Analogicos II, UFPR, 2001/2002.

function [zer, pol, T0]=ncauer(Rp, Rs, n)

  ## Cutoff frequency = 1:
  wp=1;

  ## Stop band edge ws:
  ws=__ellip_ws(n, Rp, Rs);

  k=wp/ws;
  k1=sqrt(1-k^2);
  q0=(1/2)*((1-sqrt(k1))/(1+sqrt(k1)));
  q= q0 + 2*q0^5 + 15*q0^9 + 150*q0^13; #(....)
  D=(10^(0.1*Rs)-1)/(10^(0.1*Rp)-1);

  ## Filter order maybe this, but not used now:
  ## n=ceil(log10(16*D)/log10(1/q))

  l=(1/(2*n))*log((10^(0.05*Rp)+1)/(10^(0.05*Rp)-1));
  sig01=0; sig02=0;
  for m=0 : 30
    sig01=sig01+(-1)^m * q^(m*(m+1)) * sinh((2*m+1)*l);
  endfor
  for m=1 : 30
    sig02=sig02+(-1)^m * q^(m^2) * cosh(2*m*l);
  endfor
  sig0=abs((2*q^(1/4)*sig01)/(1+2*sig02));

  w=sqrt((1+k*sig0^2)*(1+sig0^2/k));
  ##
  if rem(n,2)
    r=(n-1)/2;
  else
    r=n/2;
  endif
  ##
  wi=zeros(1,r);
  for ii=1 : r
    if rem(n,2)
      mu=ii;
    else
      mu=ii-1/2;
    endif
    soma1=0;
    for m=0 : 30
      soma1 = soma1 + 2*q^(1/4) * ((-1)^m * q^(m*(m+1)) * sin(((2*m+1)*pi*mu)/n));
    endfor
    soma2=0;
    for m=1 : 30
      soma2 = soma2 + 2*((-1)^m * q^(m^2) * cos((2*m*pi*mu)/n));
    endfor
    wi(ii)=(soma1/(1+soma2));
  endfor
  ##
  Vi=sqrt((1-(k.*(wi.^2))).*(1-(wi.^2)/k));
  A0i=1./(wi.^2);
  sqrA0i=1./(wi);
  B0i=((sig0.*Vi).^2 + (w.*wi).^2)./((1+sig0^2.*wi.^2).^2);
  B1i=(2 * sig0.*Vi)./(1 + sig0^2 * wi.^2);

  ## Gain T0:
  if rem(n,2)
    T0=sig0*prod(B0i./A0i)*sqrt(ws);
  else
    T0=10^(-0.05*Rp)*prod(B0i./A0i);
  endif

  ## zeros:
  zer=[i*sqrA0i, -i*sqrA0i];

  ## poles:
  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))];

  ## If n odd, there is a real pole  -sig0:
  if rem(n,2)
    pol=[pol, -sig0];
  endif

  ##
  pol=(sqrt(ws)).*pol;
  zer=(sqrt(ws)).*zer;

endfunction

## usage: ws = __ellip_ws(n, rp, rs)
## Calculate the stop band edge for the Cauer filter.

function ws=__ellip_ws(n, rp, rs)

  kl0 = ((10^(0.1*rp)-1)/(10^(0.1*rs)-1));
  k0  = (1-kl0);
  int = ellipke([kl0 ; k0]);
  ql0 = int(1);
  q0  = int(2);
  x   = n*ql0/q0;
  kl  = fminbnd(@(y) __ellip_ws_min(y,x) ,eps, 1-eps);
  ws  = sqrt(1/kl);

endfunction

## usage: err = __ellip_ws_min(kl, x)

function err=__ellip_ws_min(kl, x)

  int=ellipke([kl; 1-kl]);
  ql=int(1);
  q=int(2);
  err=abs((ql/q)-x);

endfunction
octave-signal 1.3.2-1 / usr / share / octave / packages / signal-1.3.2 / ncauer.m
