octave_packages/signal-1.1.3/cheb1ord.m

   1 ## Copyright (C) 2000 Paul Kienzle <pkienzle@users.sf.net>
   2 ## Copyright (C) 2000 Laurent S. Mazet
   3 ##
   4 ## This program is free software; you can redistribute it and/or modify it under
   5 ## the terms of the GNU General Public License as published by the Free Software
   6 ## Foundation; either version 3 of the License, or (at your option) any later
   7 ## version.
   8 ##
   9 ## This program is distributed in the hope that it will be useful, but WITHOUT
  10 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  11 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  12 ## details.
  13 ##
  14 ## You should have received a copy of the GNU General Public License along with
  15 ## this program; if not, see <http://www.gnu.org/licenses/>.
  16
  17 ## Compute chebyshev type I filter order and cutoff for the desired response
  18 ## characteristics. Rp is the allowable decibels of ripple in the pass
  19 ## band. Rs is the minimum attenuation in the stop band.
  20 ##
  21 ## [n, Wc] = cheb1ord(Wp, Ws, Rp, Rs)
  22 ##     Low pass (Wp<Ws) or high pass (Wp>Ws) filter design.  Wp is the
  23 ##     pass band edge and Ws is the stop band edge.  Frequencies are
  24 ##     normalized to [0,1], corresponding to the range [0,Fs/2].
  25 ##
  26 ## [n, Wc] = cheb1ord([Wp1, Wp2], [Ws1, Ws2], Rp, Rs)
  27 ##     Band pass (Ws1<Wp1<Wp2<Ws2) or band reject (Wp1<Ws1<Ws2<Wp2)
  28 ##     filter design. Wp gives the edges of the pass band, and Ws gives
  29 ##     the edges of the stop band.
  30 ##
  31 ## See also: cheby1
  32
  33 function [n, Wc] = cheb1ord(Wp, Ws, Rp, Rs)
  34
  35   if nargin != 4
  36     print_usage;
  37   elseif length(Wp) != length(Ws)
  38     error("cheb1ord: Wp and Ws must have the same length");
  39   elseif length(Wp) != 1 && length(Wp) != 2
  40     error("cheb1ord: Wp,Ws must have length 1 or 2");
  41   elseif length(Wp) == 2 && ...
  42         (all(Wp>Ws) || all(Ws>Wp) || diff(Wp)<=0 || diff(Ws)<=0)
  43     error("cheb1ord: Wp(1)<Ws(1)<Ws(2)<Wp(2) or Ws(1)<Wp(1)<Wp(2)<Ws(2)");
  44   end
  45
  46   T = 2;
  47
  48   ## returned frequency is the same as the input frequency
  49   Wc = Wp;
  50
  51   ## warp the target frequencies according to the bilinear transform
  52   Ws = (2/T)*tan(pi*Ws./T);
  53   Wp = (2/T)*tan(pi*Wp./T);
  54
  55   if (Wp(1) < Ws(1))
  56     ## low pass
  57     if (length(Wp) == 1)
  58       Wa = Ws/Wp;
  59     else
  60       ## band reject
  61       error ("band reject is not implement yet.");
  62     endif;
  63   else
  64    ## if high pass, reverse the sense of the test
  65    if (length(Wp) == 1)
  66       Wa = Wp/Ws;
  67     else
  68       ## band pass
  69       Wa=(Ws.^2 - Wp(1)*Wp(2))./(Ws*(Wp(1)-Wp(2)));
  70     endif;
  71   endif;
  72   Wa = min(abs(Wa));
  73
  74   ## compute minimum n which satisfies all band edge conditions
  75   stop_atten = 10^(abs(Rs)/10);
  76   pass_atten = 10^(abs(Rp)/10);
  77   n = ceil(acosh(sqrt((stop_atten-1)/(pass_atten-1)))/acosh(Wa));
  78
  79 endfunction
  80
  81 %!demo
  82 %! Fs = 10000;
  83 %! [n, Wc] = cheb1ord (1000/(Fs/2), 1200/(Fs/2), 0.5, 29);
  84 %!
  85 %! subplot (221);
  86 %! axis ([ 0, 1500, -1, 0]);
  87 %! title("Pass band Wp=1000 Rp=0.5");
  88 %! xlabel("Frequency (Hz)");
  89 %! ylabel("Attenuation (dB)");
  90 %! grid;
  91 %! plot ([0, 1000, 1000, 0, 0], [0, 0, -0.5, -0.5, 0], ";;");
  92 %! hold on;
  93 %! [b, a] = cheby1 (n, 0.5, Wc);
  94 %! [h, w] = freqz (b, a, [], Fs);
  95 %! plot (w, 20*log10(abs(h)), ";;");
  96 %! hold off;
  97 %!
  98 %! subplot (222);
  99 %! axis ([ 0, Fs/2, -250, 0]);
 100 %! title("Stop band Ws=1200 Rs=29");
 101 %! xlabel("Frequency (Hz)");
 102 %! ylabel("Attenuation (dB)");
 103 %! grid;
 104 %! plot ([1200, Fs/2, Fs/2, 1200, 1200], [-29, -29, -500, -500, -29], ";;");
 105 %! hold on;
 106 %! [b, a] = cheby1 (n, 0.5, Wc);
 107 %! [h, w] = freqz (b, a, [], Fs);
 108 %! plot (w, 20*log10(abs(h)), ";;");
 109 %! hold off;
 110 %!
 111 %! subplot (223);
 112 %! axis ([ 990, 1010, -0.6, -0.4]);
 113 %! title("Pass band detail Wp=1000 Rp=0.5");
 114 %! xlabel("Frequency (Hz)");
 115 %! ylabel("Attenuation (dB)");
 116 %! grid;
 117 %! plot ([0, 1000, 1000, 0, 0], [0, 0, -0.5, -0.5, 0], ";;");
 118 %! hold on;
 119 %! [b, a] = cheby1 (n, 0.5, Wc);
 120 %! [h, w] = freqz (b, a, [990:1010], Fs);
 121 %! plot (w, 20*log10(abs(h)), ";filter n;");
 122 %! [b, a] = cheby1 (n-1, 0.5, Wc);
 123 %! [h, w] = freqz (b, a, [990:1010], Fs);
 124 %! plot (w, 20*log10(abs(h)), ";filter n-1;");
 125 %! [b, a] = cheby1 (n+1, 0.5, Wc);
 126 %! [h, w] = freqz (b, a, [990:1010], Fs);
 127 %! plot (w, 20*log10(abs(h)), ";filter n+1;");
 128 %! hold off;
 129 %!
 130 %! subplot (224);
 131 %! axis ([ 1190, 1210, -40, -20]);
 132 %! title("Stop band detail Wp=1200 Rp=29");
 133 %! xlabel("Frequency (Hz)");
 134 %! ylabel("Attenuation (dB)");
 135 %! grid;
 136 %! plot ([1200, Fs/2, Fs/2, 1200, 1200], [-29, -29, -500, -500, -29], ";;");
 137 %! hold on;
 138 %! [b, a] = cheby1 (n, 0.5, Wc);
 139 %! [h, w] = freqz (b, a, [1190:1210], Fs);
 140 %! plot (w, 20*log10(abs(h)), ";filter n;");
 141 %! [b, a] = cheby1 (n-1, 0.5, Wc);
 142 %! [h, w] = freqz (b, a, [1190:1210], Fs);
 143 %! plot (w, 20*log10(abs(h)), ";filter n-1;");
 144 %! [b, a] = cheby1 (n+1, 0.5, Wc);
 145 %! [h, w] = freqz (b, a, [1190:1210], Fs);
 146 %! plot (w, 20*log10(abs(h)), ";filter n+1;");
 147 %! hold off;