/usr/share/octave/packages/signal-1.3.0/cheb2ord.m is in octave-signal 1.3.0-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 | ## Copyright (C) 2000 Paul Kienzle <pkienzle@users.sf.net>
##
## 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/>.
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{n}, @var{Wc}] =} cheb2ord (@var{Wp}, @var{Ws}, @var{Rp}, @var{Rs})
## @deftypefnx {Function File} {[@var{n}, @var{Wc}] =} cheb2ord ([@var{Wp1}, @var{Wp2}], [@var{Ws1}, @var{Ws2}], @var{Rp}, @var{Rs})
## Compute chebyshev type II 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] = cheb2ord(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] = cheb2ord([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:
##
## @seealso{cheby2}
## @end deftypefn
function [n, Wc] = cheb2ord(Wp, Ws, Rp, Rs)
if nargin != 4
print_usage;
elseif length(Wp) != length(Ws)
error("cheb2ord: Wp and Ws must have the same length");
elseif length(Wp) != 1 && length(Wp) != 2
error("cheb2ord: Wp,Ws must have length 1 or 2");
elseif length(Wp) == 2 && ...
(all(Wp>Ws) || all(Ws>Wp) || diff(Wp)<=0 || diff(Ws)<=0)
error("cheb2ord: Wp(1)<Ws(1)<Ws(2)<Wp(2) or Ws(1)<Wp(1)<Wp(2)<Ws(2)");
endif
T = 2;
## returned frequency is the same as the input frequency
Wc = Ws;
## warp the target frequencies according to the bilinear transform
Ws = (2/T)*tan(pi*Ws./T);
Wp = (2/T)*tan(pi*Wp./T);
if (Wp(1) < Ws(1))
## low pass
if (length(Wp) == 1)
Wa = Wp/Ws;
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 = Ws/Wp;
else
## band pass
Wa=(Wp.^2 - Ws(1)*Ws(2))./(Wp*(Ws(1)-Ws(2)));
endif;
endif;
Wa = min(abs(Wa));
## compute minimum n which satisfies all band edge conditions
stop_atten = 10^(abs(Rs)/10);
pass_atten = 10^(abs(Rp)/10);
n = ceil(acosh(sqrt((stop_atten-1)/(pass_atten-1)))/acosh(1/Wa));
endfunction
%!demo
%! Fs = 10000;
%! [n, Wc] = cheb2ord (1000/(Fs/2), 1200/(Fs/2), 0.5, 29);
%!
%! subplot (221);
%! plot ([0, 1000, 1000, 0, 0], [0, 0, -0.5, -0.5, 0], ";;");
%! hold on;
%! grid;
%! title("Pass band Wp=1000 Rp=0.0");
%! xlabel("Frequency (Hz)");
%! ylabel("Attenuation (dB)");
%! [b, a] = cheby2 (n, 29, Wc);
%! [h, w] = freqz (b, a, [], Fs);
%! plot (w, 20*log10(abs(h)), ";;");
%! axis ([ 0, 1500, -1, 0]);
%! hold off;
%!
%! subplot (222);
%! plot ([1200, Fs/2, Fs/2, 1200, 1200], [-29, -29, -500, -500, -29], ";;");
%! hold on;
%! axis ([ 0, Fs/2, -250, 0]);
%! title("Stop band Ws=1200 Rs=29");
%! xlabel("Frequency (Hz)");
%! ylabel("Attenuation (dB)");
%! grid;
%! [b, a] = cheby2 (n, 29, Wc);
%! [h, w] = freqz (b, a, [], Fs);
%! plot (w, 20*log10(abs(h)), ";;");
%! hold off;
%!
%! subplot (223);
%! plot ([0, 1000, 1000, 0, 0], [0, 0, -0.5, -0.5, 0], ";;");
%! hold on;
%! axis ([ 800, 1010, -0.6, -0.0]);
%! title("Pass band detail Wp=1000 Rp=0.5");
%! xlabel("Frequency (Hz)");
%! ylabel("Attenuation (dB)");
%! grid;
%! [b, a] = cheby2 (n, 29, Wc);
%! [h, w] = freqz (b, a, [800:1010], Fs);
%! plot (w, 20*log10(abs(h)), "r;filter n;");
%! [b, a] = cheby2 (n-1, 29, Wc);
%! [h, w] = freqz (b, a, [800:1010], Fs);
%! plot (w, 20*log10(abs(h)), "b;filter n-1;");
%! [b, a] = cheby2 (n+1, 29, Wc);
%! [h, w] = freqz (b, a, [800:1010], Fs);
%! plot (w, 20*log10(abs(h)), "g;filter n+1;");
%! hold off;
%!
%! subplot (224);
%! plot ([1200, Fs/2, Fs/2, 1200, 1200], [-29, -29, -500, -500, -29], ";;");
%! hold on;
%! axis ([ 1190, 1210, -40, -20]);
%! title("Stop band detail Wp=1200 Rp=29");
%! xlabel("Frequency (Hz)");
%! ylabel("Attenuation (dB)");
%! grid;
%! [b, a] = cheby2 (n, 29, Wc);
%! [h, w] = freqz (b, a, [1190:1210], Fs);
%! plot (w, 20*log10(abs(h)), "r;filter n;");
%! [b, a] = cheby2 (n-1, 29, Wc);
%! [h, w] = freqz (b, a, [1190:1210], Fs);
%! plot (w, 20*log10(abs(h)), "b;filter n-1;");
%! [b, a] = cheby2 (n+1, 29, Wc);
%! [h, w] = freqz (b, a, [1190:1210], Fs);
%! plot (w, 20*log10(abs(h)), "g;filter n+1;");
%! hold off;
|