Unexpected Convolution Results

Question

I'm trying to implement the following convolution in R, but not getting the expected result:

$$ C_{\sigma}[i]=\sum\limits_{k=-P}^P SDL_{\sigma}[i-k,i] \centerdot S[i] $$

where $S[i]$ is a vector of spectral intensities (a Lorentzian signal / NMR spectrum), and $i \in [1,N]$ where $N$ is the number of data points (in actual examples, perhaps 32K values). This is equation 1 in Jacob, Deborde and Moing, Analytical Bioanalytical Chemistry (2013) 405:5049-5061 (DOI 10.1007/s00216-013-6852-y).

$SDL_{\sigma}$ is a function to compute the 2nd derivative of a Lorentzian curve, which I have implemented as follows (based on equation 2 in the paper):

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 is the peak width at half maximum, and x0 is the center of the Lorentzian signal.

I believe that SDL works correctly (because the returned values have a shape like the empirical Savitzky-Golay 2nd derivative). My problem is with implementing $C_{\sigma}$, which I have written as:

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)
    }

Per the reference, the computation of the 2nd derivative is limited to a window of about 2000 data points to save time. I think this part works fine. It should produce only trivial distortions.

Here is a demonstration of the entire process and the problem:

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")

As you can see, the convolution product from CP using SDL (in blue) doesn't resemble the convolution product from CP using the SG method (in red, which is correct, except for scale). I expect the results from using the SDL method should have a similar shape but a different scale.

If you've stuck with me so far, a) thanks, and b) can you see what's wrong? No doubt, I have a fundamental misunderstanding.

Answer 1

There are a few issues with the manual convolution you are performing. If you look at the convolution function as defined on the Wikipedia page for "Savitzky–Golay filter" here, you'll see the y[j+i] term inside the summation which conflicts with the S[i] term in the equation you referenced. I believe your referenced equation may be incorrect / typoed.

I modified your function as follows and it seems to work now to produce the same shape as the sgolayfilt() version, although I'm not certain my implementation is totally correct. Note that the choice of sigma is important and does affect the resulting shape. If you don't get the same shape initially, trying significantly tweaking the sigma parameter.

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

There are a couple of issues with your code. In CP when you are calculating SDL it looks like you are trying to do the summation in your equation for $C_{\sigma}$ but this summation is the definition of the convolution.

When you are actually calculating SDL you are changing the value of x0, but this value is the mean of the Lorentzian and should be constant (in this case it is 0).

Finally, you can calculate the bounds of the convolution and pull out the signal with the original bounds

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)
}

Running this test.

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)

results in the expected shape:

cpSDL

I'm also not entirely sure that your SDL definition is correct. This article has a much more complex formula for the second derivative of a lorentzian.