octave_packages/signal-1.1.3/cheb2ord.m

   1 ## Copyright (C) 2000 Paul Kienzle <pkienzle@users.sf.net>
   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 ## Compute chebyshev type II filter order and cutoff for the desired response
  17 ## characteristics. Rp is the allowable decibels of ripple in the pass
  18 ## band. Rs is the minimum attenuation in the stop band.
  19 ##
  20 ## [n, Wc] = cheb2ord(Wp, Ws, Rp, Rs)
  21 ##     Low pass (Wp<Ws) or high pass (Wp>Ws) filter design.  Wp is the
  22 ##     pass band edge and Ws is the stop band edge.  Frequencies are
  23 ##     normalized to [0,1], corresponding to the range [0,Fs/2].
  24 ##
  25 ## [n, Wc] = cheb2ord([Wp1, Wp2], [Ws1, Ws2], Rp, Rs)
  26 ##     Band pass (Ws1<Wp1<Wp2<Ws2) or band reject (Wp1<Ws1<Ws2<Wp2)
  27 ##     filter design. Wp gives the edges of the pass band, and Ws gives
  28 ##     the edges of the stop band.
  29 ##
  30 ## Theory:
  31 ##
  32 ## See also: cheby2
  33
  34 function [n, Wc] = cheb2ord(Wp, Ws, Rp, Rs)
  35
  36   if nargin != 4
  37     print_usage;
  38   elseif length(Wp) != length(Ws)
  39     error("cheb2ord: Wp and Ws must have the same length");
  40   elseif length(Wp) != 1 && length(Wp) != 2
  41     error("cheb2ord: Wp,Ws must have length 1 or 2");
  42   elseif length(Wp) == 2 && ...
  43           (all(Wp>Ws) || all(Ws>Wp) || diff(Wp)<=0 || diff(Ws)<=0)
  44     error("cheb2ord: Wp(1)<Ws(1)<Ws(2)<Wp(2) or Ws(1)<Wp(1)<Wp(2)<Ws(2)");
  45   end
  46
  47   T = 2;
  48
  49   ## returned frequency is the same as the input frequency
  50   Wc = Ws;
  51
  52   ## warp the target frequencies according to the bilinear transform
  53   Ws = (2/T)*tan(pi*Ws./T);
  54   Wp = (2/T)*tan(pi*Wp./T);
  55
  56   if (Wp(1) < Ws(1))
  57     ## low pass
  58     if (length(Wp) == 1)
  59       Wa = Wp/Ws;
  60     else
  61       ## band reject
  62       error ("band reject is not implement yet.");
  63     endif;
  64   else
  65    ## if high pass, reverse the sense of the test
  66    if (length(Wp) == 1)
  67       Wa = Ws/Wp;
  68     else
  69       ## band pass
  70       Wa=(Wp.^2 - Ws(1)*Ws(2))./(Wp*(Ws(1)-Ws(2)));
  71     endif;
  72   endif;
  73   Wa = min(abs(Wa));
  74
  75   ## compute minimum n which satisfies all band edge conditions
  76   stop_atten = 10^(abs(Rs)/10);
  77   pass_atten = 10^(abs(Rp)/10);
  78   n = ceil(acosh(sqrt((stop_atten-1)/(pass_atten-1)))/acosh(1/Wa));
  79
  80 endfunction
  81
  82 %!demo
  83 %! Fs = 10000;
  84 %! [n, Wc] = cheb2ord (1000/(Fs/2), 1200/(Fs/2), 0.5, 29);
  85 %!
  86 %! subplot (221);
  87 %! plot ([0, 1000, 1000, 0, 0], [0, 0, -0.5, -0.5, 0], ";;");
  88 %! hold on;
  89 %! grid;
  90 %! title("Pass band Wp=1000 Rp=0.0");
  91 %! xlabel("Frequency (Hz)");
  92 %! ylabel("Attenuation (dB)");
  93 %! [b, a] = cheby2 (n, 29, Wc);
  94 %! [h, w] = freqz (b, a, [], Fs);
  95 %! plot (w, 20*log10(abs(h)), ";;");
  96 %! axis ([ 0, 1500, -1, 0]);
  97 %! hold off;
  98 %!
  99 %! subplot (222);
 100 %! plot ([1200, Fs/2, Fs/2, 1200, 1200], [-29, -29, -500, -500, -29], ";;");
 101 %! hold on;
 102 %! axis ([ 0, Fs/2, -250, 0]);
 103 %! title("Stop band Ws=1200 Rs=29");
 104 %! xlabel("Frequency (Hz)");
 105 %! ylabel("Attenuation (dB)");
 106 %! grid;
 107 %! [b, a] = cheby2 (n, 29, Wc);
 108 %! [h, w] = freqz (b, a, [], Fs);
 109 %! plot (w, 20*log10(abs(h)), ";;");
 110 %! hold off;
 111 %!
 112 %! subplot (223);
 113 %! plot ([0, 1000, 1000, 0, 0], [0, 0, -0.5, -0.5, 0], ";;");
 114 %! hold on;
 115 %! axis ([ 800, 1010, -0.6, -0.0]);
 116 %! title("Pass band detail Wp=1000 Rp=0.5");
 117 %! xlabel("Frequency (Hz)");
 118 %! ylabel("Attenuation (dB)");
 119 %! grid;
 120 %! [b, a] = cheby2 (n, 29, Wc);
 121 %! [h, w] = freqz (b, a, [800:1010], Fs);
 122 %! plot (w, 20*log10(abs(h)), "r;filter n;");
 123 %! [b, a] = cheby2 (n-1, 29, Wc);
 124 %! [h, w] = freqz (b, a, [800:1010], Fs);
 125 %! plot (w, 20*log10(abs(h)), "b;filter n-1;");
 126 %! [b, a] = cheby2 (n+1, 29, Wc);
 127 %! [h, w] = freqz (b, a, [800:1010], Fs);
 128 %! plot (w, 20*log10(abs(h)), "g;filter n+1;");
 129 %! hold off;
 130 %!
 131 %! subplot (224);
 132 %! plot ([1200, Fs/2, Fs/2, 1200, 1200], [-29, -29, -500, -500, -29], ";;");
 133 %! hold on;
 134 %! axis ([ 1190, 1210, -40, -20]);
 135 %! title("Stop band detail Wp=1200 Rp=29");
 136 %! xlabel("Frequency (Hz)");
 137 %! ylabel("Attenuation (dB)");
 138 %! grid;
 139 %! [b, a] = cheby2 (n, 29, Wc);
 140 %! [h, w] = freqz (b, a, [1190:1210], Fs);
 141 %! plot (w, 20*log10(abs(h)), "r;filter n;");
 142 %! [b, a] = cheby2 (n-1, 29, Wc);
 143 %! [h, w] = freqz (b, a, [1190:1210], Fs);
 144 %! plot (w, 20*log10(abs(h)), "b;filter n-1;");
 145 %! [b, a] = cheby2 (n+1, 29, Wc);
 146 %! [h, w] = freqz (b, a, [1190:1210], Fs);
 147 %! plot (w, 20*log10(abs(h)), "g;filter n+1;");
 148 %! hold off;