/usr/share/octave/packages/signal-1.3.2/cheb1ord.m is in octave-signal 1.3.2-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 153 154 155 156 157 158 159 160 161 162 163 164 165 166 | ## Copyright (C) 2000 Paul Kienzle <pkienzle@users.sf.net>
## Copyright (C) 2000 Laurent S. Mazet
##
## 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} =} cheb1ord (@var{wp}, @var{ws}, @var{rp}, @var{rs})
## @deftypefnx {Function File} {@var{n} =} cheb1ord ([@var{wp1}, @var{wp2}], [@var{ws1}, @var{ws2}], @var{rp}, @var{rs})
## @deftypefnx {Function File} {[@var{n}, @var{wc}] =} cheb1ord (@dots{})
## Compute the minimum filter order of a Chebyshev type I filter with the
## desired response characteristics. The filter frequency band edges are
## specified by the passband frequency @var{wp} and stopband frequency @var{ws}.
## Frequencies are normalized to the Nyquist frequency in the range [0,1].
## @var{rp} is the allowable passband ripple measured in decibels, and @var{rs}
## is the minimum attenuation in the stop band, also in decibels. The output
## arguments @var{n} and @var{wc} can be given as inputs to @code{cheby1}.
##
## If @var{wp} and @var{ws} are scalars, then @var{wp} is the passband cutoff
## frequency and @var{ws} is the stopband edge frequency. If @var{ws} is
## greater than @var{wp}, the filter is a low-pass filter. If @var{wp} is
## greater than @var{ws}, the filter is a high-pass filter.
##
## If @var{wp} and @var{ws} are vectors of length 2, then @var{wp} defines the
## passband interval and @var{ws} defines the stopband interval. If @var{wp}
## is contained within @var{ws} (@var{ws1} < @var{wp1} < @var{wp2} < @var{ws2}),
## the filter is a band-pass filter. If @var{ws} is contained within @var{wp}
## (@var{wp1} < @var{ws1} < @var{ws2} < @var{wp2}), the filter is a band-stop
## or band-reject filter.
##
## @seealso{buttord, cheby1, cheb2ord, ellipord}
## @end deftypefn
function [n, Wc] = cheb1ord(Wp, Ws, Rp, Rs)
if nargin != 4
print_usage;
else
validate_filter_bands ("cheb1ord", Wp, Ws);
endif
T = 2;
## returned frequency is the same as the input frequency
Wc = Wp;
## 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 = Ws/Wp;
else
## FIXME: Implement band reject filter type
error ("cheb1ord: band reject is not yet implemented");
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
stop_atten = 10^(abs(Rs)/10);
pass_atten = 10^(abs(Rp)/10);
n = ceil(acosh(sqrt((stop_atten-1)/(pass_atten-1)))/acosh(Wa));
endfunction
%!demo
%! Fs = 10000;
%! [n, Wc] = cheb1ord (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;
%! axis ([ 0, 1500, -1, 0]);
%! title("Pass band Wp=1000 Rp=0.5");
%! xlabel("Frequency (Hz)");
%! ylabel("Attenuation (dB)");
%! grid;
%! [b, a] = cheby1 (n, 0.5, Wc);
%! [h, w] = freqz (b, a, [], Fs);
%! plot (w, 20*log10(abs(h)), ";;");
%! 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] = cheby1 (n, 0.5, 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 ([ 990, 1010, -0.6, -0.4]);
%! title("Pass band detail Wp=1000 Rp=0.5");
%! xlabel("Frequency (Hz)");
%! ylabel("Attenuation (dB)");
%! grid;
%! [b, a] = cheby1 (n, 0.5, Wc);
%! [h, w] = freqz (b, a, [990:1010], Fs);
%! plot (w, 20*log10(abs(h)), ";filter n;");
%! [b, a] = cheby1 (n-1, 0.5, Wc);
%! [h, w] = freqz (b, a, [990:1010], Fs);
%! plot (w, 20*log10(abs(h)), ";filter n-1;");
%! [b, a] = cheby1 (n+1, 0.5, Wc);
%! [h, w] = freqz (b, a, [990:1010], Fs);
%! plot (w, 20*log10(abs(h)), ";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] = cheby1 (n, 0.5, Wc);
%! [h, w] = freqz (b, a, [1190:1210], Fs);
%! plot (w, 20*log10(abs(h)), ";filter n;");
%! [b, a] = cheby1 (n-1, 0.5, Wc);
%! [h, w] = freqz (b, a, [1190:1210], Fs);
%! plot (w, 20*log10(abs(h)), ";filter n-1;");
%! [b, a] = cheby1 (n+1, 0.5, Wc);
%! [h, w] = freqz (b, a, [1190:1210], Fs);
%! plot (w, 20*log10(abs(h)), ";filter n+1;");
%! hold off;
%% Test case from demo
%!assert (cheb1ord (0.2, 0.24, 0.5, 29), 8)
%% Test input validation
%!error cheb1ord ()
%!error cheb1ord (.1)
%!error cheb1ord (.1, .2)
%!error cheb1ord (.1, .2, 3)
%!error cheb1ord (.1, .2, 3, 4, 5)
%!error cheb1ord ([.1 .1], [.2 .2], 3, 4)
%!error cheb1ord ([.1 .2], [.5 .6], 3, 4)
%!error cheb1ord ([.1 .5], [.2 .6], 3, 4)
|