Como posso alinhar / sincronizar dois sinais?

21

Estou fazendo alguma pesquisa, mas fiquei paralisado na fase de análise (deveria ter prestado mais atenção às minhas aulas de estatísticas).

Eu coletei dois sinais simultâneos: vazão integrada para volume e alteração na expansão torácica. Eu gostaria de comparar os sinais e, finalmente, esperar obter volume do sinal de expansão torácica. Mas primeiro tenho que alinhar / sincronizar meus dados.

Como a gravação não inicia exatamente ao mesmo tempo e a expansão do tórax é capturada por períodos mais longos, preciso encontrar os dados que correspondem aos meus dados de volume no conjunto de dados de expansão do tórax e ter uma medida de quão bem eles estão alinhados. Não tenho muita certeza de como fazer isso se os dois sinais não iniciarem exatamente ao mesmo tempo ou entre dados em diferentes escalas e em diferentes resoluções.

Anexei um exemplo dos dois sinais ( https://docs.google.com/spreadsheet/ccc?key=0As4oZTKp4RZ3dFRKaktYWEhZLXlFbFVKNmllbGVXNHc ), informe-me se houver mais alguma coisa que eu possa fornecer.

person157
fonte
Não sei o suficiente para dar uma resposta e não tenho certeza de que isso resolva a questão, mas uma abordagem dos sinais de sincronização é chamada "registro", que é um subconjunto da análise de dados funcionais. Este tópico é discutido no livro da Ramsey e da FDA sobre Silverman. A idéia básica é que os sinais observados possam ser "distorcidos" (por exemplo, se estivéssemos interessados ​​na mecânica da maneira como as pessoas mastigam, mas temos dados sobre mastigando em velocidades diferentes - o eixo do tempo é "distorcido" neste caso) e O registro tenta definir o sinal subjacente em uma escala comum "sem distorção".
Macro
1
Você já coletou todos os seus dados? Este é um assunto piloto? Se você está apenas começando, eu gostaria de separar o sinal do seu cinto e usá-lo como gatilho (ou mesmo para marcar um carimbo de data / hora) na sua gravação de fluxo. Geralmente, os sistemas de aquisição têm essa capacidade com uma porta auxiliar ou acionadora. Tenho certeza de que existem maneiras de diferenciá-lo usando seus dados como sugerido por Macro, mas adicionar essa etapa extra exigirá muita suposição.
jonsca
1
Eu acho que você deseja estimar apenas um atraso fixo. Você pode usar a correlação cruzada, conforme descrito aqui: stats.stackexchange.com/questions/16121/…
thias th
1
Você pode fazer esta pergunta no dsp.SE, onde também pensa em sincronização de sinais.
precisa
1
@Thias está correto, mas parece que a primeira série deve ser reamostrada para ter intervalos comuns.
whuber

Respostas:

24

A pergunta pergunta como encontrar a quantidade pela qual uma série temporal ("expansão") fica atrasada para outra ("volume") quando a série é amostrada em intervalos regulares, mas diferentes .

Nesse caso, ambas as séries exibem um comportamento razoavelmente contínuo, como as figuras mostrarão. Isso implica (1) pouca ou nenhuma suavização inicial pode ser necessária e (2) a reamostragem pode ser tão simples quanto a interpolação linear ou quadrática. Quadrático pode ser um pouco melhor devido à suavidade. Após a reamostragem, o atraso é encontrado maximizando a correlação cruzada , conforme mostrado no encadeamento. Para duas séries de dados amostrados por offset, qual é a melhor estimativa do offset entre eles? .


Para ilustrar , podemos usar os dados fornecidos na pergunta, empregando Ro pseudocódigo. Vamos começar com a funcionalidade básica, correlação cruzada e reamostragem:

cor.cross <- function(x0, y0, i=0) {
  #
  # Sample autocorrelation at (integral) lag `i`:
  # Positive `i` compares future values of `x` to present values of `y`';
  # negative `i` compares past values of `x` to present values of `y`.
  #
  if (i < 0) {x<-y0; y<-x0; i<- -i}
  else {x<-x0; y<-y0}
  n <- length(x)
  cor(x[(i+1):n], y[1:(n-i)], use="complete.obs")
}

Este é um algoritmo bruto: um cálculo baseado em FFT seria mais rápido. Mas para esses dados (envolvendo cerca de 4000 valores), é bom o suficiente.

resample <- function(x,t) {
  #
  # Resample time series `x`, assumed to have unit time intervals, at time `t`.
  # Uses quadratic interpolation.
  #
  n <- length(x)
  if (n < 3) stop("First argument to resample is too short; need 3 elements.")
  i <- median(c(2, floor(t+1/2), n-1)) # Clamp `i` to the range 2..n-1
  u <- t-i
  x[i-1]*u*(u-1)/2 - x[i]*(u+1)*(u-1) + x[i+1]*u*(u+1)/2
}

Baixei os dados como um arquivo CSV separado por vírgula e retirei o cabeçalho. (O cabeçalho causou alguns problemas para o R que eu não queria diagnosticar.)

data <- read.table("f:/temp/a.csv", header=FALSE, sep=",", 
                    col.names=c("Sample","Time32Hz","Expansion","Time100Hz","Volume"))

NB Esta solução assume que cada série de dados está em ordem temporal, sem lacunas em nenhuma delas. Isso permite usar índices nos valores como proxies para o tempo e escalar esses índices pelas frequências de amostragem temporal para convertê-los em tempos.

Acontece que um ou ambos os instrumentos flutuam um pouco com o tempo. É bom remover essas tendências antes de prosseguir. Além disso, como há um afunilamento do sinal de volume no final, devemos cortá-lo.

n.clip <- 350      # Number of terminal volume values to eliminate
n <- length(data$Volume) - n.clip
indexes <- 1:n
v <- residuals(lm(data$Volume[indexes] ~ indexes))
expansion <- residuals(lm(data$Expansion[indexes] ~ indexes)

Reamostro as séries menos frequentes para obter o máximo de precisão possível no resultado.

e.frequency <- 32  # Herz
v.frequency <- 100 # Herz
e <- sapply(1:length(v), function(t) resample(expansion, e.frequency*t/v.frequency))

Agora a correlação cruzada pode ser calculada - por eficiência, buscamos apenas uma janela razoável de defasagens - e o atraso em que o valor máximo é encontrado pode ser identificado.

lag.max <- 5       # Seconds
lag.min <- -2      # Seconds (use 0 if expansion must lag volume)
time.range <- (lag.min*v.frequency):(lag.max*v.frequency)
data.cor <- sapply(time.range, function(i) cor.cross(e, v, i))
i <- time.range[which.max(data.cor)]
print(paste("Expansion lags volume by", i / v.frequency, "seconds."))

A saída nos diz que a expansão diminui o volume em 1,85 segundos. (Se os últimos 3,5 segundos de dados não foram cortados, a saída seria 1,84 segundos.)

É uma boa ideia verificar tudo de várias maneiras, de preferência visualmente. Primeiro, a função de correlação cruzada :

plot(time.range * (1/v.frequency), data.cor, type="l", lwd=2,
     xlab="Lag (seconds)", ylab="Correlation")
points(i * (1/v.frequency), max(data.cor), col="Red", cex=2.5)

gráfico de correlação cruzada

Em seguida, vamos registrar as duas séries no tempo e plotá-las juntas nos mesmos eixos .

normalize <- function(x) {
  #
  # Normalize vector `x` to the range 0..1.
  #
  x.max <- max(x); x.min <- min(x); dx <- x.max - x.min
  if (dx==0) dx <- 1
  (x-x.min) / dx
}
times <- (1:(n-i))* (1/v.frequency)
plot(times, normalize(e)[(i+1):n], type="l", lwd=2, 
     xlab="Time of volume measurement, seconds", ylab="Normalized values (volume is red)")
lines(times, normalize(v)[1:(n-i)], col="Red", lwd=2)

Parcelas registradas

Parece muito bom! No entanto, podemos ter uma noção melhor da qualidade do registro com um gráfico de dispersão . Eu vario as cores pelo tempo para mostrar a progressão.

colors <- hsv(1:(n-i)/(n-i+1), .8, .8)
plot(e[(i+1):n], v[1:(n-i)], col=colors, cex = 0.7,
     xlab="Expansion (lagged)", ylab="Volume")

Gráfico de dispersão

Procuramos os pontos a serem rastreados ao longo de uma linha: variações que refletem não linearidades na resposta com atraso de expansão do volume. Embora existam algumas variações, elas são bem pequenas. No entanto, como essas variações mudam ao longo do tempo pode ter algum interesse fisiológico. O maravilhoso das estatísticas, especialmente seu aspecto exploratório e visual, é como elas tendem a criar boas perguntas e idéias, além de respostas úteis .

whuber
fonte
1
Santo inferno, você é incrível. Correlação cruzada é exatamente o que eu estava imaginando (eu sabia que tinha que haver um nome para isso), mas sua resposta / explicação foi além. Muito obrigado!
person157
Não tenho tempo para uma explicação completa agora, mas uma ótima conta aparece nos livros "Receitas numéricas". Por exemplo, olhar para o capítulo 13.2, "Correlação e Autocorrelação Usando a FFT," em Numerical Recipes in C . Você também pode examinar a acffunção de R.
whuber
Novo em 'r', por favor, seja gentil: a função 'normalizar' usada na plotagem combinada (2ª última plotagem) não funcionará para mim. Existe uma atualização para essa função desde que esta resposta foi publicada?
CmKndy
1
@CmKndy Eu também era novo Rquando postei essa resposta e esqueci de fornecer uma definição para essa função. Aqui está o original:normalize <- function(x) { x.max <- max(x); x.min <- min(x); dx <- x.max - x.min; if (dx==0) dx <- 1; (x-x.min) / dx }
whuber
Perfeito, obrigado @whuber. Se você pudesse postar uma resposta como esta, quando você era novo para R, estou ainda mais recente do que eu pensava;)
CmKndy