bất ngờ Convolution Kết quả

Question

Tôi đang cố gắng để thực hiện chập sau trong R, nhưng không nhận được kết quả mong đợi:

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

nơi $ S [i] $ là một vector của cường độ quang phổ (một tín hiệu Lorentzian/NMR), và $ i \ in [1, N] $ trong đó $ N $ là số điểm dữ liệu (trong các ví dụ thực tế, có lẽ là giá trị 32K). Đây là phương trình 1 trong Jacob, Deborde và Moing, Hóa phân tích sinh học phân tích (2013) 405: 5049-5061 (DOI 10.1007/s00216-013-6852-y).

$ SDL _ {\ sigma} $ là một hàm để tính đạo hàm thứ 2 của một đường cong Lorentzian, mà tôi đã thực hiện như sau (dựa trên phương trình 2 trong giấy):

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 là chiều rộng đỉnh tối đa một nửa và x0 là trung tâm của tín hiệu Lorentzian.

Tôi tin rằng SDL hoạt động chính xác (vì giá trị trả về có hình dạng giống như đạo hàm thứ 2 của Savitzky-Golay). Vấn đề của tôi là với việc thực hiện $ C _ {\ sigma} $, mà tôi đã viết như sau:

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

mỗi tài liệu tham khảo, việc tính toán đạo hàm thứ 2 được giới hạn trong một cửa sổ của khoảng 2000 điểm dữ liệu để tiết kiệm thời gian. Tôi nghĩ rằng phần này hoạt động tốt. Nó chỉ tạo ra những biến dạng tầm thường.

Đây là một cuộc biểu tình của toàn bộ quá trình và các vấn đề:

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

Như bạn có thể thấy, các sản phẩm chập từ CP sử dụng SDL (màu xanh) không giống chập sản phẩm từ CP sử dụng phương pháp SG (màu đỏ, chính xác, ngoại trừ tỷ lệ). Tôi hy vọng kết quả từ việc sử dụng phương pháp SDL phải có hình dạng tương tự nhưng có quy mô khác nhau.

Nếu bạn đã mắc kẹt với tôi cho đến nay, a) cảm ơn và b) bạn có thể thấy điều gì sai? Không nghi ngờ gì, tôi có một sự hiểu lầm cơ bản.

Answer 1

Có một số vấn đề với quá trình xoay vòng thủ công bạn đang thực hiện. Nếu bạn nhìn vào hàm convolution như được định nghĩa trên trang Wikipedia cho "bộ lọc Savitzky – Golay" here, bạn sẽ thấy cụm từ y[j+i] bên trong tổng kết xung đột với cụm từ S[i] trong phương trình bạn đã tham chiếu. Tôi tin rằng phương trình tham chiếu của bạn có thể không đúng/lỗi.

Tôi đã sửa đổi chức năng của bạn như sau và nó có vẻ hoạt động ngay bây giờ để tạo ra hình dạng giống như phiên bản sgolayfilt(), mặc dù tôi không chắc chắn việc triển khai của tôi là hoàn toàn chính xác. Lưu ý rằng sự lựa chọn của sigma là quan trọng và không ảnh hưởng đến hình dạng kết quả. Nếu bạn không nhận được cùng một hình dạng ban đầu, hãy thử tinh chỉnh đáng kể thông số 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

Có một số vấn đề với mã của bạn.Trong CP khi bạn đang tính toán SDL, có vẻ như bạn đang cố gắng thực hiện tổng kết trong phương trình của bạn với $ C _ {\ sigma} $ nhưng tổng kết này là định nghĩa của convolution.

Khi bạn đang tính toán SDL, bạn đang thay đổi giá trị của x0, nhưng giá trị này là giá trị trung bình của Lorentzian và phải là hằng số (trong trường hợp này là 0).

Cuối cùng, bạn có thể tính toán giới hạn của chập và kéo ra các tín hiệu với các giới hạn ban đầu

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

Chạy thử nghiệm này.

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) 

kết quả trong hình dạng mong muốn:

cpSDL

Tôi cũng không hoàn toàn chắc chắn rằng định nghĩa SDL của bạn là đúng. This article có công thức phức tạp hơn nhiều đối với đạo hàm thứ hai của một lorentzian.