/*****************************************************************************
*                                                                            *
*   HANKEL.SPL Copyright (C) 2014-2016 DSP Development Corporation           *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Creates a Hankel matrix                                     *
*                                                                            *
*   Revisions:    6 Jun 2014  RRR  Creation                                  *
*                 1 Jun 2016  RRR  outerprod speedup                         *
*                                                                            *
*****************************************************************************/


#if @HELP_HANKEL

    HANKEL

    Purpose: Creates a Hankel matrix                         

    Syntax:  HANKEL(c, r)

              c - A series, the first column of the output array.

              r - Optional. A series, the last row of the output array.
                  If not specified, the first element of the last row
                  equals the last column value and the remaining values
                  are all 0.

    Returns: An array, a Hankel matrix.

    Example:
             c = {1, 2, 3};
             h = hankel(c);

             h == {{1, 2, 3},
                   {2, 3, 0},
                   {3, 0, 0}}

             The first row and column of H is identical to C.

    Example:
             c = {1, 2, 3};
             r = {3, 8, 9};
             h = hankel(c, r);

             h == {{1, 2, 3},
                   {2, 3, 8},
                   {3, 8, 9}}

             The last row is identical to R and the first column is identical
             to C.

    Remarks:
             For column C and row R, a Hankel matrix takes the form:

              C[1]   C[2]   C[3]  ...  R[1]
              C[2]   C[1]   R[1]  ...  R[2]
              C[3]   R[1]   R[2]  ...  R[3]
              .      .      .   .      .
              .      .      .     .    .
              .      .      .       .  .
              C[n]   R[n-2] R[n-1] ... R[n]


             If R is not specified, a square Hankel matrix is returned
             where C represents the first row and column. The elements
             below the first anti-diagonal are zero such that:

               H[i, j] = c[i + j - 1]  if i + j - 1 <= n
               H[i, j] = 0             otherwise


             If R is specified, an unsymmetric Hankel matrix is returned
             where the first column is C and the last row is R such that:

               H[i, j] = c[i + j - 1]  if i + j - 1 <= n
               H[i, j] = r[i + j - n]  otherwise

             In this case, if C[end] != R[1], C[end] is used as the last
             element.
            
    See Also:
             Kron
             Toeplitz
             Vander
#endif


/* create a symmetric or unsymmetric Hankel matrix */
hankel(c, r)
{
        local nr, nc, hk, hs, cc;

        if (argc < 1)
        {
                error("hankel - input series required");
        }

        /* cast and flatten */
        c  = unravel({c});
        nc = length(c);

        if (nc == 0) return({});

        if (argc < 2)
        {
                r  = zeros(nc, 1);
                nr = nc;
        }
        else
        {
                /* cast and flatten row */
                r  = unravel({r});
                nr = length(r);

                if (nr == 0) return({});

                /* check diagonal conflict */
                if (r[1] != c[nc])
                {
                        echo("hankel - last column value used as last diagonal value");
                }
        }

        /* subscript vector */
        hs = concat(c, extract(r, 2, nr - 1));

        /* subscript array */
        hk = outerprod(repcol(1..nc, nr), 0..nr - 1, "+");

        /* preserve Hankel shape */
        hk[..] = hs[hk];

        return(hk);
}