R 中使用 bandSparse 和矩阵乘法的线性卷积
Linear convolution in R using bandSparse & matrix multiplication
在 R 中,我可以使用
计算向量 y
与卷积核 K
的线性卷积
linconv = function (y,K) { K=K/sum(K)
out=convolve(y, K, type="open")
out=head(out,-(length(out)-length(y)) )
return(out)
}
例如如果我有一个向量 y
t = 1:1000
y = 10000*exp(sin(t*0.02))
和卷积核K
K = dexp(1:41, rate=1/5)
y
与卷积核K
的线性卷积将是linconv(y, K)
现在我想计算我的线性卷积,而不是使用矩阵乘法,方法是设置一个带状卷积矩阵,其中包含我的卷积核 K
的时移副本作为列。理想情况下,我想使用 Matrix::bandSparse
来设置这个带状卷积矩阵。有没有人碰巧知道 bandSparse
的正确语法来完成这个,并获得与使用上面的线性卷积函数相同的结果,但使用矩阵乘法代替?现在我做了
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)
}
但我想在上面设置矩阵 X
而不是使用 Matrix::bandSparse
。有人知道这样做的正确语法吗?
终于找到正确的语法了:
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)
}
在 R 中,我可以使用
计算向量y
与卷积核 K
的线性卷积
linconv = function (y,K) { K=K/sum(K)
out=convolve(y, K, type="open")
out=head(out,-(length(out)-length(y)) )
return(out)
}
例如如果我有一个向量 y
t = 1:1000
y = 10000*exp(sin(t*0.02))
和卷积核K
K = dexp(1:41, rate=1/5)
y
与卷积核K
的线性卷积将是linconv(y, K)
现在我想计算我的线性卷积,而不是使用矩阵乘法,方法是设置一个带状卷积矩阵,其中包含我的卷积核 K
的时移副本作为列。理想情况下,我想使用 Matrix::bandSparse
来设置这个带状卷积矩阵。有没有人碰巧知道 bandSparse
的正确语法来完成这个,并获得与使用上面的线性卷积函数相同的结果,但使用矩阵乘法代替?现在我做了
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)
}
但我想在上面设置矩阵 X
而不是使用 Matrix::bandSparse
。有人知道这样做的正确语法吗?
终于找到正确的语法了:
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)
}