/*****************************************************************************
*                                                                            *
*   CAS2ZP.SPL   Copyright (C) 2015 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Cascade form to zeros, poles and gain                       *
*                                                                            *
*   Revisions:   15 Jun 2015  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#if @HELP_CAS2ZP

    CAS2ZP

    Purpose: Converts cascade coefficients to zeros, poles and gain.

    Syntax:  CAS2ZP(c)

             (z, p, k) = CAS2ZP(c)

                 c - A series, the cascade coefficients.

    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) = cas2zp(c) returns the zeros, poles and gain as
             three separate arrays.

    Example:
             c = {1, 1, -2, 0, -0.7, 0.1};
             (z, p, k) = cas2zp(c);

             z = {0.0, 2.0};
             p = {0.5, 0.2};
             k = 1.0;

             The cascade coefficients represent the system:

                                 1 - 2z^(-1)
             H(z) = (1) * --------------------------
                          1 - 0.7 z^(-1) + 0.1 z^(-2)


             z contains the zeros, p contains the poles and k is the
             gain of the system.

    Remarks:
             CAS2ZP converts cascade coefficients to zeros, poles and gain
             of a discrete system where the input coefficients
             represent the following Z transform:

                               -1       -2                -1       -2
                    b10 + b11 z  + b12 z       b20 + b21 z  + b22 z
         H(z) = G * ________________________ *  ________________________ ...
                               -1       -2                -1       -2
                      1 + a11 z  + a12 z         1 + a21 z  + a22 z

             The cascade filter coefficients are represented by a single
             column series in the following order:

             G, b10, b11, b12, a11, a12, b20, b21, b22, a21, a22, ...

             See ZP2CAS to convert zeros, poles and gain to cascade
             coefficients.

    See Also:
             Cascade
             Tf2zpk
             Zp2cas
             Zp2tf
#endif


/* cascade coefficients to zeros, poles and gain */
ITERATE cas2zp(c)
{
        local z, p, k, pidx, zidx, mlen, dx, zpk;

        if (argc < 1)
        {
                c = {};
        }

        c = {c};

        if ((length(c) - 1) % 5 != 0)
        {
                error("cas2zp - improper cascade series");
        }

        /* gain term */
        k = c[1];

        /* sample interval */
        dx = deltax(c);

        /* remove gain from coefficients */
        c = extract(c, 2, -1);

        /* make each column a stage */
        c = ravel(c, 5);

        /* numerator coefficients */
        z = extract(c, 1, 3);

        /* find roots of each column */
        z = polyroot(z, 1);

        /* single series */
        z = unravel(z);

        /* denominator coefficients */
        p = extract(c, 4, 2);

        /* add leading 1's */
        p = concat(ones(1, numcols(p)), p);

        /* find roots of each stage */
        p = polyroot(p, 1);

        /* single series */
        p = unravel(p);

        /* remove cancelling 0's */
        zidx = find(z == 0.0);
        pidx = find(p == 0.0);

        mlen = min(length(zidx), length(pidx));

        if (mlen > 0)
        {
                zidx = extract(zidx, 1, mlen);
                pidx = extract(pidx, 1, mlen);

                z[zidx] = {};
                p[pidx] = {};
        }

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

        setdeltax(z, dx);
        setdeltax(p, dx);

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

                /* use insert to handle empty arrays */
                zpk = k;
                zpk = insertcol(zpk, p, -1);
                zpk = insertcol(zpk, z, -1);

                return(zpk);
        }
}