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);
}