Beklenmedik Konvolusyon

Question

ben R aşağıdaki konvolüsyonunu uygulamaya çalışırken, ancak beklenen sonucun elde edemiyorum Sonuçlar:

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

burada $ S [i] $ spektral yoğunlukları bir vektörü (Lorentz sinyal/NMR spektrum), ve $ i \ in [1, N] $ Burada $ N $ veri noktalarının sayısıdır (gerçek örneklerde, belki de 32K değerleri). Bu, Jacob, Deborde ve Moing'deki denklem 1, Analitik Biyanalitik Kimya (2013) 405: 5049-5061 (DOI 10.1007/s00216-013-6852-y).

$ SDL _ {\ Sigma} $ (kağıt denklem 2 göre) aşağıdaki gibi uygulayan bir Lorentz eğrisi, 2. türevini hesaplayan bir fonksiyonudur:

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 olduğu maksimum maksimum genişlikte ve x0 Lorentzian sinyalinin merkezidir.

SDL'un düzgün çalıştığına inanıyorum (çünkü döndürülen değerler ampirik Savitzky-Golay 2. türevi gibi bir şekle sahiptir). Benim sorunum olarak yazdım _ {\ sigma} $ $, C uygulamakla geçerli: referans Başına

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

, 2 türevi hesaplanması zaman kazanmak için yaklaşık 2000 veri noktaları pencereye sınırlıdır. Bu bölüm iyi çalışıyor. Sadece önemsiz çarpıklıklar üretmelidir.

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

Gördüğünüz gibi, konvolüsyonunu benzemez (mavi) SDL kullanılarak CP gelen büklüm ürünü: Burada

tüm sürecin bir gösteri ve sorundur SG yöntemini kullanarak CP ürününden (kırmızı hariç, ölçek hariç). SDL yöntemini kullanarak sonuçların benzer bir şekle sahip olması, ancak farklı bir ölçek olması beklenir.

Şimdiye kadar bana takıldınız, a) Teşekkürler, ve b) neyin yanlış olduğunu görebiliyor musunuz? Şüphesiz, temel bir yanlış anlama var.

Answer 1

El ile gerçekleştirdiğiniz evrişişlemeyle ilgili birkaç sorun var. "Savitzky – Golay filtresi" here için Wikipedia sayfasında tanımlandığı gibi konvolusyon işlevine bakarsanız, başvurulan denklemdeki S[i] terimiyle çakışan toplam içinde y[j+i] terimini görürsünüz. Başvurulan denkleminizin yanlış/yanlış yazılmış olabileceğine inanıyorum.

İşlevinizi aşağıdaki gibi değiştirdim ve uygulamamın tamamen doğru olduğundan emin olmadığım halde, sgolayfilt() sürümü ile aynı şekli oluşturmak için şimdi çalışıyor gibi görünüyor. sigma seçiminin önemli olduğunu ve sonuçta oluşan şekli etkilediğini unutmayın. Başlangıçta aynı şekli almazsanız, sigma parametresini önemli ölçüde değiştirmeye çalışın.

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

Kodunuzla ilgili birkaç sorun var.CPL’yi hesaplarken, CPL’yi hesaplarken, denkleminizdeki toplamı $ C _ {\ sigma} $ için yapmaya çalışıyorsunuz ama bu toplama, evrişimin tanımıdır.

Aslında SDL'yi hesaplarken x0 değerini değiştiriyorsunuz, ancak bu değer Lorentzian'ın ortalamasıdır ve sabit olmalıdır (bu durumda 0'dır).

Son olarak, konvolüsyon sınırlarını hesaplamak ve bu testi çalışıyor orijinal sınırları

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

ile sinyal çekebilir. Beklenen şeklinde

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) 

sonuçları:

cpSDL

Ben de SDL tanımı doğru olduğunu tamamen emin değilim. This article, bir lorentzinin ikinci türevi için çok daha karmaşık bir formüle sahiptir.