View Raw SPL
/*****************************************************************************
*                                                                            *
*   TF2ZPK.SPL   Copyright (C) 2015 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Discrete transfer function to zeros, poles and gain         *
*                                                                            *
*   Revisions:   15 Jun 2015  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#if @HELP_TF2ZPK

    TF2ZPK

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

    Syntax:  TF2ZPK(b, a)

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

                 b - A series, the numerator coefficients in ascending
                     powers of z^(-1).

                 a - A series, the denominator coefficients in ascending
                     powers of z^(-1).

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

    Example:
                             z          1
             Given H(z) = ------- = ----------
                          z - 0.5   1 - 0.5z^-1

             W1: tf2zpk({1}, {1, -0.5})

             W1 contains three columns, where the first column is a zero
             at 0.0, the second column is a pole at 0.5 and the third
             column is the gain of 1.0.

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

             z == {0.0}
             p == {0.5}
             k == 1.0

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

    Example:

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

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

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

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

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

                     b(z)    b[1] + b[2] z^(-1) + ... + b[N] z^(-N+1)
             H(z) = ------ = ----------------------------------------
                     a(z)    a[1] + a[2] z^(-1) + ... + a[N] z^(-M+1)


             See TF2ZP to convert a continuous S plane transfer function to
             zeros, poles and gain.

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

    See Also:
             Residuez
             Tf2zp
             Zfreq
             Zp2tf
#endif


/* discrete transfer function to zeros, poles and gain */
tf2zpk(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(tf2zpk_iter(num, den, outargc));
}


ITERATE tf2zpk_iter(num, den, ocnt)
{
        local z, p, k, mlen, nidx, didx, idx;

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

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

        if (length(num) % 2 == 0 && length(den) == 1)
        {
                /* even fir, append a 0 term */
                num = extract(num, 1, length(num) + 1);
        }

        /* same length */
        mlen = max(length(num), length(den));

        num = extract(num, 1, mlen);
        den = extract(den, 1, mlen);

        /* remove trailing zeros where present in both coefficients */
        nidx = find(num != 0);
        didx = find(den != 0);

        nidx = (length(nidx) > 0) ? max(nidx) : 1;
        didx = max(didx);

        idx  = max(nidx, didx);

        num = extract(num, 1, idx);
        den = extract(den, 1, idx);

        /* zeros, poles, gain */
        (z, p, k) = tf2zp(num, den);

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

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