R/polyprod.R

In fda: Functional Data Analysis

#  --------------------------------------------------------------

polyprod <- function(Coeff1, Coeff2){
# POLYCONV computes products of polynomials defined by columns of 
#   coefficient matrices Coeff1 and Coeff2

#  Last modified 6 June 2005

polyorder1 <- dim(Coeff1)[1] 
norder1    <- dim(Coeff1)[2] 
polyorder2 <- dim(Coeff2)[1] 
norder2    <- dim(Coeff2)[2] 
ndegree1 <- polyorder1 - 1 
ndegree2 <- polyorder2 - 1 

#  if the degrees are not equal, pad out the smaller matrix with 0s

if(ndegree1 != ndegree2){
    if(ndegree1 > ndegree2)  Coeff2 <- rbind(Coeff2, matrix(0,ndegree1-ndegree2,norder2))
    else                     Coeff1 <- rbind(Coeff1, matrix(0,ndegree2-ndegree1,norder1))
}

#  find order of the product
D <- max(c(ndegree1,ndegree2))   # maximum degree
N <- 2*D+1                       # order of product

#  compute the coefficients for the products
convmat <- array(0,c(norder1,norder2,N)) 
for (i in 0:(D-1)){
    ind <- c(0:i) + 1 
    if (length(ind) == 1) {
        convmat[,,i+1] <-     outer(Coeff1[ind,    ],Coeff2[i-ind+2,]) 
        convmat[,,N-i] <-     outer(Coeff1[D-ind+2,],Coeff2[D-i+ind,])	
    } else {
        convmat[,,i+1] <- crossprod(Coeff1[ind,    ],Coeff2[i-ind+2,]) 
        convmat[,,N-i] <- crossprod(Coeff1[D-ind+2,],Coeff2[D-i+ind,])
    }
}
ind <- c(0:D) + 1 
convmat[,,D+1] <-     crossprod(Coeff1[ind,    ],Coeff2[D-ind+2,])

if (ndegree1 != ndegree2) {
	convmat <- convmat[,,1:(ndegree1+ndegree2+1)] 
	convmat <- array(convmat,c(norder1,norder2,ndegree1+ndegree2+1))
}

return(convmat)
}