意外的卷积结果

我正在尝试实现以下卷积R,但没有得到预期的结果:

$$ C _ {\ sigma} [i] =\sum\limits_ {k = -P} ^ P SDL _ {\ sigma} [ik,i]\centerdot S [i] $$

其中$ S [i] $是光谱强度的矢量(洛伦兹信号/核磁共振光谱),而$ i\in [1,N] $其中$ N $是数据点的数量(在实际例子中,可能是32K)值).这是Jacob,Deborde和Moing,Analytical Bioanalytical Chemistry(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)你能看出什么是错的吗？毫无疑问,我有一个根本的误解.