Script 'mail_helper' called by obssrc Hello community, here is the log from the commit of package octave-forge-signal for openSUSE:Factory checked in at 2022-11-09 12:57:52 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Comparing /work/SRC/openSUSE:Factory/octave-forge-signal (Old) and /work/SRC/openSUSE:Factory/.octave-forge-signal.new.1597 (New) ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Package is "octave-forge-signal" Wed Nov 9 12:57:52 2022 rev:6 rq:1034677 version:1.4.3 Changes: -------- --- /work/SRC/openSUSE:Factory/octave-forge-signal/octave-forge-signal.changes 2022-07-19 17:19:17.336360804 +0200 +++ /work/SRC/openSUSE:Factory/.octave-forge-signal.new.1597/octave-forge-signal.changes 2022-11-09 12:58:18.680644196 +0100 @@ -1,0 +2,10 @@ +Tue Nov 8 17:08:13 UTC 2022 - Stefan Br??ns <stefan.bru...@rwth-aachen.de> + +- Update to version 1.4.3: + * Minor bug fixes and documentation improvements have been made to + the following functions: + fir2 remez residuez + sos2tf sos2zp xcorr + zp2sos + +------------------------------------------------------------------- Old: ---- signal-1.4.2.tar.gz New: ---- signal-1.4.3.tar.gz ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Other differences: ------------------ ++++++ octave-forge-signal.spec ++++++ --- /var/tmp/diff_new_pack.KVO9ju/_old 2022-11-09 12:58:19.140646790 +0100 +++ /var/tmp/diff_new_pack.KVO9ju/_new 2022-11-09 12:58:19.148646835 +0100 @@ -18,12 +18,12 @@ %define octpkg signal Name: octave-forge-%{octpkg} -Version: 1.4.2 +Version: 1.4.3 Release: 0 Summary: Signal processing tools for Octave License: GPL-3.0-or-later AND SUSE-Public-Domain Group: Productivity/Scientific/Math -URL: https://octave.sourceforge.io +URL: https://gnu-octave.github.io/packages/%{octpkg}/ Source0: https://downloads.sourceforge.net/octave/%{octpkg}-%{version}.tar.gz BuildRequires: gcc-c++ BuildRequires: hdf5-devel ++++++ signal-1.4.2.tar.gz -> signal-1.4.3.tar.gz ++++++ diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/signal-1.4.2/DESCRIPTION new/signal-1.4.3/DESCRIPTION --- old/signal-1.4.2/DESCRIPTION 2022-04-23 13:21:25.000000000 +0200 +++ new/signal-1.4.3/DESCRIPTION 2022-10-28 13:26:05.000000000 +0200 @@ -1,6 +1,6 @@ Name: signal -Version: 1.4.2 -Date: 2022-04-22 +Version: 1.4.3 +Date: 2022-10-28 Author: various authors Maintainer: Mike Miller <mtmil...@octave.org> Title: Signal Processing diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/signal-1.4.2/NEWS new/signal-1.4.3/NEWS --- old/signal-1.4.2/NEWS 2022-04-23 13:21:25.000000000 +0200 +++ new/signal-1.4.3/NEWS 2022-10-28 13:26:05.000000000 +0200 @@ -1,3 +1,13 @@ +Summary of important user-visible changes for signal 1.4.3: +---------------------------------------------------------- + + ** Minor bug fixes and documentation improvements have been made to the + following functions: + + fir2 remez residuez + sos2tf sos2zp xcorr + zp2sos + Summary of important user-visible changes for signal 1.4.2: ---------------------------------------------------------- diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/signal-1.4.2/inst/fir2.m new/signal-1.4.3/inst/fir2.m --- old/signal-1.4.2/inst/fir2.m 2022-04-23 13:21:25.000000000 +0200 +++ new/signal-1.4.3/inst/fir2.m 2022-10-28 13:26:05.000000000 +0200 @@ -74,6 +74,13 @@ if nargin < 4, grid_n=[]; endif if nargin < 5, ramp_n=[]; endif + if (! (isscalar (n) && (n == fix (n)) && (n >= 0))) + error ("fir2: n must be a non negative integer"); + else + # ensure n is used as a double type + n = double(n); + endif + ## find the window parameter, or default to hamming w=[]; if length(grid_n)>1, w=grid_n; grid_n=[]; endif @@ -189,6 +196,14 @@ %! assert (h(4) <= 1/sqrt (2)) %! assert (h(5), 1, 2e-3) +%!test #bug 59066 +%! f = [0, 0.45, 0.45, 0.55, 0.55, 1]; m = [1, 1, 0, 0, 1, 1]; +%! b = fir2 (int32(50), f, m); +%! assert(numel(b), 51) +%! +%! fail ("fir2 (50.1, f, m)", "fir2: n must be a non negative integer") +%! fail ("fir2 (-1, f, m)", "fir2: n must be a non negative integer") + %!demo %! f=[0, 0.3, 0.3, 0.6, 0.6, 1]; m=[0, 0, 1, 1/2, 0, 0]; %! [h, w] = freqz(fir2(100,f,m)); diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/signal-1.4.2/inst/residuez.m new/signal-1.4.3/inst/residuez.m --- old/signal-1.4.2/inst/residuez.m 2022-04-23 13:21:25.000000000 +0200 +++ new/signal-1.4.3/inst/residuez.m 2022-10-28 13:26:05.000000000 +0200 @@ -76,7 +76,7 @@ [r,p,f,m]=residue(conj(fliplr(NUM)),conj(fliplr(DEN))); p = 1 ./ p; r = r .* ((-p) .^m); - if f, f = conj(fliplr(f)); endif + f = conj(fliplr(f)); endfunction @@ -170,3 +170,10 @@ %! assert(p(is),pise,100*eps); %! assert(f,[],100*eps); %! assert(m,[1;1;1;1;1],100*eps); + +%!test # bug 57359 +%! [r,p,k] = residuez([1 1 1.5 .5],[1 1.5 .5]); +%! [rs,is] = sort(r); +%! assert(r(is), [-1; 2], 100*eps); +%! assert(p(is), [-0.5; -1], 100*eps); +%! assert(k, [0 1], 100*eps); diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/signal-1.4.2/inst/sos2tf.m new/signal-1.4.3/inst/sos2tf.m --- old/signal-1.4.2/inst/sos2tf.m 2022-04-23 13:21:25.000000000 +0200 +++ new/signal-1.4.3/inst/sos2tf.m 2022-10-28 13:26:05.000000000 +0200 @@ -1,4 +1,5 @@ ## Copyright (C) 2005 Julius O. Smith III <j...@ccrma.stanford.edu> +## Copyright (C) 2021 Charles Praplan ## ## 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 @@ -17,7 +18,7 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{b}, @var{a}] =} sos2tf (@var{sos}) ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} sos2tf (@var{sos}, @var{g}) -## Convert series second-order sections to direct form @math{H(z) = B(z)/A(z)}. +## Convert series second-order sections to transfer function. ## ## INPUTS: ## @itemize @@ -28,8 +29,13 @@ ## @var{sos} = [@var{B1}.' @var{A1}.'; ...; @var{BN}.' @var{AN}.'] ## @end example ## where -## @code{@var{B1}.' = [b0 b1 b2] and @var{A1}.' = [1 a1 a2]} for -## section 1, etc. The b0 entry must be nonzero for each section. +## @code{@var{B1}.' = [b0 b1 b2] and @var{A1}.' = [a0 a1 a2]} for +## section 1, etc. +## +## a0 is usually equal to 1 because all 2nd order transfer functions +## can be scaled so that a0 = 1. +## However, this is not mandatory for this implementation, which supports +## all kinds of transfer functions, including first order transfer functions. ## See @code{filter} for documentation of the second-order direct-form filter ## coefficients @var{B}i and @var{A}i. ## @@ -40,7 +46,8 @@ ## @end itemize ## ## RETURNED: -## @var{b} and @var{a} are vectors specifying the digital filter @math{H(z) = B(z)/A(z)}. +## @var{b} and @var{a} are vectors specifying the analog or digital filter +## @math{H(s) = B(s)/A(s)} or @math{H(z) = B(z)/A(z)}. ## See @code{filter} for further details. ## ## @seealso{tf2sos, zp2sos, sos2pz, zp2tf, tf2zp} @@ -66,16 +73,26 @@ A = conv(A, sos(i,4:6)); endfor - nB = length(B); - while nB && B(nB)==0 - B=B(1:nB-1); - nB=length(B); + nA = length (A); + nB = length (B); + if nA ~= nB + error('Internal error: length (A) not equal to length (B)'); + # This error cannot occur (if the above calls to the conv function + # are not modified and if the conv function is also not modified) + endif; + + # Removing trailing zeros if present in numerator and denominator + while nB && B(nB)==0 && A(nB)==0 + B = B(1:nB-1); + A = A(1:nA-1); + nB = length (B); endwhile - nA = length(A); - while nA && A(nA)==0 - A=A(1:nA-1); - nA=length(A); + # Removing leading zeros if present in numerator and denominator + while nB && B(1)==0 && A(1)==0 + A = A(2:end); + B = B(2:end); + nB--; endwhile B = B .* prod (g); @@ -120,3 +137,44 @@ %! assert (Bh, 8 * B, 10*eps); %! assert (Ah, A, 10*eps); +## Test with trailing zero in numerator +%!test +%! sos = [1, 1, 0, 0, 1, 0.5]; +%! [Bh, Ah] = sos2tf (sos); +%! assert (Bh, sos(1,1:3) , 10*eps); +%! assert (Ah, sos(1,4:6), 10*eps); + +## Test with trailing zero in denominator +%!test +%! sos = [0, 1, 1, 1, 0.5, 0]; +%! [Bh, Ah] = sos2tf (sos); +%! assert (Bh, sos(1,1:3) , 10*eps); +%! assert (Ah, sos(1,4:6), 10*eps); + +## Test with trailing zero both in numerator and in denominator +%!test +%! sos = [1, 1, 0, 1, 0.5, 0]; +%! [Bh, Ah] = sos2tf (sos); +%! assert (Bh, [1, 1] , 10*eps); +%! assert (Ah, [1, 0.5], 10*eps); + +## Test with leading zero in numerator +%!test +%! sos = [0, 1, 1, 1, 1, 0.5]; +%! [Bh, Ah] = sos2tf (sos); +%! assert (Bh, sos(1,1:3) , 10*eps); +%! assert (Ah, sos(1,4:6), 10*eps); + +## Test with leading zero in denominator +%!test +%! sos = [1, 1, 0, 0, 1, 0.5]; +%! [Bh, Ah] = sos2tf (sos); +%! assert (Bh, sos(1,1:3) , 10*eps); +%! assert (Ah, sos(1,4:6), 10*eps); + +## Test with leading zero both in numerator and in denominator +%!test +%! sos = [0, 1, 1, 0, 1, 0.5]; +%! [Bh, Ah] = sos2tf (sos); +%! assert (Bh, [1, 1] , 10*eps); +%! assert (Ah, [1, 0.5], 10*eps); diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/signal-1.4.2/inst/sos2zp.m new/signal-1.4.3/inst/sos2zp.m --- old/signal-1.4.2/inst/sos2zp.m 2022-04-23 13:21:25.000000000 +0200 +++ new/signal-1.4.3/inst/sos2zp.m 2022-10-28 13:26:05.000000000 +0200 @@ -1,4 +1,5 @@ ## Copyright (C) 2005 Julius O. Smith III <j...@ccrma.stanford.edu> +## Copyright (C) 2021 Charles Praplan ## ## 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 @@ -29,8 +30,13 @@ ## @var{sos} = [@var{B1}.' @var{A1}.'; ...; @var{BN}.' @var{AN}.'] ## @end example ## where -## @code{@var{B1}.' = [b0 b1 b2] and @var{A1}.' = [1 a1 a2]} for -## section 1, etc. The b0 entry must be nonzero for each section. +## @code{@var{B1}.' = [b0 b1 b2] and @var{A1}.' = [a0 a1 a2]} for +## section 1, etc. +## +## a0 is usually equal to 1 because all 2nd order transfer functions can +## be scaled so that a0 = 1. +## However, this is not mandatory for this implementation, which supports +## all kinds of transfer functions, including first order transfer functions. ## See @code{filter} for documentation of the second-order direct-form filter ## coefficients @var{B}i and @var{A}i. ## @@ -75,22 +81,43 @@ print_usage; endif - gains = sos(:,1); # All b0 coeffs - k = prod(gains)*g; # pole-zero gain - if k==0, error('sos2zp: one or more section gains is zero'); endif - sos(:,1:3) = sos(:,1:3)./ [gains gains gains]; - [N,m] = size(sos); if m~=6, error('sos2zp: sos matrix should be N by 6'); endif - z = zeros(2*N,1); - p = zeros(2*N,1); + k = g; + for i = 1:N + b = sos(i,1:3); + N2 = 3; + while N2 && b(1)==0 # Removing leading zeros + b = b(2:end); + N2 = N2 - 1; + endwhile + + a = sos(i,4:6); + N2 = 3; + while N2 && a(1)==0 # Removing leading zeros + a = a(2:end); + N2--; + endwhile + + if isempty (a) + warning ('Infinite gain detected'); + k = Inf; + elseif a(1) == 0 + warning ('Infinite gain detected'); + k = Inf; + elseif isempty (b) + k = 0; + else + k = k * b(1) / a(1); + endif + endfor + + z = []; + p = []; for i=1:N - ndx = [2*i-1:2*i]; - zi = roots(sos(i,1:3)); - z(ndx) = zi; - pi = roots(sos(i,4:6)); - p(ndx) = pi; + z=[z; roots(sos(i,1:3))]; + p=[p; roots(sos(i,4:6))]; endfor endfunction @@ -106,3 +133,45 @@ %! k = 4; %! [z2,p2,k2] = sos2zp(sos,1); %! assert({cplxpair(z2),cplxpair(p2),k2},{z,p,k},100*eps); + +## Test with trailing zero in numerator +%!test +%! sos = [1, 1, 0, 1, 1, 0.5]; +%! [Z, P] = sos2zp (sos); +%! assert (Z, roots (sos(1,1:3)), 10*eps); +%! assert (P, roots (sos(1,4:6)), 10*eps); + +## Test with trailing zero in denominator +%!test +%! sos = [0, 1, 1, 1, 0.5, 0]; +%! [Z, P] = sos2zp (sos); +%! assert (Z, roots (sos(1,1:3)), 10*eps); +%! assert (P, roots (sos(1,4:6)), 10*eps); + +## Test with trailing zero both in numerator and in denominator +%!test +%! sos = [1, 1, 0, 1, 0.5, 0]; +%! [Z, P] = sos2zp (sos); +%! assert (Z, roots (sos(1,1:3)), 10*eps); +%! assert (P, roots (sos(1,4:6)), 10*eps); + +## Test with leading zero in numerator +%!test +%! sos = [0, 1, 1, 1, 1, 0.5]; +%! [Z, P] = sos2zp (sos); +%! assert (Z, roots (sos(1,1:3)), 10*eps); +%! assert (P, roots (sos(1,4:6)), 10*eps); + +## Test with leading zero in denominator +%!test +%! sos = [1, 1, 0, 0, 1, 0.5]; +%! [Z, P] = sos2zp (sos); +%! assert (Z, roots (sos(1,1:3)), 10*eps); +%! assert (P, roots (sos(1,4:6)), 10*eps); + +## Test with leading zero both in numerator and in denominator +%!test +%! sos = [0, 1, 1, 0, 1, 0.5]; +%! [Z, P] = sos2zp (sos); +%! assert (Z, roots (sos(1,1:3)), 10*eps); +%! assert (P, roots (sos(1,4:6)), 10*eps); diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/signal-1.4.2/inst/xcorr.m new/signal-1.4.3/inst/xcorr.m --- old/signal-1.4.2/inst/xcorr.m 2022-04-23 13:21:25.000000000 +0200 +++ new/signal-1.4.3/inst/xcorr.m 2022-10-28 13:26:05.000000000 +0200 @@ -1,6 +1,7 @@ ## Copyright (C) 1999-2001 Paul Kienzle <pkien...@users.sf.net> ## Copyright (C) 2004 <asbjorn.s...@broadpark.no> ## Copyright (C) 2008, 2010 Peter Lanspeary <peter.lanspe...@.adelaide.edu.au> +## Copyright (C) 2022 Octave-Forge community <maintain...@octave.org> ## ## 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 @@ -77,7 +78,7 @@ ## return the biased average, R/N, ## @item unbiased ## return the unbiased average, R(k)/(N-|k|), -## @item coeff +## @item coeff or normalized ## return the correlation coefficient, R/(rms(x).rms(y)), ## where "k" is the lag, and "N" is the length of @var{X}. ## If omitted, the default value is "none". @@ -135,26 +136,36 @@ ## If length(x) == length(y) + k, then you can use the simpler ## ( hankel(x(1:k),x(k:N-k)) * y ) ./ N -function [R, lags] = xcorr (X, Y, maxlag, scale) +function [R, lags] = xcorr (X, varargin) if (nargin < 1 || nargin > 4) print_usage; endif - ## assign arguments that are missing from the list - ## or reassign (right shift) them according to data type - if nargin==1 - Y=[]; maxlag=[]; scale=[]; - elseif nargin==2 - maxlag=[]; scale=[]; - if ischar(Y), scale=Y; Y=[]; - elseif isscalar(Y), maxlag=Y; Y=[]; + ## assign optional arguments according to data type + maxlag=[]; scale=[]; Y=[]; + + for idx=1:length(varargin) + arg = varargin{idx}; + if ischar(arg) + if isempty(scale) + scale = arg; + else + error ("xcorr: unexpected char value '%s' when scale is already set to '%s'", arg, scale); + endif + elseif isscalar(arg) + if isempty(maxlag) + maxlag = arg; + else + error ("xcorr: unexpected scalar value '%f' when maxlag is already set to '%f'", arg, maxlag); + endif + elseif idx == 1 + Y = arg; + else + error ("xcorr: unknown optional input variable at position %d", idx); endif - elseif nargin==3 - scale=[]; - if ischar(maxlag), scale=maxlag; maxlag=[]; endif - if isscalar(Y), maxlag=Y; Y=[]; endif - endif + + endfor ## assign defaults to missing arguments if isvector(X) @@ -167,7 +178,7 @@ if isempty(scale), scale='none'; endif ## check argument values - if isempty(X) || isscalar(X) || ischar(Y) || ! ismatrix(X) + if isempty(X) || isscalar(X) || ! ismatrix(X) error("xcorr: X must be a vector or matrix"); endif if isscalar(Y) || ischar(Y) || (!isempty(Y) && !isvector(Y)) @@ -179,6 +190,7 @@ if !isscalar(maxlag) || !isreal(maxlag) || maxlag<0 || fix(maxlag)!=maxlag error("xcorr: maxlag must be a single non-negative integer"); endif + ## ## sanity check on number of requested lags ## Correlations for lags in excess of +/-(N-1) @@ -248,7 +260,7 @@ R = R ./ N; elseif strcmp(scale, 'unbiased') R = R ./ ( [ N-maxlag:N-1, N, N-1:-1:N-maxlag ]' * ones(1,columns(R)) ); - elseif strcmp(scale, 'coeff') + elseif strcmp(scale, 'coeff') || strcmp(scale, 'normalized') ## R = R ./ R(maxlag+1) works only for autocorrelation ## For cross correlation coeff, divide by rms(X)*rms(Y). if !isvector(X) @@ -328,3 +340,91 @@ ## endfor ##endif ##-------------------------------------------------------------- + +%!shared x, y +%! x = 0.5.^(0:15); +%! y = circshift(x,5); + +## Test input validation +%!error xcorr () +%!error xcorr (1) +%!error xcorr (x, 1, x) +%!error xcorr (x, 'none', x) +%!error xcorr (x, x, 'invalid') +%!error xcorr (x, 'invalid') + +%!test +%! [c,lags] = xcorr(x); +%! # largest spike at 0 lag, where X matches itself - ie the center +%! [m, im] = max(c); +%! assert(m, 4/3, 1e-6) +%! assert(im, (numel(lags)+1)/2); +%! +%! [c1,lags1] = xcorr(x, x); +%! [m, im] = max(c1); +%! assert(m, 4/3, 1e-6) +%! assert(im, (numel(lags1)+1)/2); +%! assert(c1, c, 2*eps); +%! assert(lags1, lags); + +%!test +%! [c,lags] = xcorr(x,y); +%! # largest spike at 0 lag, where X matches Y +%! [m, im] = max(c); +%! assert(m, 4/3, 1e-6) +%! assert(lags(im), -5); + +%!test +%! [c0,lags0] = xcorr(x,y); +%! [c1,lags1] = xcorr(x,y, 'none'); +%! assert(c0, c1); +%! assert(lags0, lags1); + +%!test +%! [c0,lags0] = xcorr(x,y); +%! [c1,lags1] = xcorr(x,y, 'normalized'); +%! assert(lags0, lags1); +%! [m, im] = max(c1); +%! # at 0 lag, should be 1 +%! assert(m, 1, 1e-6); +%! [c2,lags2] = xcorr(x,y, 'coeff'); +%! assert(c1, c2); +%! assert(lags1, lags2); + +%!test +%! [c0,lags0] = xcorr(x,y); +%! [c1,lags1] = xcorr(x,y, 'biased'); +%! assert(lags0, lags1); +%! [m, im] = max(c1); +%! assert(m, 1/12, 1e-6); +%! +%! [c1,lags1] = xcorr(x, 'biased'); +%! assert(lags0, lags1); +%! [m, im] = max(c1); +%! assert(m, 1/12, 1e-6); + +%!test +%! [c0,lags0] = xcorr(x,y); +%! [c1,lags1] = xcorr(x,y, 'unbiased'); +%! assert(lags0, lags1); +%! [m, im] = max(c1); +%! assert(m, 1/8.25, 1e-6); + +%!test +%! [c,lags] = xcorr(x,y, 10); +%! [m, im] = max(c); +%! assert(lags(im), -5); +%! assert(lags(1), -10); +%! assert(lags(end), 10); +%! +%! [c,lags] = xcorr(x,10); +%! [m, im] = max(c); +%! assert(lags(1), -10); +%! assert(lags(end), 10); + +%!test +%! [c0,lags0] = xcorr(x,y, 'normalized', 10); +%! [c1,lags1] = xcorr(x,y, 10, 'normalized'); +%! assert(c0, c1); +%! assert(lags0, lags1); + diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/signal-1.4.2/inst/zp2sos.m new/signal-1.4.3/inst/zp2sos.m --- old/signal-1.4.2/inst/zp2sos.m 2022-04-23 13:21:25.000000000 +0200 +++ new/signal-1.4.3/inst/zp2sos.m 2022-10-28 13:26:05.000000000 +0200 @@ -1,154 +1,213 @@ -## Copyright (C) 2005 Julius O. Smith III <j...@ccrma.stanford.edu> -## -## 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; see the file COPYING. If not, see -## <https://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{sos}, @var{g}] =} zp2sos (@var{z}) -## @deftypefnx {Function File} {[@var{sos}, @var{g}] =} zp2sos (@var{z}, @var{p}) -## @deftypefnx {Function File} {[@var{sos}, @var{g}] =} zp2sos (@var{z}, @var{p}, @var{k}) -## @deftypefnx {Function File} {@var{sos} =} zp2sos (@dots{}) -## Convert filter poles and zeros to second-order sections. -## -## INPUTS: -## @itemize -## @item -## @var{z} = column-vector containing the filter zeros -## @item -## @var{p} = column-vector containing the filter poles -## @item -## @var{k} = overall filter gain factor -## If not given the gain is assumed to be 1. -## @end itemize -## -## RETURNED: -## @itemize -## @item -## @var{sos} = matrix of series second-order sections, one per row: -## @example -## @var{sos} = [@var{B1}.' @var{A1}.'; ...; @var{BN}.' @var{AN}.'] -## @end example -## where -## @code{@var{B1}.' = [b0 b1 b2] and @var{A1}.' = [1 a1 a2]} for -## section 1, etc. The b0 entry must be nonzero for each section. -## See @code{filter} for documentation of the second-order direct-form filter -## coefficients @var{B}i and %@var{A}i, i=1:N. -## -## @item -## @var{g} is the overall gain factor that effectively scales -## any one of the @var{B}i vectors. -## @end itemize -## -## If called with only one output argument, the overall filter gain is -## applied to the first second-order section in the matrix @var{sos}. -## -## EXAMPLE: -## @example -## [z, p, k] = tf2zp ([1 0 0 0 0 1], [1 0 0 0 0 .9]); -## [sos, g] = zp2sos (z, p, k) -## -## sos = -## 1.0000 0.6180 1.0000 1.0000 0.6051 0.9587 -## 1.0000 -1.6180 1.0000 1.0000 -1.5843 0.9587 -## 1.0000 1.0000 0 1.0000 0.9791 0 -## -## g = -## 1 -## @end example -## -## @seealso{sos2pz, sos2tf, tf2sos, zp2tf, tf2zp} -## @end deftypefn - -function [sos,g] = zp2sos(z,p,k) - - if nargin<3, k=1; endif - if nargin<2, p=[]; endif - - [zc,zr] = cplxreal(z(:)); - [pc,pr] = cplxreal(p(:)); - - ## zc,zr,pc,pr - - nzc=length(zc); - npc=length(pc); - - nzr=length(zr); - npr=length(pr); - - ## Pair up real zeros: - if nzr - if mod(nzr,2)==1, zr=[zr;0]; nzr=nzr+1; endif - nzrsec = nzr/2; - zrms = -zr(1:2:nzr-1)-zr(2:2:nzr); - zrp = zr(1:2:nzr-1).*zr(2:2:nzr); - else - nzrsec = 0; - endif - - ## Pair up real poles: - if npr - if mod(npr,2)==1, pr=[pr;0]; npr=npr+1; endif - nprsec = npr/2; - prms = -pr(1:2:npr-1)-pr(2:2:npr); - prp = pr(1:2:npr-1).*pr(2:2:npr); - else - nprsec = 0; - endif - - nsecs = max(nzc+nzrsec,npc+nprsec); - - ## Convert complex zeros and poles to real 2nd-order section form: - zcm2r = -2*real(zc); - zca2 = abs(zc).^2; - pcm2r = -2*real(pc); - pca2 = abs(pc).^2; - - sos = zeros(nsecs,6); - sos(:,1) = ones(nsecs,1); # all 2nd-order polynomials are monic - sos(:,4) = ones(nsecs,1); - - nzrl=nzc+nzrsec; # index of last real zero section - nprl=npc+nprsec; # index of last real pole section - - for i=1:nsecs - - if i<=nzc # lay down a complex zero pair: - sos(i,2:3) = [zcm2r(i) zca2(i)]; - elseif i<=nzrl # lay down a pair of real zeros: - sos(i,2:3) = [zrms(i-nzc) zrp(i-nzc)]; - endif - - if i<=npc # lay down a complex pole pair: - sos(i,5:6) = [pcm2r(i) pca2(i)]; - elseif i<=nprl # lay down a pair of real poles: - sos(i,5:6) = [prms(i-npc) prp(i-npc)]; - endif - endfor - - ## If no output argument for the overall gain, combine it into the - ## first section. - if (nargout < 2) - sos(1,1:3) *= k; - else - g = k; - endif - -endfunction - -%!test -%! B=[1 0 0 0 0 1]; A=[1 0 0 0 0 .9]; -%! [z,p,k] = tf2zp(B,A); -%! [sos,g] = zp2sos(z,p,k); -%! [Bh,Ah] = sos2tf(sos,g); -%! assert({Bh,Ah},{B,A},100*eps); +## Copyright (C) 2005 Julius O. Smith III <j...@ccrma.stanford.edu> +## Copyright (C) 2021 Charles Praplan +## +## 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; see the file COPYING. If not, see +## <https://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{sos}, @var{g}] =} zp2sos (@var{z}) +## @deftypefnx {Function File} {[@var{sos}, @var{g}] =} zp2sos (@var{z}, @var{p}) +## @deftypefnx {Function File} {[@var{sos}, @var{g}] =} zp2sos (@var{z}, @var{p}, @var{k}) +## @deftypefnx {Function File} {@var{sos} =} zp2sos (@dots{}) +## Convert filter poles and zeros to second-order sections. +## +## INPUTS: +## @itemize +## @item +## @var{z} = column-vector containing the filter zeros +## @item +## @var{p} = column-vector containing the filter poles +## @item +## @var{k} = overall filter gain factor. If not given the gain is assumed to be 1. +## @end itemize +## +## RETURNED: +## @itemize +## @item +## @var{sos} = matrix of series second-order sections, one per row: +## @example +## @var{sos} = [@var{B1}.' @var{A1}.'; ...; @var{BN}.' @var{AN}.'] +## @end example +## where +## @code{@var{B1}.' = [b0 b1 b2] and @var{A1}.' = [a0 a1 a2]} for +## section 1, etc. +## See @code{filter} for documentation of the second-order direct-form filter +## coefficients @var{B}i and %@var{A}i, i=1:N. +## +## @item +## @var{g} is the overall gain factor that effectively scales +## any one of the @var{B}i vectors. +## @end itemize +## +## If called with only one output argument, the overall filter gain is +## applied to the first second-order section in the matrix @var{sos}. +## +## EXAMPLE: +## @example +## [z, p, k] = tf2zp ([1 0 0 0 0 1], [1 0 0 0 0 .9]); +## [sos, g] = zp2sos (z, p, k) +## +## sos = +## 1.0000 0.6180 1.0000 1.0000 0.6051 0.9587 +## 1.0000 -1.6180 1.0000 1.0000 -1.5843 0.9587 +## 0 1.0000 1.0000 0 1.0000 0.9791 +## +## g = +## 1 +## @end example +## +## @seealso{sos2zp, sos2tf, tf2sos, zp2tf, tf2zp} +## @end deftypefn + +function [SOS,G] = zp2sos(z,p,k,DoNotCombineReal) + + if nargin<3, k=1; endif + if nargin<2, p=[]; endif + + DoNotCombineReal = 1; + + [zc, zr] = cplxreal (z(:)); + [pc, pr] = cplxreal (p(:)); + + nzc = length (zc); + npc = length (pc); + + nzr = length (zr); + npr = length (pr); + + if DoNotCombineReal + + # Handling complex conjugate poles + for count = 1:npc + SOS(count, 4:6) = [1, -2 * real(pc(count)), abs(pc(count))^2]; + endfor + + # Handling real poles + for count = 1:npr + SOS(count + npc, 4:6) = [0, 1, -pr(count)]; + endfor + + # Handling complex conjugate zeros + for count = 1:nzc + SOS(count, 1:3) = [1, -2 * real(zc(count)), abs(zc(count))^2]; + endfor + + # Handling real zeros + for count = 1:nzr + SOS(count+nzc, 1:3) = [0, 1, -zr(count)]; + endfor + + # Completing SOS if needed (sections without pole or zero) + if npc + npr > nzc + nzr + for count = nzc + nzr + 1 : npc + npr % sections without zero + SOS(count, 1:3) = [0, 0, 1]; + end + else + for count = npc + npr + 1 : nzc + nzr % sections without pole + SOS(count, 4:6) = [0, 0, 1]; + endfor + endif + + else + + # Handling complex conjugate poles + for count = 1:npc + SOS(count, 4:6) = [1, -2 * real(pc(count)), abs(pc(count))^2]; + endfor + + # Handling pair of real poles + for count = 1:floor(npr/2) + SOS(count+npc, 4:6) = [1, - pr(2 * count - 1) - pr(2 * count), pr(2 * count - 1) * pr(2 * count)]; + endfor + + # Handling last real pole (if any) + if mod (npr,2) == 1 + SOS(npc + floor (npr / 2) + 1, 4:6)= [0, 1, -pr(end)]; + endif + + + # Handling complex conjugate zeros + for count = 1:nzc + SOS(count, 1:3)= [1, -2 * real(zc(count)), abs(zc(count))^2]; + endfor + + # Handling pair of real zeros + for count = 1:floor(nzr / 2) + SOS(count+nzc, 1:3)= [1, - zr(2 * count - 1) - zr(2 * count), zr(2 * count - 1) * zr(2 * count)]; + endfor + + # Handling last real zero (if any) + if mod (nzr, 2) == 1 + SOS(nzc + floor (nzr / 2) + 1, 1:3) = [0, 1, -zr(end)]; + endif + + # Completing SOS if needed (sections without pole or zero) + if npc + ceil(npr / 2) > nzc + ceil(nzr / 2) + for count = nzc + ceil (nzr / 2) + 1 : npc + ceil (npr / 2) % sections without zero + SOS(count, 1:3) = [0, 0, 1]; + endfor + else + for count = npc + ceil(npr / 2) + 1:nzc + ceil (nzr / 2) % sections without pole + SOS(count, 4:6) = [0, 0, 1]; + endfor + endif + endif + + if ~exist ('SOS') + SOS=[0, 0, 1, 0, 0, 1]; + endif + + ## If no output argument for the overall gain, combine it into the + ## first section. + if (nargout < 2) + SOS(1,1:3) = k * SOS(1,1:3); + else + G = k; + endif + +%!test +%! B=[1 0 0 0 0 1]; A=[1 0 0 0 0 .9]; +%! [z,p,k] = tf2zp(B,A); +%! [sos,g] = zp2sos(z,p,k); +%! [Bh,Ah] = sos2tf(sos,g); +%! assert({Bh,Ah},{B,A},100*eps); + +%!test +%! sos = zp2sos ([]); +%! assert (sos, [0, 0, 1, 0, 0, 1], 100*eps); + +%!test +%! sos = zp2sos ([], []); +%! assert (sos, [0, 0, 1, 0, 0, 1], 100*eps); + +%!test +%! sos = zp2sos ([], [], 2); +%! assert (sos, [0, 0, 2, 0, 0, 1], 100*eps); + +%!test +%! [sos, g] = zp2sos ([], [], 2); +%! assert (sos, [0, 0, 1, 0, 0, 1], 100*eps); +%! assert (g, 2, 100*eps); + +%!test +%! sos = zp2sos([], [0], 1); +%! assert (sos, [0, 0, 1, 0, 1, 0], 100*eps); + +%!test +%! sos = zp2sos([0], [], 1); +%! assert (sos, [0, 1, 0, 0, 0, 1], 100*eps); + +%!test +%! sos = zp2sos([-1-j -1+j], [-1-2j -1+2j], 10); +%! assert (sos, [10, 20, 20, 1, 2, 5], 100*eps); diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/signal-1.4.2/src/remez.cc new/signal-1.4.3/src/remez.cc --- old/signal-1.4.2/src/remez.cc 2022-04-23 13:21:25.000000000 +0200 +++ new/signal-1.4.3/src/remez.cc 2022-10-28 13:26:05.000000000 +0200 @@ -763,7 +763,7 @@ Parks-McClellan optimal FIR filter design.\n\ @table @var\n\ @item n\n\ -gives the number of taps in the returned filter\n\ +gives the filter order, where the generated filter length taps is n+1\n\ @item f\n\ gives frequency at the band edges [b1 e1 b2 e2 b3 e3 @dots{}]\n\ @item a\n\