/*****************************************************************************
*                                                                            *
*   TF2ZP.SPL    Copyright (C) 2015 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Converts transfer function form to zeros, poles gain        *
*                                                                            *
*   Revisions:   18 Jun 2015  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#if @HELP_TF2ZP

    TF2ZP

    Purpose: Converts S plane transfer function form to zeros, poles and gain.

    Syntax:  TF2ZP(b, a)

             (z, p, k) = TF2ZP(b, a)

                 b - A series, the numerator coefficients in descending
                     powers of s.

                 a - A series, the denominator coefficients in descending
                     powers of s.

    Returns: A Nx3 array where the first column contains the zeros, the second
             column contains the poles and the third column contains the gain.
             
             (z, p, k) = tf2zp(b, a) returns the zeros, poles and gain as
             three separate arrays.

    Example:
                             s
             Given H(s) = -------
                          s + 0.5

             W1: tf2zp({1}, {1, 0.5})

             W1 contains two columns, where the first column is a pole at
             s = -0.5 and the second column is the gain of 1.0. The system
             has no zeros.

    Example:
             (z, p, k) = tf2zp({1}, {1, 0.5})

             z == {}
             p == {-0.5}
             k == 1

             Same as above except the values are returned in three separate
             variables.

    Example:

                              s^2 + 2s
             Given H(z) = ----------------
                          s^2 + 0.7s + 0.1

             b = {1, 2, 0};
             a = {1, 0.7, 0.1};
             (z, p, k) = tf2zp(b, a);

             z == {0, -2.0)
             p == {-0.5, -0.2}
             k == 1

    Remarks:
             For TF2ZP(b, a), the input series represent the
             terms of the rational polynomial H(s) = b(s)/a(s) where:

             N = length(b) and M = length(a):

                     b(s)    b[1] s^(N-1) + b[1] s^(N-1} + ... + b[N]
             H(s) = ------ = ----------------------------------------
                     a(s)    a[1] s^(M-1) + a[1] s^(M-1} + ... + a[M]


             See TF2ZPK to convert a discrete Z plane transfer function
             to zeros, poles and gain.

             See ZP2TF to convert zeros, poles and gain to transfer function
             form.

    See Also:
             Residue
             Sfreq
             Splane
             Tf2zpk
             Zp2tf
             Zplane
#endif


/* transfer function (direct form) to zeros, poles and gain */
tf2zp(num, den)
{
        if (argc < 2)
        {
                den = {};

                if (argc == 1)
                {
                        if (numcols(num) == 2)
                        {
                                /* assume raveled b, a */
                                den = col(num, 2);
                                num = col(num, 1);
                        }
                        else
                        {
                                /* assume FIR form */
                                den = ones(1, numcols(num));
                        }
                }
        }

        return(tf2zp_iter(num, den, outargc));
}


ITERATE tf2zp_iter(num, den, ocnt)
{
        local nz, z, p, k;

        if (isempty(den))
        {
                error("tf2zp - input denominator cannot be empty");
        }

        if (den[1] == 0)
        {
                error("tf2zp - leading denominator term cannot be zero");
        }

        /* strip leading numerator zeros */
        nz  = find(not(num == 0));

        if (length(nz) > 0 && not(isempty(num)))
        {
                num = num[nz[1]..end];
                k   = num[1] / den[1];
        }
        else
        {
                num = {};
                k   = 0;
        }
                
        z = roots(num);
        p = roots(den);

        setcomment(z, "Zeros");
        setcomment(p, "Poles");

        if (ocnt > 0)
        {
                return(z, p, k);
        }
        else
        {
                k = {k};
                setcomment(k, "Gain");

                return(ravel(z, p, {k}));
        }
}