我正在尝试在R中实现以下卷积，但未获得预期的结果：
$$
C _ {\ sigma} [i] = \ sum \ limits_ {k = -P} ^ P SDL _ {\ sigma} [i-k，i] \ centerdot S [i]
$$
其中$ S [i] $是频谱强度的向量（洛伦兹信号/ NMR谱），而$ i \ in [1，N] $，其中$ N $是数据点的数量（在实际示例中，可能是32K值）。这是Jacob，Deborde and Moing，分析生物分析化学（2013）405：5049-5061（DOI 10.1007 / s00216-013-6852-y）中的方程式1。
$ SDL _ {\ sigma} $是用于计算洛伦兹曲线的二阶导数的函数，我已如下实现（基于本文的公式2）：
SDL <- function(x, x0, sigma = 0.0005){
if (!sigma > 0) stop("sigma must be greater than zero.")
num <- 16 * sigma * ((12 * (x-x0)^2) - sigma^2)
denom <- pi * ((4 * (x - x0)^2) + sigma^2)^3
sdl <- num/denom
return(sdl)
}
sigma是最大宽度的一半，而x0是洛伦兹信号的中心。我相信
SDL可以正常工作（因为返回值的形状类似于经验Savitzky-Golay二阶导数）。我的问题是实现$ C _ {\ sigma} $，我将其编写为：CP <- function(S = NULL, X = NULL, method = "SDL", W = 2000, sigma = 0.0005) {
# S is the spectrum, X is the frequencies, W is the window size (2*P in the eqn above)
# Compute the requested 2nd derivative
if (method == "SDL") {
P <- floor(W/2)
sdl <- rep(NA_real_, length(X)) # initialize a vector to store the final answer
for(i in 1:length(X)) {
# Shrink window if necessary at each extreme
if ((i + P) > length(X)) P <- (length(X) - i + 1)
if (i < P) P <- i
# Assemble the indices corresponding to the window
idx <- seq(i - P + 1, i + P - 1, 1)
# Now compute the sdl
sdl[i] <- sum(SDL(X[idx], X[i], sigma = sigma))
P <- floor(W/2) # need to reset at the end of each iteration
}
}
if (method == "SG") {
sdl <- sgolayfilt(S, m = 2)
}
# Now convolve! There is a built-in function for this!
cp <- convolve(S, sdl, type = "open")
# The convolution has length 2*(length(S)) - 1 due to zero padding
# so we need rescale back to the scale of S
# Not sure if this is the right approach, but it doesn't affect the shape
cp <- c(cp, 0.0)
cp <- colMeans(matrix(cp, ncol = length(cp)/2)) # stackoverflow.com/q/32746842/633251
return(cp)
}
根据参考，二阶导数的计算限于大约2000个数据点的窗口以节省时间。我认为这部分工作正常。它只会产生微小的失真。
这是整个过程和问题的演示：
require("SpecHelpers")
require("signal")
# Create a Lorentzian curve
loren <- data.frame(x0 = 0, area = 1, gamma = 0.5)
lorentz1 <- makeSpec(loren, plot = FALSE, type = "lorentz", dd = 100, x.range = c(-10, 10))
#
# Compute convolution
x <- lorentz1[1,] # Frequency values
y <- lorentz1[2,] # Intensity values
sig <- 100 * 0.0005 # per the reference
cpSDL <- CP(S = y, X = x, sigma = sig)
sdl <- sgolayfilt(y, m = 2)
cpSG <- CP(S = y, method = "SG")
#
# Plot the original data, compare to convolution product
ylabel <- "data (black), Conv. Prod. SDL (blue), Conv. Prod. SG (red)"
plot(x, y, type = "l", ylab = ylabel, ylim = c(-0.75, 0.75))
lines(x, cpSG*100, col = "red")
lines(x, cpSDL/2e5, col = "blue")
如您所见，使用
CP的SDL的卷积乘积（蓝色）与使用CP方法的SG的卷积乘积（红色，除了比例尺，这是正确的）不相似。我希望使用SDL方法得到的结果应具有相似的形状但比例不同。如果您到目前为止一直对我不满意，a）谢谢，b）您能看到出什么问题了吗？毫无疑问，我有一个基本的误解。
最佳答案
您正在执行的手动卷积存在一些问题。如果查看Wikipedia页面上为“ Savitzky-Golay滤波器” here定义的卷积函数，您将在求和内看到与所引用方程式中的y[j+i]项冲突的S[i]项。我相信您引用的等式可能不正确/有错字。
我修改了您的函数，如下所示，现在看来可以生成与sgolayfilt()版本相同的形状，尽管我不确定我的实现是完全正确的。请注意，sigma的选择很重要，并且确实会影响最终的形状。如果最初没有得到相同的形状，请尝试大幅调整sigma参数。
CP <- function(S = NULL, X = NULL, method = "SDL", W = 2000, sigma = 0.0005) {
# S is the spectrum, X is the frequencies, W is the window size (2*P in the eqn above)
# Compute the requested 2nd derivative
if (method == "SDL") {
sdl <- rep(NA_real_, length(X)) # initialize a vector to store the final answer
for(i in 1:length(X)) {
bound1 <- 2*i - 1
bound2 <- 2*length(X) - 2*i + 1
P <- min(bound1, bound2)
# Assemble the indices corresponding to the window
idx <- seq(i-(P-1)/2, i+(P-1)/2, 1)
# Now compute the sdl
sdl[i] <- sum(SDL(X[idx], X[i], sigma = sigma) * S[idx])
}
}
if (method == "SG") {
sdl <- sgolayfilt(S, m = 2)
}
# Now convolve! There is a built-in function for this!
cp <- convolve(S, sdl, type = "open")
# The convolution has length 2*(length(S)) - 1 due to zero padding
# so we need rescale back to the scale of S
# Not sure if this is the right approach, but it doesn't affect the shape
cp <- c(cp, 0.0)
cp <- colMeans(matrix(cp, ncol = length(cp)/2)) # stackoverflow.com/q/32746842/633251
return(cp)
}
关于r - 意外的卷积结果，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/34032500/