Mercurial > hg > octave-nkf
annotate scripts/polynomial/residue.m @ 14363:f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
* wavread.m, acosd.m, acot.m, acotd.m, acoth.m, acsc.m, acscd.m, acsch.m,
asec.m, asecd.m, asech.m, asind.m, atand.m, cosd.m, cot.m, cotd.m, coth.m,
csc.m, cscd.m, csch.m, sec.m, secd.m, sech.m, sind.m, tand.m, accumarray.m,
accumdim.m, bitcmp.m, bitget.m, bitset.m, blkdiag.m, cart2pol.m, cart2sph.m,
celldisp.m, chop.m, circshift.m, colon.m, common_size.m, cplxpair.m,
cumtrapz.m, curl.m, dblquad.m, deal.m, divergence.m, flipdim.m, fliplr.m,
flipud.m, genvarname.m, gradient.m, idivide.m, int2str.m, interp1.m,
interp1q.m, interp2.m, interp3.m, interpft.m, interpn.m, isa.m, isdir.m,
isequal.m, isequalwithequalnans.m, issquare.m, logspace.m, nargchk.m,
narginchk.m, nargoutchk.m, nextpow2.m, nthargout.m, num2str.m, pol2cart.m,
polyarea.m, postpad.m, prepad.m, profile.m, profshow.m, quadgk.m, quadv.m,
randi.m, rat.m, repmat.m, rot90.m, rotdim.m, shift.m, shiftdim.m, sph2cart.m,
structfun.m, trapz.m, triplequad.m, convhull.m, dsearch.m, dsearchn.m,
griddata3.m, griddatan.m, rectint.m, tsearchn.m, __makeinfo__.m, doc.m,
get_first_help_sentence.m, help.m, type.m, unimplemented.m, which.m, imread.m,
imwrite.m, dlmwrite.m, fileread.m, is_valid_file_id.m, strread.m, textread.m,
textscan.m, commutation_matrix.m, cond.m, condest.m, cross.m,
duplication_matrix.m, expm.m, housh.m, isdefinite.m, ishermitian.m,
issymmetric.m, logm.m, normest.m, null.m, onenormest.m, orth.m, planerot.m,
qzhess.m, rank.m, rref.m, trace.m, vech.m, ans.m, bincoeff.m, bug_report.m,
bzip2.m, comma.m, compare_versions.m, computer.m, edit.m, fileparts.m,
fullfile.m, getfield.m, gzip.m, info.m, inputname.m, isappdata.m, isdeployed.m,
ismac.m, ispc.m, isunix.m, list_primes.m, ls.m, mexext.m, namelengthmax.m,
news.m, orderfields.m, paren.m, recycle.m, rmappdata.m, semicolon.m,
setappdata.m, setfield.m, substruct.m, symvar.m, ver.m, version.m,
warning_ids.m, xor.m, fminbnd.m, fsolve.m, fzero.m, lsqnonneg.m, optimset.m,
pqpnonneg.m, sqp.m, matlabroot.m, __gnuplot_drawnow__.m,
__plt_get_axis_arg__.m, ancestor.m, cla.m, clf.m, close.m, colorbar.m,
colstyle.m, comet3.m, contourc.m, figure.m, gca.m, gcbf.m, gcbo.m, gcf.m,
ginput.m, graphics_toolkit.m, gtext.m, hggroup.m, hist.m, hold.m, isfigure.m,
ishghandle.m, ishold.m, isocolors.m, isonormals.m, isosurface.m, isprop.m,
legend.m, line.m, loglog.m, loglogerr.m, meshgrid.m, ndgrid.m, newplot.m,
orient.m, patch.m, plot3.m, plotyy.m, __print_parse_opts__.m, quiver3.m,
refreshdata.m, ribbon.m, semilogx.m, semilogxerr.m, semilogy.m, stem.m,
stem3.m, subplot.m, title.m, uigetfile.m, view.m, whitebg.m, compan.m, conv.m,
deconv.m, mkpp.m, mpoles.m, pchip.m, poly.m, polyaffine.m, polyder.m,
polyfit.m, polygcd.m, polyint.m, polyout.m, polyval.m, polyvalm.m, ppder.m,
ppint.m, ppjumps.m, ppval.m, residue.m, roots.m, spline.m, intersect.m,
ismember.m, powerset.m, setdiff.m, setxor.m, union.m, unique.m,
autoreg_matrix.m, bartlett.m, blackman.m, detrend.m, fftconv.m, fftfilt.m,
fftshift.m, freqz.m, hamming.m, hanning.m, ifftshift.m, sinc.m, sinetone.m,
sinewave.m, unwrap.m, bicg.m, bicgstab.m, gmres.m, gplot.m, nonzeros.m, pcg.m,
pcr.m, spaugment.m, spconvert.m, spdiags.m, speye.m, spfun.m, spones.m,
sprand.m, sprandsym.m, spstats.m, spy.m, svds.m, treelayout.m, bessel.m,
beta.m, betaln.m, factor.m, factorial.m, isprime.m, lcm.m, legendre.m,
nchoosek.m, nthroot.m, perms.m, pow2.m, primes.m, reallog.m, realpow.m,
realsqrt.m, hadamard.m, hankel.m, hilb.m, invhilb.m, magic.m, rosser.m,
vander.m, __finish__.m, center.m, cloglog.m, corr.m, cov.m, gls.m, histc.m,
iqr.m, kendall.m, kurtosis.m, logit.m, mahalanobis.m, mean.m, meansq.m,
median.m, mode.m, moment.m, ols.m, ppplot.m, prctile.m, probit.m, quantile.m,
range.m, ranks.m, run_count.m, runlength.m, skewness.m, spearman.m,
statistics.m, std.m, table.m, var.m, zscore.m, betacdf.m, betainv.m, betapdf.m,
betarnd.m, binocdf.m, binoinv.m, binopdf.m, binornd.m, cauchy_cdf.m,
cauchy_inv.m, cauchy_pdf.m, cauchy_rnd.m, chi2cdf.m, chi2inv.m, chi2pdf.m,
chi2rnd.m, discrete_cdf.m, discrete_inv.m, discrete_pdf.m, discrete_rnd.m,
empirical_cdf.m, empirical_inv.m, empirical_pdf.m, empirical_rnd.m, expcdf.m,
expinv.m, exppdf.m, exprnd.m, fcdf.m, finv.m, fpdf.m, frnd.m, gamcdf.m,
gaminv.m, gampdf.m, gamrnd.m, geocdf.m, geoinv.m, geopdf.m, geornd.m,
hygecdf.m, hygeinv.m, hygepdf.m, hygernd.m, kolmogorov_smirnov_cdf.m,
laplace_cdf.m, laplace_inv.m, laplace_pdf.m, laplace_rnd.m, logistic_cdf.m,
logistic_inv.m, logistic_pdf.m, logistic_rnd.m, logncdf.m, logninv.m,
lognpdf.m, lognrnd.m, nbincdf.m, nbininv.m, nbinpdf.m, nbinrnd.m, normcdf.m,
norminv.m, normpdf.m, normrnd.m, poisscdf.m, poissinv.m, poisspdf.m,
poissrnd.m, stdnormal_cdf.m, stdnormal_inv.m, stdnormal_pdf.m, stdnormal_rnd.m,
tcdf.m, tinv.m, tpdf.m, trnd.m, unidcdf.m, unidinv.m, unidpdf.m, unidrnd.m,
unifcdf.m, unifinv.m, unifpdf.m, unifrnd.m, wblcdf.m, wblinv.m, wblpdf.m,
wblrnd.m, kolmogorov_smirnov_test.m, kruskal_wallis_test.m, base2dec.m,
bin2dec.m, blanks.m, cstrcat.m, deblank.m, dec2base.m, dec2bin.m, dec2hex.m,
findstr.m, hex2dec.m, index.m, isletter.m, mat2str.m, rindex.m, str2num.m,
strcat.m, strjust.m, strmatch.m, strsplit.m, strtok.m, strtrim.m, strtrunc.m,
substr.m, validatestring.m, demo.m, example.m, fail.m, speed.m, addtodate.m,
asctime.m, clock.m, ctime.m, date.m, datenum.m, datetick.m, datevec.m,
eomday.m, etime.m, is_leap_year.m, now.m:
Use Octave coding conventions in all m-file %!test blocks
| author | Rik <octave@nomad.inbox5.com> |
|---|---|
| date | Mon, 13 Feb 2012 07:29:44 -0800 |
| parents | 4d917a6a858b |
| children | 86854d032a37 |
| rev | line source |
|---|---|
|
14138
72c96de7a403
maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents:
14104
diff
changeset
|
1 ## Copyright (C) 1994-2012 John W. Eaton |
| 6964 | 2 ## Copyright (C) 2007 Ben Abbott |
| 2313 | 3 ## |
| 4 ## This file is part of Octave. | |
| 5 ## | |
| 6 ## Octave is free software; you can redistribute it and/or modify it | |
| 7 ## under the terms of the GNU General Public License as published by | |
| 7016 | 8 ## the Free Software Foundation; either version 3 of the License, or (at |
| 9 ## your option) any later version. | |
| 2313 | 10 ## |
| 11 ## Octave is distributed in the hope that it will be useful, but | |
| 12 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
| 14 ## General Public License for more details. | |
| 15 ## | |
| 16 ## You should have received a copy of the GNU General Public License | |
| 7016 | 17 ## along with Octave; see the file COPYING. If not, see |
| 18 ## <http://www.gnu.org/licenses/>. | |
| 904 | 19 |
| 3368 | 20 ## -*- texinfo -*- |
|
13929
9cae456085c2
Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents:
13128
diff
changeset
|
21 ## @deftypefn {Function File} {[@var{r}, @var{p}, @var{k}, @var{e}] =} residue (@var{b}, @var{a}) |
|
13127
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
22 ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}) |
|
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
23 ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}, @var{e}) |
|
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
24 ## The first calling form computes the partial fraction expansion for the |
|
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
25 ## quotient of the polynomials, @var{b} and @var{a}. |
| 3368 | 26 ## @tex |
| 27 ## $$ | |
| 6978 | 28 ## {B(s)\over A(s)} = \sum_{m=1}^M {r_m\over (s-p_m)^e_m} |
| 3368 | 29 ## + \sum_{i=1}^N k_i s^{N-i}. |
| 30 ## $$ | |
| 31 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
32 ## @ifnottex |
| 3426 | 33 ## |
| 3368 | 34 ## @example |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
35 ## @group |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
36 ## B(s) M r(m) N |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
37 ## ---- = SUM ------------- + SUM k(i)*s^(N-i) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
38 ## A(s) m=1 (s-p(m))^e(m) i=1 |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
39 ## @end group |
| 3368 | 40 ## @end example |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
41 ## |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
42 ## @end ifnottex |
| 3368 | 43 ## @noindent |
| 3499 | 44 ## where @math{M} is the number of poles (the length of the @var{r}, |
| 6978 | 45 ## @var{p}, and @var{e}), the @var{k} vector is a polynomial of order @math{N-1} |
| 46 ## representing the direct contribution, and the @var{e} vector specifies | |
| 8494 | 47 ## the multiplicity of the m-th residue's pole. |
| 3426 | 48 ## |
| 6964 | 49 ## For example, |
| 3426 | 50 ## |
| 3368 | 51 ## @example |
| 52 ## @group | |
| 6964 | 53 ## b = [1, 1, 1]; |
| 54 ## a = [1, -5, 8, -4]; | |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
55 ## [r, p, k, e] = residue (b, a) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
56 ## @result{} r = [-2; 7; 3] |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
57 ## @result{} p = [2; 2; 1] |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
58 ## @result{} k = [](0x0) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
59 ## @result{} e = [1; 2; 1] |
| 3368 | 60 ## @end group |
| 61 ## @end example | |
| 3426 | 62 ## |
| 3368 | 63 ## @noindent |
| 6978 | 64 ## which represents the following partial fraction expansion |
| 3368 | 65 ## @tex |
| 66 ## $$ | |
| 67 ## {s^2+s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} | |
| 68 ## $$ | |
| 69 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
70 ## @ifnottex |
| 3426 | 71 ## |
| 3368 | 72 ## @example |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
73 ## @group |
| 3368 | 74 ## s^2 + s + 1 -2 7 3 |
| 2311 | 75 ## ------------------- = ----- + ------- + ----- |
| 76 ## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1) | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
77 ## @end group |
| 3368 | 78 ## @end example |
| 6964 | 79 ## |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
80 ## @end ifnottex |
| 6978 | 81 ## |
|
13127
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
82 ## The second calling form performs the inverse operation and computes |
|
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
83 ## the reconstituted quotient of polynomials, @var{b}(s)/@var{a}(s), |
|
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
84 ## from the partial fraction expansion; represented by the residues, |
|
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
85 ## poles, and a direct polynomial specified by @var{r}, @var{p} and |
|
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
86 ## @var{k}, and the pole multiplicity @var{e}. |
| 7011 | 87 ## |
| 88 ## If the multiplicity, @var{e}, is not explicitly specified the multiplicity is | |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
89 ## determined by the function @code{mpoles}. |
| 6978 | 90 ## |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
91 ## For example: |
| 6964 | 92 ## |
| 93 ## @example | |
| 94 ## @group | |
| 7011 | 95 ## r = [-2; 7; 3]; |
| 96 ## p = [2; 2; 1]; | |
| 97 ## k = [1, 0]; | |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
98 ## [b, a] = residue (r, p, k) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
99 ## @result{} b = [1, -5, 9, -3, 1] |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
100 ## @result{} a = [1, -5, 8, -4] |
| 7011 | 101 ## |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
102 ## where mpoles is used to determine e = [1; 2; 1] |
| 7011 | 103 ## @end group |
| 104 ## @end example | |
| 105 ## | |
| 106 ## Alternatively the multiplicity may be defined explicitly, for example, | |
| 107 ## | |
| 108 ## @example | |
| 109 ## @group | |
| 110 ## r = [7; 3; -2]; | |
| 111 ## p = [2; 1; 2]; | |
| 112 ## k = [1, 0]; | |
| 113 ## e = [2; 1; 1]; | |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
114 ## [b, a] = residue (r, p, k, e) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
115 ## @result{} b = [1, -5, 9, -3, 1] |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
116 ## @result{} a = [1, -5, 8, -4] |
| 6964 | 117 ## @end group |
| 118 ## @end example | |
| 119 ## | |
| 120 ## @noindent | |
| 6978 | 121 ## which represents the following partial fraction expansion |
| 6964 | 122 ## @tex |
| 123 ## $$ | |
| 6978 | 124 ## {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} + s = {s^4-5s^3+9s^2-3s+1\over s^3-5s^2+8s-4} |
| 6964 | 125 ## $$ |
| 126 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
127 ## @ifnottex |
| 6964 | 128 ## |
| 129 ## @example | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
130 ## @group |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
131 ## -2 7 3 s^4 - 5s^3 + 9s^2 - 3s + 1 |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
132 ## ----- + ------- + ----- + s = -------------------------- |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
133 ## (s-2) (s-2)^2 (s-1) s^3 - 5s^2 + 8s - 4 |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
134 ## @end group |
| 6964 | 135 ## @end example |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
136 ## |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
137 ## @end ifnottex |
|
14104
614505385171
doc: Overhaul docstrings for polynomial functions.
Rik <octave@nomad.inbox5.com>
parents:
13963
diff
changeset
|
138 ## @seealso{mpoles, poly, roots, conv, deconv} |
| 3368 | 139 ## @end deftypefn |
| 1025 | 140 |
| 3202 | 141 ## Author: Tony Richardson <arichard@stark.cc.oh.us> |
| 6964 | 142 ## Author: Ben Abbott <bpabbott@mac.com> |
| 2312 | 143 ## Created: June 1994 |
| 144 ## Adapted-By: jwe | |
| 559 | 145 |
| 6978 | 146 function [r, p, k, e] = residue (b, a, varargin) |
| 559 | 147 |
| 7011 | 148 if (nargin < 2 || nargin > 4) |
| 6046 | 149 print_usage (); |
| 559 | 150 endif |
| 151 | |
| 6964 | 152 toler = .001; |
| 153 | |
| 7011 | 154 if (nargin >= 3) |
| 155 if (nargin >= 4) | |
| 156 e = varargin{2}; | |
| 157 else | |
| 158 e = []; | |
| 159 endif | |
| 6964 | 160 ## The inputs are the residue, pole, and direct part. Solve for the |
| 161 ## corresponding numerator and denominator polynomials | |
| 7011 | 162 [r, p] = rresidue (b, a, varargin{1}, toler, e); |
| 6964 | 163 return |
| 7011 | 164 endif |
| 559 | 165 |
| 2303 | 166 ## Make sure both polynomials are in reduced form. |
| 1025 | 167 |
| 168 a = polyreduce (a); | |
| 169 b = polyreduce (b); | |
| 559 | 170 |
| 1025 | 171 b = b / a(1); |
| 172 a = a / a(1); | |
| 559 | 173 |
| 1025 | 174 la = length (a); |
| 175 lb = length (b); | |
| 559 | 176 |
| 2303 | 177 ## Handle special cases here. |
| 1025 | 178 |
| 179 if (la == 0 || lb == 0) | |
| 559 | 180 k = r = p = e = []; |
| 181 return; | |
| 182 elseif (la == 1) | |
| 1025 | 183 k = b / a; |
| 559 | 184 r = p = e = []; |
| 185 return; | |
| 186 endif | |
| 187 | |
| 2303 | 188 ## Find the poles. |
| 1025 | 189 |
| 190 p = roots (a); | |
| 191 lp = length (p); | |
| 559 | 192 |
| 6964 | 193 ## Sort poles so that multiplicity loop will work. |
| 194 | |
| 195 [e, indx] = mpoles (p, toler, 1); | |
| 196 p = p (indx); | |
| 559 | 197 |
| 7398 | 198 ## For each group of pole multiplicity, set the value of each |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
199 ## pole to the average of the group. This reduces the error in |
| 7398 | 200 ## the resulting poles. |
| 201 | |
| 202 p_group = cumsum (e == 1); | |
| 203 for ng = 1:p_group(end) | |
| 204 m = find (p_group == ng); | |
| 205 p(m) = mean (p(m)); | |
| 206 endfor | |
| 207 | |
| 2303 | 208 ## Find the direct term if there is one. |
| 1025 | 209 |
| 210 if (lb >= la) | |
| 6964 | 211 ## Also return the reduced numerator. |
| 1025 | 212 [k, b] = deconv (b, a); |
| 213 lb = length (b); | |
| 559 | 214 else |
| 215 k = []; | |
| 216 endif | |
| 217 | |
| 7398 | 218 ## Determine if the poles are (effectively) zero. |
| 219 | |
| 220 small = max (abs (p)); | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
221 if (isa (a, "single") || isa (b, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
222 small = max ([small, 1]) * eps ("single") * 1e4 * (1 + numel (p))^2; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
223 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
224 small = max ([small, 1]) * eps * 1e4 * (1 + numel (p))^2; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
225 endif |
| 7398 | 226 p(abs (p) < small) = 0; |
| 227 | |
| 228 ## Determine if the poles are (effectively) real, or imaginary. | |
| 229 | |
| 230 index = (abs (imag (p)) < small); | |
| 231 p(index) = real (p(index)); | |
| 232 index = (abs (real (p)) < small); | |
| 233 p(index) = 1i * imag (p(index)); | |
| 234 | |
| 235 ## The remainder determines the residues. The case of one pole | |
| 236 ## is trivial. | |
| 237 | |
| 1025 | 238 if (lp == 1) |
| 239 r = polyval (b, p); | |
| 559 | 240 return; |
| 241 endif | |
| 242 | |
| 6964 | 243 ## Determine the order of the denominator and remaining numerator. |
| 244 ## With the direct term removed the potential order of the numerator | |
| 245 ## is one less than the order of the denominator. | |
| 1025 | 246 |
| 6964 | 247 aorder = numel (a) - 1; |
| 248 border = aorder - 1; | |
| 1025 | 249 |
| 6964 | 250 ## Construct a system of equations relating the individual |
| 251 ## contributions from each residue to the complete numerator. | |
| 559 | 252 |
| 6964 | 253 A = zeros (border+1, border+1); |
| 254 B = prepad (reshape (b, [numel(b), 1]), border+1, 0); | |
| 255 for ip = 1:numel(p) | |
| 256 ri = zeros (size (p)); | |
| 257 ri(ip) = 1; | |
| 258 A(:,ip) = prepad (rresidue (ri, p, [], toler), border+1, 0).'; | |
| 259 endfor | |
| 559 | 260 |
| 2303 | 261 ## Solve for the residues. |
| 1025 | 262 |
| 263 r = A \ B; | |
| 559 | 264 |
| 265 endfunction | |
| 6964 | 266 |
| 7011 | 267 function [pnum, pden, e] = rresidue (r, p, k, toler, e) |
| 6964 | 268 |
| 269 ## Reconstitute the numerator and denominator polynomials from the | |
| 270 ## residues, poles, and direct term. | |
| 271 | |
| 7011 | 272 if (nargin < 2 || nargin > 5) |
| 6964 | 273 print_usage (); |
| 274 endif | |
| 275 | |
| 7011 | 276 if (nargin < 5) |
| 277 e = []; | |
| 278 endif | |
| 279 | |
| 6964 | 280 if (nargin < 4) |
| 281 toler = []; | |
| 282 endif | |
| 283 | |
| 284 if (nargin < 3) | |
| 285 k = []; | |
| 286 endif | |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
287 |
| 7011 | 288 if numel (e) |
| 289 indx = 1:numel(p); | |
| 290 else | |
| 291 [e, indx] = mpoles (p, toler, 0); | |
| 292 p = p (indx); | |
| 293 r = r (indx); | |
| 294 endif | |
| 6964 | 295 |
| 296 indx = 1:numel(p); | |
| 297 | |
| 298 for n = indx | |
| 299 pn = [1, -p(n)]; | |
| 300 if n == 1 | |
| 301 pden = pn; | |
| 302 else | |
| 303 pden = conv (pden, pn); | |
| 304 endif | |
| 305 endfor | |
| 306 | |
| 307 ## D is the order of the denominator | |
| 308 ## K is the order of the direct polynomial | |
| 309 ## N is the order of the resulting numerator | |
| 310 ## pnum(1:(N+1)) is the numerator's polynomial | |
| 311 ## pden(1:(D+1)) is the denominator's polynomial | |
| 312 ## pm is the multible pole for the nth residue | |
| 313 ## pn is the numerator contribution for the nth residue | |
| 314 | |
| 315 D = numel (pden) - 1; | |
| 316 K = numel (k) - 1; | |
| 317 N = K + D; | |
| 318 pnum = zeros (1, N+1); | |
| 7011 | 319 for n = indx(abs (r) > 0) |
| 6964 | 320 p1 = [1, -p(n)]; |
| 7011 | 321 for m = 1:e(n) |
| 322 if (m == 1) | |
| 6964 | 323 pm = p1; |
| 324 else | |
| 325 pm = conv (pm, p1); | |
| 326 endif | |
| 327 endfor | |
| 328 pn = deconv (pden, pm); | |
| 329 pn = r(n) * pn; | |
| 7183 | 330 pnum = pnum + prepad (pn, N+1, 0, 2); |
| 6964 | 331 endfor |
| 332 | |
| 333 ## Add the direct term. | |
| 334 | |
| 335 if (numel (k)) | |
| 336 pnum = pnum + conv (pden, k); | |
| 337 endif | |
| 338 | |
| 339 ## Check for leading zeros and trim the polynomial coefficients. | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
340 if (isa (r, "single") || isa (p, "single") || isa (k, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
341 small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps ("single"); |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
342 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
343 small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
344 endif |
| 6964 | 345 |
| 7011 | 346 pnum(abs (pnum) < small) = 0; |
| 347 pden(abs (pden) < small) = 0; | |
| 6964 | 348 |
| 349 pnum = polyreduce (pnum); | |
| 350 pden = polyreduce (pden); | |
| 351 | |
| 352 endfunction | |
| 6968 | 353 |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
354 |
| 6968 | 355 %!test |
| 356 %! b = [1, 1, 1]; | |
| 357 %! a = [1, -5, 8, -4]; | |
| 358 %! [r, p, k, e] = residue (b, a); | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
359 %! assert (r, [-2; 7; 3], 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
360 %! assert (p, [2; 2; 1], 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
361 %! assert (isempty (k)); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
362 %! assert (e, [1; 2; 1]); |
| 6994 | 363 %! k = [1 0]; |
| 7011 | 364 %! b = conv (k, a) + prepad (b, numel (k) + numel (a) - 1, 0); |
| 365 %! a = a; | |
| 366 %! [br, ar] = residue (r, p, k); | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
367 %! assert (br, b, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
368 %! assert (ar, a, 1e-12); |
| 7011 | 369 %! [br, ar] = residue (r, p, k, e); |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
370 %! assert (br, b, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
371 %! assert (ar, a, 1e-12); |
| 6994 | 372 |
| 373 %!test | |
| 374 %! b = [1, 0, 1]; | |
| 375 %! a = [1, 0, 18, 0, 81]; | |
| 7398 | 376 %! [r, p, k, e] = residue (b, a); |
| 6994 | 377 %! r1 = [-5i; 12; +5i; 12]/54; |
| 378 %! p1 = [+3i; +3i; -3i; -3i]; | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
379 %! assert (r, r1, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
380 %! assert (p, p1, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
381 %! assert (isempty (k)); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
382 %! assert (e, [1; 2; 1; 2]); |
| 6994 | 383 %! [br, ar] = residue (r, p, k); |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
384 %! assert (br, b, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
385 %! assert (ar, a, 1e-12); |
| 7011 | 386 |
| 387 %!test | |
| 388 %! r = [7; 3; -2]; | |
| 389 %! p = [2; 1; 2]; | |
| 390 %! k = [1 0]; | |
| 391 %! e = [2; 1; 1]; | |
| 392 %! [b, a] = residue (r, p, k, e); | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
393 %! assert (b, [1, -5, 9, -3, 1], 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
394 %! assert (a, [1, -5, 8, -4], 1e-12); |
| 7011 | 395 %! [rr, pr, kr, er] = residue (b, a); |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
396 %! [jnk, n] = mpoles (p); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
397 %! assert (rr, r(n), 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
398 %! assert (pr, p(n), 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
399 %! assert (kr, k, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
400 %! assert (er, e(n), 1e-12); |
| 7011 | 401 |
| 7188 | 402 %!test |
| 403 %! b = [1]; | |
| 404 %! a = [1, 10, 25]; | |
| 7398 | 405 %! [r, p, k, e] = residue (b, a); |
| 7188 | 406 %! r1 = [0; 1]; |
| 407 %! p1 = [-5; -5]; | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
408 %! assert (r, r1, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
409 %! assert (p, p1, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
410 %! assert (isempty (k)); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
411 %! assert (e, [1; 2]); |
| 7188 | 412 %! [br, ar] = residue (r, p, k); |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
413 %! assert (br, b, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
414 %! assert (ar, a, 1e-12); |
|
13128
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
415 |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
416 ## The following test is due to Bernard Grung (bug #34266) |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
417 %!xtest |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
418 %! z1 = 7.0372976777e6; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
419 %! p1 = -3.1415926536e9; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
420 %! p2 = -4.9964813512e8; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
421 %! r1 = -(1 + z1/p1)/(1 - p1/p2)/p2/p1; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
422 %! r2 = -(1 + z1/p2)/(1 - p2/p1)/p2/p1; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
423 %! r3 = (1 + (p2 + p1)/p2/p1*z1)/p2/p1; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
424 %! r4 = z1/p2/p1; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
425 %! r = [r1; r2; r3; r4]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
426 %! p = [p1; p2; 0; 0]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
427 %! k = []; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
428 %! e = [1; 1; 1; 2]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
429 %! b = [1, z1]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
430 %! a = [1, -(p1 + p2), p1*p2, 0, 0]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
431 %! [br, ar] = residue (r, p, k, e); |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
432 %! assert (br, b, 1e-8); |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
433 %! assert (ar, a, 1e-8); |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
434 |