In R I can calculate a linear convolution of a vector y with a convolution kernel K using
linconv = function (y,K) { K=K/sum(K)
out=convolve(y, K, type="open")
out=head(out,-(length(out)-length(y)) )
return(out)
}
E.g. if I had a vector y
t = 1:1000
y = 10000*exp(sin(t*0.02))
and convolution kernel K
K = dexp(1:41, rate=1/5)
the linear convolution of y with convolution kernel K would be linconv(y, K)
Now I would like to calculate my linear convolution instead using matrix multiplication, by setting up a banded convolution matrix containing as columns time-shifted copies of my convolution kernel K. Ideally I would like to use Matrix::bandSparse to set up this banded convolution matrix. Would anyone happen to know the right syntax of bandSparse though to accomplish this, and obtain the same result as using the linear convolution function above, but obtained using matrix multiplication instead? Right now I did
linconv2 = function (y,K) { K=K/sum(K)
require(Hmisc)
lags = 0:(N-1)
X = as.matrix(data.frame(do.call(cbind, lapply(lags, function(lag) Lag(c(K, rep(0, length(y))), shift=lag)))))
X[is.na(X)] = 0
ntoclip = (nrow(X)-length(y))/2
X = X[head(1:nrow(X), -(nrow(X)-length(y))),]
out = as.vector(y %*% X)
return(out)
}
but I would like to set up the matrix X above instead using Matrix::bandSparse. Would anyone know the correct syntax to do that?
Ha found the correct syntax after all:
linconv2 = function (y,K) { K=K/sum(K)
require(Matrix)
X = as.matrix(bandSparse(length(y),
k = seq(-(length(K)-1),0,1),
diag = t(replicate(length(y), rev(K))), symm=FALSE))
out = X %*% as.matrix(y, ncol=1)
return(out)
}