Неожиданные результаты свертки

Я пытаюсь реализовать следующую свертку в R, но не получая ожидаемого результата:

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

где $S [i] $- вектор спектральных интенсивностей (Лоренцевый спектр сигнала/ЯМР) и $i\in [1, N] $, где $N $- количество точек данных (в реальных примерах, возможно 32K значений). Это уравнение 1 в Jacob, Deborde and Moing, Analytical Bioanalytical Chemistry (2013) 405: 5049-5061 (DOI 10.1007/s00216-013-6852-y).

$SDL _ {\ sigma} $- функция для вычисления 2-й производной лоренцевой кривой, которую я выполнил следующим образом (на основе уравнения 2 в статье):

sigma - ширина пика в половине максимума, а x0 - центр лоренцевого сигнала.

Я считаю, что SDL работает правильно (потому что возвращаемые значения имеют форму, подобную эмпирической второй производной Савицки-Голея). Моя проблема заключается в реализации $C_ {\ sigma} $, который я написал как:

Для справки, вычисление 2-й производной ограничено окном около 2000 точек данных для экономии времени. Я думаю, что эта часть отлично работает. Он должен давать только тривиальные искажения.

Вот демонстрация всего процесса и проблемы:

Как вы можете видеть, продукт свертки из CP с использованием SDL (синим цветом) не похож на продукт свертки из CP с использованием метода SG (красным, что является правильным, за исключением масштаб). Я ожидаю, что результаты использования метода SDL должны иметь одинаковую форму, но в другом масштабе.

Если вы застряли со мной до сих пор, а) спасибо, и б) можете ли вы понять, что случилось? Несомненно, у меня есть фундаментальное недоразумение.

Ответы

Ответ 1

Есть несколько проблем с ручной сверткой, которую вы выполняете. Если вы посмотрите на функцию свертки, определенную на странице Википедии "Фильтр Савицкого-Голея" здесь, вы увидите термин 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)
}

Ответ 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 верное. Эта статья содержит гораздо более сложную формулу для второй производной лоренциана.