意外的卷积结果

Question

我正在尝试在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 和 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) 你能看出有什么问题吗？毫无疑问，我有一个根本性的误解。

Answer 1

您正在执行的手动卷积存在一些问题。如果您查看“Savitzky–Golay 滤波器”here 的 Wikipedia 页面上定义的卷积函数，您会在总和中看到 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)
}

Answer 2

您的代码存在几个问题。在 CP 中，当您计算 SDL 时，看起来您正试图在方程中对 $C_{\sigma}$ 进行求和，但这个求和是卷积的定义。

当您实际计算 SDL 时，您正在更改 x0 的值，但该值是洛伦兹的平均值并且应该是常数（在本例中为 0）。

最后，你可以计算卷积的边界，并用原始边界拉出信号

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)) {
        sdl[i] <- SDL(X[i], 0, sigma = sigma)
        }
    }

if (method == "SG") {
    sdl <- sgolayfilt(S, m = 2)     
    }

# Now convolve!  There is a built-in function for this!
cp <- convolve(S, sdl, type = "open")
shift <- floor(length(S)/2) #both signals are the same length and symmetric around zero
                             #so the fist point of the convolution is half the signal 
                             #before the first point of the signal   
print(shift)      
cp <- cp[shift:(length(cp)-shift)]
return (cp)
}

运行此测试。

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)

#
# Plot the original data, compare to convolution product

plot(x, cpSDL)

产生预期的形状：

cpSDL

我也不完全确定您的 SDL 定义是否正确。 This article 有一个更复杂的洛伦兹二阶导数公式。