octave_packages/signal-1.1.3/sos2zp.m

   1 %% Copyright (C) 2005 Julius O. Smith III <jos@ccrma.stanford.edu>
   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 %% -*- texinfo -*-
  17 %% @deftypefn {Function File} {[@var{z}, @var{p}, @var{g}] =} sos2zp (@var{sos}, @var{Bscale})
  18 %% Convert series second-order sections to zeros, poles, and gains
  19 %% (pole residues).
  20 %%
  21 %% INPUTS:@*
  22 %% @itemize
  23 %%
  24 %% @item
  25 %% @var{sos} = matrix of series second-order sections, one per row:@*
  26 %% @var{sos} = [@var{B1}.' @var{A1}.'; ...; @var{BN}.' @var{AN}.'], where@*
  27 %% @code{@var{B1}.'==[b0 b1 b2] and @var{A1}.'==[1 a1 a2]} for
  28 %% section 1, etc.@*
  29 %% b0 must be nonzero for each section.
  30 %% See @code{filter()} for documentation of the
  31 %% second-order direct-form filter coefficients @var{B}i and @var{A}i.
  32 %%
  33 %% @item
  34 %% @var{Bscale} is an overall gain factor that effectively scales
  35 %% any one of the input @var{B}i vectors.
  36 %% @end itemize
  37 %%
  38 %% RETURNED:
  39 %%   @itemize
  40 %%   @item
  41 %%   @var{z} = column-vector containing all zeros (roots of B(z))@*
  42 %%   @item
  43 %%   @var{p} = column-vector containing all poles (roots of A(z))@*
  44 %%   @item
  45 %%   @var{g} = overall gain = @var{B}(Inf)
  46 %%   @end itemize
  47 %%
  48 %% EXAMPLE:
  49 %% @example
  50 %%   [z,p,g] = sos2zp([1 0 1, 1 0 -0.81; 1 0 0, 1 0 0.49])
  51 %%   => z = [i; -i; 0; 0], p = [0.9, -0.9, 0.7i, -0.7i], g=1
  52 %% @end example
  53 %%
  54 %% @seealso{zp2sos sos2tf tf2sos zp2tf tf2zp}
  55 %% @end deftypefn
  56
  57 function [z,p,g] = sos2zp (sos, Bscale = 1)
  58
  59   if (nargin < 1 || nargin > 2)
  60     print_usage;
  61   endif
  62
  63   gains = sos(:,1); % All b0 coeffs
  64   g = prod(gains)*Bscale; % pole-zero gain
  65   if g==0, error('sos2zp: one or more section gains is zero'); end
  66   sos(:,1:3) = sos(:,1:3)./ [gains gains gains];
  67
  68   [N,m] = size(sos);
  69   if m~=6, error('sos2zp: sos matrix should be N by 6'); end
  70
  71   z = zeros(2*N,1);
  72   p = zeros(2*N,1);
  73   for i=1:N
  74     ndx = [2*i-1:2*i];
  75     zi = roots(sos(i,1:3));
  76     z(ndx) = zi;
  77     pi = roots(sos(i,4:6));
  78     p(ndx) = pi;
  79   end
  80 end
  81
  82 %!test
  83 %! b1t=[1 2 3]; a1t=[1 .2 .3];
  84 %! b2t=[4 5 6]; a2t=[1 .4 .5];
  85 %! sos=[b1t a1t; b2t a2t];
  86 %! z = [-1-1.41421356237310i;-1+1.41421356237310i;...
  87 %!     -0.625-1.05326872164704i;-0.625+1.05326872164704i];
  88 %! p = [-0.2-0.678232998312527i;-0.2+0.678232998312527i;...
  89 %!      -0.1-0.538516480713450i;-0.1+0.538516480713450i];
  90 %! k = 4;
  91 %! [z2,p2,k2] = sos2zp(sos,1);
  92 %! assert({cplxpair(z2),cplxpair(p2),k2},{z,p,k},100*eps);