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

12

Gravei 2 sinais de um osciloscópio. Eles se parecem com isso: insira a descrição da imagem aqui

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: insira a descrição da imagem aqui

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:
insira a descrição da imagem aqui

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

Nick Sinas
fonte
Como, exatamente, você está fazendo a correlação cruzada? Em resposta à sua pergunta direta, você não precisa fazer nada com seus sinais antes da correlação cruzada, embora em alguns casos a filtragem primeiro ajude a eliminar ruídos que possam distorcer os resultados.
Jim Clay
1
Por favor, poste o código que você usou e, mais importante, um gráfico do sinal de correlação cruzada. Algumas ferramentas / bibliotecas colocam a pontuação (lag = 0) no meio do gráfico; Não me lembro se o Matlab faz isso.
Pichenettes
@pichenettes: post atualizado
Nick Sinas
@JimClay: post atualizado
Nick Sinas
@NickS. Se os seus sinais estiverem perfeitamente alinhados, você obterá um pico no meio do gráfico cc. Portanto, o pico em 2000 significa sem demora. Agora, digamos que você tenha um atraso de 10 amostras, o que significa que o sinal2 está a 10 amostras do sinal1. Isso apenas moverá seu pico no CC de 2000 para 2010 (ou 1990). Portanto, o seu atraso corresponde tempo para a sua localização de pico real, menos 2000.
Spacey

Respostas:

11

@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.

insira a descrição da imagem aqui

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:

insira a descrição da imagem aqui

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:

insira a descrição da imagem aqui

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:

insira a descrição da imagem aqui

δt=241

Espero que isso tenha ajudado.

Spacey
fonte
Uau. Muito obrigado. Todas essas são novas técnicas que eu começarei a pesquisar. Existe alguma maneira de eu dar uma olhada no código do matlab que você usou?
precisa saber é o seguinte
Portanto, eu tenho as etapas 1 e 2 executadas no Matlab e meus resultados correspondem aos seus, mas estou tendo problemas com a etapa 3. Que funções você usou?
precisa saber é o seguinte
@NickS. Peça, .. e você receberá, envie-me um e-mail e posso enviar o código que usei.
Spacey
@ Mohammed Você poderia postar seu código para estimar o atraso. Enviei um e-mail a respeito desta ajuda importa então por favor
6

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.

Hilmar
fonte
4

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.

Jim Clay
fonte
Significando que os sinais não estão atrasados?
Nick SINAS
Não tenho certeza - onde está o pico? Percebo que está perto do meio, mas não está claro se está realmente no meio. Além disso, há um problema de normalização de energia que abordarei em uma edição da minha resposta.
Jim Clay
O parâmetro 'imparcial' definitivamente faz com que pareça melhor. mais do que eu esperaria. Vou continuar investigando isso. Obrigado.
precisa
@ JimClay Talvez Nick S, esteja correlacionando os envelopes de seus sinais e não os sinais reais (Nick é verdade?). Isso produziria (aproximadamente) essa forma triangular que eu imagino.
Spacey
2
@NickS. O comentário de Mohammad me fez olhar e perceber que você tem um enorme deslocamento de DC que está atrapalhando seus resultados. Subtraia a média de ambos os sinais e execute xcorr neles. Eu tentaria sem a opção "imparcial" primeiro.
Jim Clay
4

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)

insira a descrição da imagem aqui

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');

insira a descrição da imagem aqui

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.

Jason R
fonte
Alguma idéia para o que mais eu poderia fazer? Talvez outra técnica que não seja a correlação cruzada ou talvez algum tipo de filtro? Obrigado.
9788 Nick Sinas
@NickS. Eu também olhei para ele e eles não são cópias atrasadas um do outro. Dito isto, você deseja estimar o atraso de energia ? Eu acho que faria mais sentido, neste caso, VS atraso dos sinais . Se você nos contar mais sobre o canal / experimento subjacente que está executando, podemos saber mais sobre os possíveis caminhos.
Spacey
@Mohammad Thanks. O canal subjacente é o aço. Alguma idéia de como estimar o atraso de energia?
Nick Sinas
@ Mohammed, você acha que a distorção dos sinais pode ser algum tipo de reverberação que pode ser limpa com a filtragem?
precisa
@NickS. Pode haver alguns truques de limpeza de reverberação por aí (não sei como isso seria realizado), mas juntei algo simples que será um estimador de energia, se você quiser dar uma olhada.
Spacey
0

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')
AshG
fonte