View Raw SPL
/*****************************************************************************
*                                                                            *
*   VANDER.SPL   Copyright (C) 2014 DSP Development Corporation              *
*                                 All Rights Reserved                        *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Creates a Vandermonde matrix                                *
*                                                                            *
*   Revisions:    8 Jun 2014  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/


#if @HELP_VANDER

    VANDER

    Purpose: Creates a Vandermonde matrix

    Syntax:  VANDER(c, n)

              c - A series, the base of the Vandermonde matrix.

              n - Optional. An integer, the degree. Defaults to length(c).

    Returns: An array, a Vandermonde matrix.

    Example:
             c = {1, 2, 3, 4};
             v = vander(c);

             v == {{ 1,  1, 1, 1},
                   { 8,  4, 2, 1},
                   {27,  9, 3, 1},
                   {64, 16, 4, 1}}

             The second to last column of V is identical to C.

    Example:
             c = {1, 2, 3, 4};
             v = vander(c, 3);

             v == {{ 1, 1, 1},
                   { 4, 2, 1},
                   { 9, 3, 1},
                   {16, 4, 1}}

             A third degree Vandermonde matrix is returned where the
             result is the same as the last 3 columns of the full
             Vandermonde matrix.

    Remarks:
             For column C, a Vandermonde matrix takes the form:

              C[1]^(n-1) ...  C[1]^3  C[1]^2   C[1]  1
              C[2]^(n-2) ...  C[2]^2  C[2]^2   C[2]  1
              C[3]^(n-3) ...  C[3]^3  C[3]^2   C[3]  1
              .            .  .       .        .     .
              .           .   .       .        .     .
              .          .    .       .        .     .
              C[n]^(n-1) ...  C[n]^3  C[n]^2   C[n]  1


             If the degree N is not specified, a square Vandermonde
             matrix is returned where the degree is set to length(C)-1.

             If N is specified, a non-square Vandermonde matrix of size
             length(C) x N is returned.
            
    See Also:
             Hankel
             Kron
             Toeplitz
#endif


/* create a Vandermonde matrix */
vander(c, n)
{
        local m;

        if (argc < 2)
        {
                if (argc < 1) error("vander - input series required");

                n = -1;
        }

        /* force series */        
        c = {c};

        if (numcols(c) > 1)
        {
                error("vander - input must be a single column series")
        }

        if (n < 0)
        {
                n = length(c);
        }

        /*
         *  Note: the quick and dirty:
         *
         *        m = repcol(c, n) ^ transpose((n - 1)..-1..0);
         *
         *        is inaccurate for high powers of c, so we compute the power
         *        operation by cumulative multiplication
         */

        /* extend columns and transpose */
        m = repcol(c, n)';

        /* set first row to all ones */
        m[1, ..] = 1;

        /* cumulative multiply each element and transpose back */
        m = interpose(m, "*")';

        /* reverse columns */
        m = fliplr(m);

        return(m);
}