Estimation du temps de propagation des signaux de l'oscilloscope à l'aide de la corrélation croisée

J'ai enregistré 2 signaux d'un oscope. Ils ressemblent à ceci:

Je veux mesurer le délai entre eux dans Matlab. Chaque signal a 2000 échantillons avec une fréquence d'échantillonnage de 2001000,5.

Les données sont dans un fichier csv. C'est ce que j'ai jusqu'à présent.
J'ai effacé les données temporelles du fichier csv afin que seuls les niveaux de tension soient dans le fichier csv.

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

Cela donne ce résultat:

D'après ce que j'ai lu, je dois prendre la corrélation croisée de ces signaux et cela devrait me donner un pic relatif au retard. Cependant, quand je prends la corrélation croisée de ces signaux, j'obtiens un pic à 2000 que je sais n'est pas correct. Que dois-je faire pour ces signaux avant de les croiser? Je cherche juste une direction.

EDIT: après avoir supprimé le décalage DC, voici le résultat que j'obtiens maintenant:

Existe-t-il un moyen de nettoyer cela pour obtenir un délai plus défini?

EDIT 2: Voici les fichiers:
http://dl.dropbox.com/u/10147354/scope1col.csv
http://dl.dropbox.com/u/10147354/scope2col.csv

@NickS

Puisqu'il est loin d'être certain que le deuxième signal dans les tracés soit en fait une version uniquement retardée du premier, d'autres méthodes en plus de la corrélation croisée classique doivent être tentées. En effet, la corrélation croisée (CC) est simplement un estimateur du maximum de vraisemblance si vos signaux sont des versions retardées les uns des autres. Dans ce cas, ils ne le sont clairement pas, pour ne rien dire non plus de leur non-stationnarité.

Dans ce cas, je pense que ce qui peut fonctionner est une estimation du temps de l' énergie significative des signaux. Certes, «significatif» peut ou ne peut pas être quelque peu subjectif, mais je crois qu'en examinant vos signaux d'un point de vue statistique, nous serons en mesure de quantifier «significatif» et de partir de là.

À cette fin, j'ai fait ce qui suit:

ÉTAPE 1: Calculez les enveloppes de signal:

Cette étape est simple, car la valeur absolue de sortie de la transformée de Hilbert de chacun de vos signaux est calculée. Il existe d'autres méthodes pour calculer les enveloppes, mais c'est assez simple. Cette méthode calcule essentiellement la forme analytique de votre signal, en d'autres termes, la représentation des phaseurs. Lorsque vous prenez la valeur absolue, vous détruisez la phase et seulement après l'énergie.

De plus, puisque nous recherchons une estimation temporelle de l'énergie de vos signaux, cette approche est justifiée.

ÉTAPE 2: Supprime le bruit avec des filtres médiaux non linéaires préservant les bords:

Ceci est une étape importante. L'objectif ici est de lisser vos enveloppes énergétiques, mais sans destruction ni lissage de vos bords et des temps de montée rapides. Il y a en fait un champ entier consacré à cela, mais pour nos besoins ici, nous pouvons simplement utiliser un filtre Medial non linéaire facile à implémenter . (Filtrage médian). Il s'agit d'une technique puissante car contrairement au filtrage moyen , le filtrage médial n'annulera pas vos bords, mais en même temps `` lissera '' votre signal sans dégradation significative des bords importants, car à aucun moment aucune arithmétique n'est effectuée sur votre signal (à condition que la longueur de la fenêtre soit impaire). Pour notre cas ici, j'ai sélectionné un filtre médian d'échantillons de taille de fenêtre 25:

ÉTAPE 3: Suppression du temps: construction des fonctions d'estimation de la densité du noyau gaussien:

Que se passerait-il si vous regardiez l'intrigue ci-dessus de côté au lieu de la manière normale? Mathématiquement parlant, cela signifie, qu'obtiendriez-vous si vous projetiez chaque échantillon de nos signaux débruités sur l'axe d'amplitude y? Ce faisant, nous parviendrons à gagner du temps, pour ainsi dire, et à pouvoir étudier uniquement les statistiques du signal.

Intuitivement, qu'est-ce qui ressort de la figure ci-dessus? Bien que l'énergie sonore soit faible, elle a l'avantage d'être plus «populaire». En revanche, alors que l'enveloppe de signal qui a de l'énergie est plus énergétique que le bruit, elle est fragmentée entre les seuils. Et si nous considérions la «popularité» comme une mesure d'énergie? C'est ce que nous ferons avec (mon grossier) l'implémentation d'une fonction de densité de noyau , (KDE), avec un noyau gaussien.

Pour ce faire, chaque échantillon est prélevé et une fonction gaussienne construite en utilisant sa valeur comme moyenne, et une largeur de bande (variance) prédéfinie sélectionnée a priori. Le réglage de la variance de votre gaussien est un paramètre important, mais vous pouvez le définir en fonction des statistiques de bruit en fonction de votre application et des signaux typiques. (Je n'ai que vos 2 fichiers pour démarrer). Si nous construisons ensuite l'estimation de KDE, nous obtenons le graphique suivant:

Vous pouvez considérer le KDE comme une forme continue d'histogramme pour ainsi dire, et la variance comme la largeur de votre bac. Cependant, il a l'avantage de garantir un PDF fluide sur lequel nous pouvons ensuite effectuer des calculs de première et deuxième dérivées. Maintenant que nous avons les KDE gaussiens, nous pouvons voir où les échantillons de bruit atteignent un sommet de popularité. N'oubliez pas que l'axe des x représente ici les projections de nos données sur l'espace d'amplitude. Ainsi, nous pouvons voir dans quels seuils le bruit est le plus «énergétique» et ceux-ci nous indiquent quels seuils éviter.

Dans le deuxième graphique, la dérivée première des KDE gaussiens est prise, et nous choisissons les abscisses du premier échantillon après la dérivée première après le pic du mélange de gaussiens pour atteindre une certaine valeur proche de zéro. (Ou premier passage à zéro). Nous pouvons utiliser cette méthode et être «sûrs» car notre KDE a été construit avec des gaussiens lisses de bande passante raisonnable, et la première dérivée de cette fonction lisse et sans bruit a été prise. (Généralement, les dérivées premières peuvent être problématiques dans tout sauf les signaux SNR élevés car elles amplifient le bruit).

Les lignes noires montrent alors à quels seuils il serait sage de «segmenter» l'image, de manière à éviter tout le bruit de fond. Si nous appliquons ensuite à nos signaux d'origine, nous atteignons les graphiques suivants, avec les lignes noires indiquant le début de l'énergie de nos signaux:

$\delta{t} = 241$

J'espère que cela a aidé.

Il y a quelques problèmes avec l'autocorrélation

1. Énorme décalage CC (déjà fixé)
2. Fenêtre temporelle: xcorr () de Matlab a la convention gênante de "mettre à zéro" essentiellement le signal aux deux extrémités lorsque vous faites glisser le décalage temporel. C'est-à-dire que la fenêtre de données est fonction du décalage temporel. Cela créera une forme triangulaire pour tout signal stationnaire (y compris les ondes sinusoïdales). Les meilleurs choix sont de choisir une fenêtre de corrélation de sorte que la taille de la fenêtre plus le décalage temporel maximum s'inscrivent dans votre fenêtre de données totale, ou de normaliser la corrélation croisée par le nombre d'échantillons non remplis.
3. Les deux signaux ne me semblent pas particulièrement corrélés. La forme est quelque peu similaire, mais l'espacement spécifique des pics et des creux est assez différent, donc je doute que même une auto-corrélation appropriée produise beaucoup d'informations ici.

Une approche beaucoup plus simple serait d'utiliser un détecteur de seuil pour trouver les points de départ et d'utiliser simplement la différence entre ces points comme retard.

Comme l'indiquent les pichenettes, dans ce cas un pic au milieu de la sortie indique 0 retard. Le décalage du pic par rapport au point central correspond à votre décalage temporel.

EDIT: Cela me préoccupe que la corrélation est donc presque un triangle parfait. Cela m'indique que la corrélation croisée ne fait aucune normalisation de puissance. Cela donne un biais injuste aux petits décalages par rapport aux décalages plus importants. Je voudrais modifier votre appel xcorr à "cc = xcorr (x1, x2, 'impartial');".

Ce n'est pas, pensez-vous, une solution parfaite car les résultats à grand décalage sont désormais plus instables que les résultats à faible décalage car ils sont basés sur moins de données. Un grand pic aux extrémités peut être faux pour la même raison que vous pouvez obtenir 100% de têtes et pas de queue sur seulement quelques lancers de pièces, alors qu'il est très peu probable qu'il se produise sur de nombreux lancers.

Comme les autres l'ont souligné, et il semble que vous ayez réalisé sur la base de votre dernière modification de la question, il ne semble pas que la corrélation croisée vous donnera une bonne estimation du délai pour les ensembles de données affichés. La corrélation mesure la similitude de forme entre deux séries temporelles en glissant l'une sur l'autre pour une plage de décalages temporels et en calculant un produit interne entre les deux séries à chaque décalage. Le résultat aura une grande ampleur lorsque les deux séries sont qualitativement similaires ou qu'elles sont "corrélées" l'une avec l'autre. Ceci est similaire à la façon dont un produit intérieur de deux vecteurs est plus grand lorsque les deux vecteurs sont pointés dans la même direction.

Le problème avec les données que vous avez montrées est que (du moins pour les extraits que nous pouvons voir) il ne semble pas y avoir beaucoup de similitudes dans la forme. Il n'y a aucun délai que vous pouvez appliquer à l'un des signaux pour le faire ressembler à l'autre, c'est exactement ce que vous faites en calculant leur corrélation croisée.

Il existe cependant des cas où la corrélation croisée est utile. Supposons que votre deuxième signal était vraiment une version décalée dans le temps de l'original, même avec du bruit supplémentaire ajouté:

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)

Il n'est maintenant pas immédiatement clair que les deux signaux sont liés par un retard. Cependant, si nous prenons la corrélation croisée, nous obtenons:

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

qui montre un pic au décalage correct de 200 échantillons. La corrélation peut être un outil utile pour déterminer le délai, lorsqu'elle est appliquée à des ensembles de données qui contiennent le bon type de similitude.

Sur la base de la suggestion de Muhammad, j'ai essayé de faire un script Matlab. Cependant, je ne suis pas en mesure de déduire s'il construit une distribution gaussienne basée sur les variances et prend ensuite une estimation KDE ou s'il effectue une estimation KDE avec une hypothèse gaussienne.

Il est également difficile de déduire comment il traduit le temps de décalage KDE dans le domaine temporel. Voici ma tentative. Tout utilisateur qui souhaite utiliser le script est libre de mettre à jour la version améliorée si possible.

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