From octave-sources-request at bevo dot che dot wisc dot edu  Mon Jan 15 08:47:56 2001
Subject: Complement of cheb1ord.m
From: Laurent Mazet <mazet at crm dot mot dot com>
To: octave-sources at bevo dot che dot wisc dot edu
Cc: pkienzle at kienzle dot powernet dot co dot uk
Date: Mon, 15 Jan 2001 15:41:34 +0100 (MET)


Hi,

I've just completed the cheb1ord function from the Paul Kienzle's package. It
is not fully functionnal (can't do reject band fiter design) but it can be
use for other Chebychev filter design.

      Laurent

## Copyright (C) 2000 Paul Kienzle
##
## This program is free software and may be used for any purpose.  This
## copyright notice must be maintained.  Paul Kienzle is not responsible
## for the consequences of using this software.
##
## Completed by: Laurent S. Mazet

## Compute chebyshev type I filter order and cutoff for the desired response
## characteristics. Rp is the allowable decibels of ripple in the pass 
## band. Rs is the minimum attenuation in the stop band.
##
## [n, Wc] = cheb1ord(Wp, Ws, Rp, Rs)
##     Low pass (Wp<Ws) or high pass (Wp>Ws) filter design.  Wp is the
##     pass band edge and Ws is the stop band edge.  Frequencies are
##     normalized to [0,1], corresponding to the range [0,Fs/2].
## 
## [n, Wc] = cheb1ord([Wp1, Wp2], [Ws1, Ws2], Rp, Rs)
##     Band pass (Ws1<Wp1<Wp2<Ws2) or band reject (Wp1<Ws1<Ws2<Wp2)
##     filter design. Wp gives the edges of the pass band, and Ws gives
##     the edges of the stop band.
##
## Theory: 
##
## See also: cheby1

function [n, Wc] = cheb1ord(Wp, Ws, Rp, Rs)

  if nargin != 4
    usage("[n, Wn] = cheb1ord(Wp, Ws, Rp, Rs)");
  end
  if length(Wp) != length(Ws)
    error("cheb1ord: Wp and Ws must have the same length");
  end
  if length(Wp) != 1 && length(Wp) != 2
    error("cheb1ord: Wp,Ws must have length 1 or 2");
  end
  if length(Wp) == 2 && (all(Wp>Ws) || all(Ws>Wp) || diff(Wp)<=0 || diff(Ws)<=0)
    error("cheb1ord: Wp(1)<Ws(1)<Ws(2)<Wp(2) or Ws(1)<Wp(1)<Wp(2)<Ws(2)");
  end

  T = 2;

  ## test for low pass
  low_pass = 0;
  if (Wp(1) < Ws(1))
    low_pass=1;
  endif;
  
  ## warp the target frequencies according to the bilinear transform
  Ws = (2/T)*tan(pi*Ws./T);
  Wp = (2/T)*tan(pi*Wp./T);

  if (low_pass)
    ## low pass
    if (length(Wp) == 1)
      Wa = Ws/Wp;
    else
      ## band reject
      error ("band reject is not implement yet.");
    endif;
  else
   ## if high pass, reverse the sense of the test
   if (length(Wp) == 1)
      Wa = Wp/Ws;
    else
      ## band pass 
      Wa=(Ws.^2 - Wp(1)*Wp(2))./(Ws*(Wp(1)-Wp(2)));
    endif;
  endif;
  Wa = min(abs(Wa));
  
  ## compute minimum n which satisfies all band edge conditions
  ## the factor 1/length(Wp) is an artificial correction for the
  ## band pass/stop case, which otherwise significantly overdesigns.
  stop_atten = 10^(abs(Rs)/10);
  pass_atten = 10^(abs(Rp)/10);

  ## compute -3dB cutoff given Wp, Rp and n
  n = ceil(acosh(sqrt((stop_atten-1)/(pass_atten-1)))/acosh(Wa));
  
  ## unwarp the returned frequency
  Wc = atan(T/2*Wp)*T/pi;
      
endfunction

-- 
Dr. Laurent Mazet  === Research Engineer === Centre de Recherche de MOTOROLA
Tel: +33 (0)1 69 35 48 30 ========================= Email: mazet at crm dot mot dot com



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------