imprevisto Convoluzione Risultati

Question

che sto cercando di implementare il seguente convoluzione in R, ma non ottenere il risultato atteso:

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

dove $ S [i] $ è un vettore di intensità spettrali (un segnale Lorentzian/spettro NMR), e $ i \ in [1, N] $ dove $ N $ è il numero di punti dati (in esempi reali, forse 32K valori). Questa è l'equazione 1 in Jacob, Deborde e Moing, Chimica bioanalitica analitica (2013) 405: 5049-5061 (DOI 10.1007/s00216-013-6852-y).

$ SDL _ {\ sigma} $ è una funzione per calcolare la seconda derivata di una curva di Lorentz, che ho implementato come segue (sulla base degli equazione 2 nel documento):

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 è la larghezza di picco a metà massimo e x0 è il centro del segnale di Lorentzian.

Credo che SDL funzioni correttamente (perché i valori restituiti hanno una forma simile alla derivata empirica di Savitzky-Golay 2). Il mio problema è di attuare $ C _ {\ sigma} $, che ho scritto come:

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 il riferimento, il calcolo della seconda derivata è limitata a una finestra di circa 2000 punti di dati per risparmiare tempo. Penso che questa parte funzioni bene. Dovrebbe produrre solo distorsioni banali.

Ecco una dimostrazione di tutto il processo e il problema:

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

Come si può vedere, il prodotto di convoluzione da CP utilizzando SDL (in blu) non assomiglia la convoluzione prodotto da CP utilizzando il metodo SG (in rosso, che è corretto, fatta eccezione per la scala). Mi aspetto che i risultati dell'utilizzo del metodo SDL debbano avere una forma simile ma una scala diversa.

Se sei rimasto con me finora, a) grazie, e b) riesci a vedere cosa c'è che non va? Senza dubbio, ho un malinteso fondamentale.

Answer 1

Ci sono alcuni problemi con la convoluzione manuale che si sta eseguendo. Se si osserva la funzione di convoluzione come definita nella pagina di Wikipedia per "filtro Savitzky-Golay" here, verrà visualizzato il termine y[j+i] all'interno della somma che è in conflitto con il termine S[i] nell'equazione a cui si fa riferimento. Credo che la tua equazione di riferimento potrebbe essere errata/tipizzata.

Ho modificato la funzione come segue e sembra che ora funzioni per produrre la stessa forma della versione sgolayfilt(), anche se non sono sicuro che la mia implementazione sia completamente corretta. Si noti che la scelta di sigma è importante e influisce sulla forma risultante. Se inizialmente non si ottiene la stessa forma, provare a modificare in modo significativo il parametro 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

Ci sono un paio di problemi con il codice.Nel CP quando si calcola SDL sembra che si stia tentando di eseguire la somma nella propria equazione per $ C _ {\ sigma} $ ma questa somma è la definizione della convoluzione.

Quando si sta effettivamente calcolando SDL si sta cambiando il valore di x0, ma questo valore è la media del Lorentzian e dovrebbe essere costante (in questo caso è 0).

Infine, è possibile calcolare i limiti della circonvoluzione ed estrarre il segnale con i limiti originali

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

eseguire questo 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) 

risultati nella forma prevista:

cpSDL

Io non sono anche del tutto sicuro che la tua definizione SDL è corretta. This article ha una formula molto più complessa per la derivata seconda di un lorenziano.