Estimativa de atraso de tempo dos sinais do osciloscópio usando correlação cruzada

Gravei 2 sinais de um osciloscópio. Eles se parecem com isso:

Quero medir o atraso de tempo entre eles no Matlab. Cada sinal possui 2000 amostras com uma frequência de amostragem de 2001000,5.

Os dados estão em um arquivo csv. É isso que eu tenho até agora.
Apaguei os dados de tempo do arquivo csv para que apenas os níveis de tensão estejam no arquivo csv.

x1 = csvread('C://scope1.csv');
x2 = csvread('C://scope2.csv');  
cc = xcorr(x1,x2);
plot(cc);  

Isso fornece este resultado:

Pelo que li, preciso fazer a correlação cruzada desses sinais e isso deve me dar um pico em relação ao atraso de tempo. No entanto, quando tomo a correlação cruzada desses sinais, obtenho um pico em 2000, que sei que não está correto. O que devo fazer com esses sinais antes de cruzá-los? Apenas procurando alguma direção.

EDIT: depois de remover o deslocamento DC, este é o resultado que agora estou obtendo:

Existe uma maneira de limpar isso para obter um atraso de tempo mais definido?

EDIT 2: Aqui estão os arquivos:
http://dl.dropbox.com/u/10147354/scope1col.csv
http://dl.dropbox.com/u/10147354/scope2col.csv

@NickS

Como está longe de ser certo que o segundo sinal nas parcelas é de fato uma versão apenas atrasada do primeiro, outros métodos além da correlação cruzada clássica devem ser tentados. Isso ocorre porque a correlação cruzada (CC) é apenas um estimador de probabilidade máxima se o (s) seu (s) sinal (es) forem atrasados ​​versões um do outro. Nesse caso, eles claramente não são, para não dizer nada sobre a não estacionariedade deles.

Nesse caso, acredito que o que pode funcionar é uma estimativa de tempo da energia significativa dos sinais. É verdade que 'significativo' pode ou não pode ser um pouco subjetivo, mas acredito que, observando seus sinais do ponto de vista estatístico, seremos capazes de quantificar 'significativo' e partir daí.

Para esse fim, fiz o seguinte:

PASSO 1: Calcule os envelopes de sinal:

Este passo é simples, pois o valor absoluto da saída da Hilbert-Transform de cada um de seus sinais é calculado. Existem outros métodos para calcular envelopes, mas isso é bastante simples. Este método calcula essencialmente a forma analítica do seu sinal, ou seja, a representação fasorial. Quando você toma o valor absoluto, está destruindo a fase e somente depois da energia.

Além disso, como estamos buscando uma estimativa de atraso de tempo da energia de seus sinais, essa abordagem é garantida.

PASSO 2: Ruído com filtros mediais não lineares que preservam as bordas:

Este é um passo importante. O objetivo aqui é suavizar os envelopes de energia, mas sem destruir ou suavizar as bordas e aumentar os tempos de subida. Na verdade, existe um campo inteiro dedicado a isso, mas para nossos propósitos aqui, podemos simplesmente usar um filtro Medial não linear fácil de implementar . (Filtragem mediana). Essa é uma técnica poderosa porque, diferentemente da filtragem média , a filtragem medial não anulará suas bordas, mas ao mesmo tempo 'suavizará' o seu sinal sem degradação significativa das bordas importantes, pois em nenhum momento qualquer aritmética está sendo realizada em seu sinal. (desde que o comprimento da janela seja ímpar). Para o nosso caso aqui, selecionei um filtro medial de amostras de tamanho de janela 25:

PASSO 3: Remover Tempo: Construir Funções de Estimativa de Densidade de Kernel Gaussiana:

O que aconteceria se você olhasse para o gráfico acima de lado e não da maneira normal? Matematicamente falando, o que você obteria se projetasse todas as amostras de nossos sinais denoisados ​​no eixo da amplitude y? Ao fazer isso, conseguiremos remover o tempo, por assim dizer, e poderemos estudar apenas as estatísticas do sinal.

Intuitivamente, o que aparece na figura acima? Embora a energia do ruído seja baixa, ela tem a vantagem de ser mais 'popular'. Por outro lado, enquanto o envelope de sinal que possui energia é mais energético que o ruído, ele é fragmentado entre os limites. E se considerarmos a 'popularidade' como uma medida de energia? Isto é o que faremos com a implementação (básica) de uma Função de Densidade do Kernel (KDE), com um Kernel Gaussiano.

Para fazer isso, cada amostra é coletada e uma função gaussiana é construída usando seu valor como média, e uma largura de banda predefinida (variação) é selecionada a priori. Definir a variação do seu gaussiano é um parâmetro importante, mas você pode configurá-lo com base em estatísticas de ruído com base em sua aplicação e sinais típicos. (Eu só tenho dois arquivos para ativar). Se construirmos a Estimação do KDE, obteremos o seguinte gráfico:

Você pode pensar no KDE como uma forma contínua de um histograma, por assim dizer, e a variação como sua largura de caixa. No entanto, ele tem a vantagem de garantir um PDF suave, no qual podemos executar o primeiro e o segundo cálculos derivados. Agora que temos os KDEs gaussianos, podemos ver onde as amostras de ruído atingem seu pico de popularidade. Lembre-se de que o eixo x aqui representa as projeções de nossos dados no espaço de amplitude. Assim, podemos ver em quais limiares o ruído é o mais 'energético', e aqueles nos dizem quais limiares evitar.

No segundo gráfico, a primeira derivada dos KDEs gaussianos é obtida, e escolhemos a abcissa da primeira amostra após a primeira derivada após o pico da mistura de gaussianos para atingir um determinado valor próximo de zero. (Ou primeiro cruzamento de zero). Podemos usar esse método e ser 'seguros' porque nosso KDE foi construído com gaussianos suaves, com largura de banda razoável, e a primeira derivada dessa função suave e sem ruído foi usada. (Normalmente, as primeiras derivadas podem ser problemáticas em qualquer coisa, exceto nos sinais SNR altos, pois aumentam o ruído).

As linhas pretas mostram, então, em que limites seria sensato segmentar a imagem, de modo a evitar todo o ruído. Se aplicarmos nossos sinais originais, obteremos os seguintes gráficos, com as linhas pretas indicando o início da energia de nossos sinais:

$\delta{t} = 241$

Espero que isso tenha ajudado.

Existem alguns problemas ao fazer isso com autocorrelação

1. Deslocamento CC enorme (já corrigido)
2. Janela de tempo: o xcorr () do Matlab possui a convenção irritante de essencialmente "zerar" o sinal nas duas extremidades enquanto você desliza o intervalo de tempo. Ou seja, a janela de dados é uma função do intervalo de tempo. Isso criará uma forma triangular para qualquer sinal estacionário (incluindo ondas senoidais). As melhores opções são escolher uma janela de correlação para que o tamanho da janela mais o intervalo de tempo máximo se ajustem à sua janela total de dados ou normalizar a correlação cruzada pelo número de amostras não preenchidas.
3. Os dois sinais não parecem particularmente correlacionados comigo. A forma é um pouco semelhante, mas o espaçamento específico de picos e quedas é bem diferente, por isso duvido que mesmo uma correlação automática adequada traga muitas informações aqui.

Uma abordagem muito mais simples seria usar um detector de limiar para encontrar os pontos de partida e simplesmente usar a diferença entre esses pontos como atraso.

Como as pichenettes indicaram, neste caso, um pico no meio da saída indica 0 atraso. O deslocamento do pico do ponto médio é o seu atraso de tempo.

EDIT: Me preocupa que a correlação seja quase um triângulo perfeito. Isso indica para mim que a correlação cruzada não está normalizando a energia. Isso fornece um viés injusto para atrasos menores em relação a atrasos maiores. Eu modificaria sua chamada xcorr para "cc = xcorr (x1, x2, 'imparcial');".

Essa não é, em sua opinião, uma solução perfeita, porque os resultados de grande atraso agora são mais instáveis ​​do que os resultados de baixo atraso, porque são baseados em menos dados. Um grande pico nas extremidades pode ser falso pela mesma razão que você pode obter 100% de cara e sem coroa em apenas alguns lançamentos de moedas, embora seja extremamente improvável que ocorra em muitos lançamentos.

Como os outros apontaram, e parece que você percebeu com base na sua última edição da pergunta, não parece que a correlação cruzada forneça uma boa estimativa de atraso de tempo para os conjuntos de dados mostrados. A correlação mede a similaridade na forma entre duas séries temporais, deslizando uma pela outra por um intervalo de tempo e calculando um produto interno entre as duas séries a cada atraso. O resultado terá uma grande magnitude quando as duas séries forem qualitativamente semelhantes ou forem "correlacionadas" umas com as outras. É semelhante à forma como um produto interno de dois vetores é maior quando os dois vetores são apontados na mesma direção.

O problema com os dados que você mostrou é que (pelo menos para os trechos que podemos ver) não parece haver muita semelhança na forma. Não há atraso que você possa aplicar a um dos sinais para que se pareça com o outro, que é exatamente o que você está fazendo calculando a correlação cruzada entre eles.

Há casos em que a correlação cruzada é útil, no entanto. Digamos que seu segundo sinal foi realmente uma versão com alteração de tempo do original, mesmo com algum ruído adicional adicionado:

a = csvread('scope1col.csv');
a = a - mean(a);               % to remove DC offset
b = a(200:end) + sqrt(0.05)*randn(1801,1);
figure; subplot(211); plot(a); subplot(212); plot(b)

Agora não está claro imediatamente que os dois sinais estão relacionados por um atraso de tempo. No entanto, se fizermos a correlação cruzada, obtemos:

[c,lags] = xcorr(a,b);
igure; plot(lags,c); grid on; xlabel('Lag'); ylabel('Cross-correlation');

que mostra um pico no atraso correto de 200 amostras. A correlação pode ser uma ferramenta útil para determinar o atraso, quando aplicada a conjuntos de dados que contêm o tipo certo de similaridade.

Com base na sugestão de Muhammad, tentei fazer um script do Matlab. No entanto, não sou capaz de deduzir se ele constrói uma distribuição gaussiana com base nas variações e, em seguida, faz uma estimativa do KDE ou executa uma estimativa do KDE com suposição gaussiana.

Também é difícil inferir como ele traduz o tempo de compensação do KDE no domínio do tempo. Aqui está a minha tentativa. Qualquer usuário interessado em usar o script é livre e atualiza a versão aprimorada, se possível.

%% Initialising data

Ws1 = data1;
Ws2 = data2;
mWs1 = nanmean(Ws1);
mWs2 = nanmean(Ws2);
sdWs1 = nanstd(Ws1);
sdWs2 = nanstd(Ws2);

%% Computing the signal envelopes
Ws1d = Ws1 - mWs1;
Ws2d = Ws2 - mWs2;
h1 = abs(hilbert(Ws1d));
h2 = abs(hilbert(Ws2d));
figure();
subplot(211)
plot([Ws1d, h1])
subplot(212)
plot([Ws2d, h2])

%% Denoise the signal with edge preserving nonlinear medial filtering
w = 25;
mf1 = medfilt1(h1, w);
mf2 = medfilt1(h2, w);
figure();
subplot(211)
plot(mf1)
subplot(212)
plot(mf2)

<%% Remove time: construct the gaussian kernel density estimation functions>
% Using the kde from Matlab central directly on the filtered data
data1 = mf1;
[bw1, den1, xmesh1, cdf1] = kde(data1, 2^14);
der1 = diff(den1);
data2 = mf2;
[bw2, den2, xmesh2, cdf2] = kde(data2, 2^14);
der2 = diff(den2);
figure();
plot([der1, der2]);
legend('Sig1', 'Sig2')

% the other method as explained in Muhammad's post
for i = 1:length(mf1)
gf1(:,i) = mf1(i) + sdWs1*randn(1000,1);
gf2(:,i) = mf2(i) + sdWs2*randn(1000,1);
end
[bwM1, denM1, xmeshM1, cdfM1] = kde(gf1(:,1), 2.^11);
dd1 = diff(denM1);
[bwM2, denM2, xmeshM2, cdfM2] = kde(gf2(:,1), 2.^11);
dd2 = diff(denM2);
figure();
plot([dd1, dd2]);
legend('Sig1', 'Sig2')