Unerwartete Convolution Ergebnisse

Question

Ich versuche, die folgende Faltung in R, zu implementieren, aber nicht das erwartete Ergebnis bekommen:

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

wobei $ S [i] $ ein Vektor der spektralen Intensitäten ist (ein Lorentz-Signal/NMR-Spektrum), und $ i \ in [1, N] $ wobei $ N $ die Anzahl der Datenpunkte ist (in tatsächlichen Beispielen vielleicht 32K Werte). Dies ist Gleichung 1 in Jacob, Deborde und Moing, Analytische Bioanalytische Chemie (2013) 405: 5049-5061 (DOI 10.1007/s00216-013-6852-y).

$ SDL _ {\ sigma} $ ist eine Funktion, die 2. Ableitung einer Lorentz-Kurve zu berechnen, die I umgesetzt habe, wie folgt (basierend auf der Gleichung 2 in dem Papier):

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 ist Peak-Breite bei der Hälfte des Maximums und x0 ist das Zentrum des Lorentz-Signals.

Ich glaube, dass SDL korrekt funktioniert (weil die zurückgegebenen Werte eine Form wie die empirische Savitzky-Golay 2. Ableitung haben). Mein Problem ist, _ $ C mit der Implementierung {\ sigma} $, die ich geschrieben habe:

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

Pro der Referenz, die Berechnung der zweiten Ableitung zu einem Fenster von etwa 2000 Datenpunkten begrenzte Zeit zu sparen. Ich denke, dieser Teil funktioniert gut. Es sollte nur triviale Verzerrungen erzeugen.

Hier ist eine Demonstration des gesamten Prozesses und das 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") 

Wie Sie sehen können, das Faltungsprodukt von CPSDL mit (in blau) nicht die Faltung nicht ähnelt Produkt von CP mit der SG Methode (in rot, das ist richtig, außer für die Skala). Ich erwarte, dass die Ergebnisse der Verwendung der SDL Methode eine ähnliche Form aber eine andere Skala haben sollten.

Wenn Sie bis jetzt bei mir geblieben sind, a) danke, und b) können Sie sehen, was los ist? Zweifellos habe ich ein grundlegendes Missverständnis.

Answer 1

Es gibt einige Probleme mit der manuellen Faltung, die Sie durchführen. Wenn Sie die Faltungsfunktion wie auf der Wikipedia-Seite für "Savitzky-Golay-Filter" here definiert betrachten, sehen Sie den y[j+i] Begriff innerhalb der Summierung, der mit dem S[i] Begriff in der von Ihnen referenzierten Gleichung in Konflikt steht. Ich glaube, dass Ihre referenzierte Gleichung falsch/falsch ist.

Ich modifizierte Ihre Funktion wie folgt und es scheint jetzt zu funktionieren, um die gleiche Form wie die sgolayfilt() Version zu produzieren, obwohl ich nicht sicher bin, ob meine Implementierung völlig korrekt ist. Beachten Sie, dass die Auswahl von sigma wichtig ist und die resultierende Form beeinflusst. Wenn Sie anfangs nicht die gleiche Form erhalten, versuchen Sie, den Parameter sigma zu optimieren.

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

Es gibt ein paar Probleme mit Ihrem Code.In CP, wenn Sie SDL berechnen, sieht es so aus, als würden Sie versuchen, die Summe in Ihrer Gleichung für $ C _ {\ sigma} $ zu berechnen, aber diese Summe ist die Definition der Faltung.

Wenn Sie SDL berechnen, ändern Sie den Wert von x0, aber dieser Wert ist der Mittelwert des Lorentz und sollte konstant sein (in diesem Fall ist es 0).

Schließlich können Sie die Grenzen der Faltung berechnen und

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

Ausführen diesen Tests das Signal mit den ursprünglichen Grenzen herausziehen.

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) 

Ergebnisse in der erwarteten Form:

cpSDL

Ich bin auch nicht ganz sicher, dass Ihre SDL Definition korrekt ist. This article hat eine viel komplexere Formel für die zweite Ableitung eines Lorentz.