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

12

J'ai enregistré 2 signaux d'un oscope. Ils ressemblent à ceci: entrez la description de l'image ici

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: entrez la description de l'image ici

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:
entrez la description de l'image ici

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

Nick Sinas
la source
Comment faites-vous exactement la corrélation croisée? En réponse à votre question directe, vous ne devriez pas avoir besoin de faire quoi que ce soit à vos signaux avant la corrélation croisée, bien que dans certains cas, le filtrage aide d'abord à se débarrasser du bruit qui peut fausser les résultats.
Jim Clay
1
Veuillez poster le code que vous avez utilisé et surtout un tracé du signal de corrélation croisée. Certains outils / bibliothèques placent le score (lag = 0) au milieu du graphique; Je ne me souviens pas si Matlab fait ça.
pichenettes
@pichenettes: article mis à jour
Nick Sinas
@JimClay: article mis à jour
Nick Sinas
@NickS. Si vos signaux sont parfaitement alignés, vous obtiendrez un pic au milieu de votre tracé cc. Un pic à 2000 signifie donc pas de retard. Supposons maintenant que vous ayez un retard de 10 échantillons, ce qui signifie que signal2 est à 10 échantillons de signal1. Cela ne fera que déplacer votre pic dans le cc de 2000 à 2010 (ou 1990). Votre délai correspond donc à votre emplacement de pointe réel, MOINS 2000.
Spacey

Réponses:

11

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

entrez la description de l'image ici

É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:

entrez la description de l'image ici

É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:

entrez la description de l'image ici

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:

entrez la description de l'image ici

δt=241

J'espère que cela a aidé.

Spacey
la source
sensationnel. Merci beaucoup. Ce sont toutes de nouvelles techniques pour moi que je vais commencer à rechercher. Existe-t-il un moyen de jeter un œil au code matlab que vous avez utilisé?
Nick Sinas
J'ai donc fait les étapes # 1 et # 2 dans Matlab, et mes résultats correspondent aux vôtres, mais j'ai des problèmes avec l'étape # 3. Quelles fonctions avez-vous utilisées?
Nick Sinas
@NickS. Demandez, .. et vous recevrez, envoyez-moi un e-mail et je pourrai vous envoyer le code que j'ai utilisé.
Spacey
@Mohammed Pourriez-vous s'il vous plaît poster votre code pour estimer le délai. Je vous ai envoyé un e-mail à ce sujet, alors veuillez aider
6

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.

Hilmar
la source
4

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.

Jim Clay
la source
Cela signifie que les signaux ne sont pas retardés?
Nick Sinas
Je ne sais pas - où est le pic? Je peux voir qu'il est près du milieu, mais il n'est pas clair qu'il se trouve réellement au milieu. En outre, il y a un problème de normalisation de l'alimentation que j'aborderai dans une modification de ma réponse.
Jim Clay
Le paramètre «impartial» le rend définitivement meilleur. plus de ce à quoi je m'attendrais. Je continuerai de regarder cela. Merci.
Nick Sinas
@ JimClay Peut-être que Nick S corrèle les enveloppes de ses signaux et non les signaux réels (Nick est-ce vrai?). Cela donnerait (à peu près) cette forme triangulaire, j'imagine.
Spacey
2
@NickS. Le commentaire de Mohammad m'a fait regarder et réaliser que vous avez un énorme décalage DC qui gâche vos résultats. Soustrayez la moyenne de vos deux signaux, puis exécutez xcorr sur eux. Je l'essayerais sans l'option "impartiale" en premier.
Jim Clay
4

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)

entrez la description de l'image ici

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

entrez la description de l'image ici

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.

Jason R
la source
Des idées pour quoi d'autre je pourrais faire? Peut-être une autre technique que la corrélation croisée ou peut-être un type de filtre? Merci.
Nick Sinas
@NickS. Je l'ai également regardé et ce ne sont pas des copies différées les unes des autres. Cela étant dit, voulez-vous estimer le retard de l' énergie ? Je pense que cela aurait plus de sens dans ce cas, VS retard des signaux . Si vous nous en dites plus sur le canal / expérience sous-jacent que vous exécutez, nous pouvons vous en dire plus sur les chemins possibles.
Spacey
@Mohammad Merci. Le canal sous-jacent est en acier. Une idée de comment estimer le retard de l'énergie?
Nick Sinas
@Mohammad pensez-vous que la distorsion des signaux pourrait être une sorte de réverbération qui pourrait être nettoyée avec un filtrage?
Nick Sinas
@NickS. Il peut y avoir des astuces de nettoyage de réverbération (je ne sais pas comment elles seraient accomplies), mais j'ai concocté quelque chose de simple qui sera un estimateur d'énergie si vous voulez y jeter un œil.
Spacey
0

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')
AshG
la source