Zeitverzögerungsschätzung von Oszilloskopsignalen unter Verwendung von Kreuzkorrelation

12

Ich habe 2 Signale von einem oscope aufgezeichnet. Sie sehen so aus: Bildbeschreibung hier eingeben

Ich möchte die Zeitverzögerung zwischen ihnen in Matlab messen. Jedes Signal hat 2000 Abtastwerte mit einer Abtastfrequenz von 2001000.5.

Die Daten befinden sich in einer CSV-Datei. Das habe ich bisher.
Ich habe die Zeitdaten aus der CSV-Datei gelöscht, sodass nur die Spannungspegel in der CSV-Datei enthalten sind.

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

Dies ergibt folgendes Ergebnis: Bildbeschreibung hier eingeben

Nach dem, was ich gelesen habe, muss ich die Kreuzkorrelation dieser Signale nehmen und dies sollte mir einen Peak in Bezug auf die Zeitverzögerung geben. Wenn ich jedoch die Kreuzkorrelation dieser Signale nehme, erhalte ich einen Peak bei 2000, von dem ich weiß, dass er nicht korrekt ist. Was soll ich mit diesen Signalen tun, bevor ich sie kreuzweise korreliere? Ich suche nur nach einer Richtung.

BEARBEITEN: Nach dem Entfernen des DC-Offsets ist dies das Ergebnis, das ich jetzt erhalte:
Bildbeschreibung hier eingeben

Gibt es eine Möglichkeit, dies zu bereinigen, um eine definierte Verzögerung zu erhalten?

EDIT 2: Hier sind die Dateien:
http://dl.dropbox.com/u/10147354/scope1col.csv
http://dl.dropbox.com/u/10147354/scope2col.csv

Nick Sinas
quelle
Wie genau machen Sie die Kreuzkorrelation? Bei der Beantwortung Ihrer direkten Frage sollten Sie vor der Kreuzkorrelation nichts mit Ihren Signalen zu tun haben. In einigen Fällen hilft das Filtern jedoch, Rauschen zu beseitigen, das die Ergebnisse verfälschen kann.
Jim Clay
1
Bitte geben Sie den von Ihnen verwendeten Code und vor allem eine grafische Darstellung des Kreuzkorrelationssignals an. Einige Tools / Bibliotheken setzen den (lag = 0) Score in die Mitte des Graphen. Ich erinnere mich nicht, ob Matlab das tut.
Pichenettes
@pichenettes: aktualisierter Beitrag
Nick Sinas
@ JimClay: aktualisierter Beitrag
Nick Sinas
@NickS. Wenn Ihre Signale perfekt ausgerichtet sind, erhalten Sie eine Spitze in der Mitte Ihres CC-Diagramms. Ein Peak bei 2000 bedeutet also keine Verzögerung. Nehmen wir nun an, Sie haben eine Verzögerung von 10 Samples, was bedeutet, dass Signal2 10 Samples von Signal1 entfernt ist. Dies verschiebt nur Ihren Höchststand im cc von 2000 bis 2010 (oder 1990). Ihre Zeitverzögerung entspricht also Ihrem tatsächlichen Höchststandort, MINUS 2000.
Spacey

Antworten:

11

@NickS

Da es sich bei dem zweiten Signal in den Diagrammen bei weitem nicht um eine ausschließlich verzögerte Version des ersten Signals handelt, müssen neben der klassischen Kreuzkorrelation andere Methoden versucht werden. Dies liegt daran, dass die Kreuzkorrelation (CC) nur ein Maximum-Likelihood-Schätzer ist, wenn Ihre Signale verzögerte Versionen voneinander sind. In diesem Fall sind sie eindeutig nicht, ganz zu schweigen von ihrer Nichtstationarität.

In diesem Fall glaube ich, dass eine Zeitschätzung der signifikanten Energie der Signale funktionieren könnte . Zugegeben, "signifikant" kann oder kann nicht subjektiv sein, aber ich glaube, wenn wir Ihre Signale von einem statistischen Standpunkt aus betrachten, können wir "signifikant" quantifizieren und von dort aus fortfahren.

Zu diesem Zweck habe ich Folgendes getan:

SCHRITT 1: Berechnen Sie die Signalhüllkurven:

Dieser Schritt ist einfach, da der absolute Ausgabewert der Hilbert-Transformation jedes Ihrer Signale berechnet wird. Es gibt andere Methoden zum Berechnen von Umschlägen, dies ist jedoch recht einfach. Diese Methode berechnet im Wesentlichen die analytische Form Ihres Signals, dh die Zeigerdarstellung. Wenn Sie den absoluten Wert nehmen, zerstören Sie Phase und nur nach Energie.

Da wir außerdem eine Zeitverzögerungsschätzung der Energie Ihrer Signale durchführen, ist dieser Ansatz gerechtfertigt.

Bildbeschreibung hier eingeben

SCHRITT 2: Rauschunterdrückung mit kantenerhaltenden nichtlinearen Medialfiltern:

Dies ist ein wichtiger Schritt. Das Ziel hierbei ist es, Ihre Energiehüllen zu glätten, ohne dabei die Kanten und die schnellen Anstiegszeiten zu zerstören oder zu glätten. Es gibt tatsächlich ein ganzes Feld, das diesem Thema gewidmet ist, aber für unsere Zwecke können wir hier einfach einen einfach zu implementierenden nichtlinearen Medial-Filter verwenden . (Medianfilterung). Dies ist eine leistungsfähige Technik , denn im Gegensatz zu mittlerer Filterung, medial - Filterung wird die Kanten nicht null, aber zugleich ‚glatt‘ out Ihr Signal ohne signifikante Verschlechterung der wichtigen Kanten, zu keiner Zeit , da jedes arithmetisches Wesen auf dem Signal durchgeführt (vorausgesetzt die Fensterlänge ist ungerade). Für unseren Fall hier habe ich einen Medialfilter mit einer Fenstergröße von 25 Samples ausgewählt:

Bildbeschreibung hier eingeben

SCHRITT 3: Entfernen der Zeit: Konstruieren der Dichteschätzfunktionen des Gaußschen Kernels:

Was würde passieren, wenn Sie die obige Handlung von der Seite betrachten, anstatt wie üblich? Mathematisch gesehen bedeutet das, was würden Sie erhalten, wenn Sie jedes Sample unserer entrauschten Signale auf die y-Amplitudenachse projizieren würden? Auf diese Weise schaffen wir es, sozusagen Zeit zu sparen und nur die Signalstatistik zu studieren.

Intuitiv, was aus der obigen Abbildung hervorgeht? Die Geräuschenergie ist zwar gering, hat aber den Vorteil, dass sie „beliebter“ ist. Im Gegensatz dazu ist die Signalhüllkurve, die Energie enthält, energiereicher als das Rauschen, sie ist jedoch über Schwellenwerte hinweg fragmentiert. Was wäre, wenn wir die Popularität als Maß für die Energie betrachten würden? Dies ist, was wir mit (meiner groben) Implementierung einer Kerneldichtefunktion (KDE) mit einem Gaußschen Kernel tun werden .

Zu diesem Zweck wird jeder Abtastwert entnommen und eine Gauß-Funktion unter Verwendung ihres Wertes als Mittelwert und einer voreingestellten Bandbreite (Varianz), die von vornherein ausgewählt wurde, konstruiert. Das Einstellen der Varianz Ihres Gaußschen ist ein wichtiger Parameter. Sie können sie jedoch basierend auf der Rauschstatistik basierend auf Ihrer Anwendung und typischen Signalen einstellen. (Ich habe nur deine 2 Dateien, auf die ich mich begeben kann). Wenn wir dann die KDE-Schätzung erstellen, erhalten wir das folgende Diagramm:

Bildbeschreibung hier eingeben

Sie können sich KDE sozusagen als kontinuierliche Form eines Histogramms und die Varianz als Ihre Bin-Breite vorstellen. Es hat jedoch den Vorteil, dass ein reibungsloses PDF gewährleistet ist, mit dem wir dann eine erste und eine zweite Ableitungsberechnung durchführen können. Jetzt, da wir die Gaußschen KDEs haben, können wir sehen, wo die Rauschproben an Popularität gewinnen. Denken Sie daran, dass die X-Achse hier die Projektionen unserer Daten auf den Amplitudenraum darstellt. Auf diese Weise können wir sehen, bei welchen Schwellenwerten das Rauschen am energiereichsten ist und welche Schwellenwerte zu vermeiden sind.

In der zweiten Darstellung wird die erste Ableitung der Gaußschen KDEs genommen, und wir wählen die Abszisse der ersten Probe nach der ersten Ableitung nach dem Peak der Mischung der Gaußschen, um einen bestimmten Wert nahe Null zu erreichen. (Oder erster Nulldurchgang). Wir können diese Methode anwenden und sind "sicher", da unser KDE aus glatten Gaußschen mit angemessener Bandbreite aufgebaut ist und die erste Ableitung dieser glatten und rauschfreien Funktion verwendet wurde. (In der Regel können erste Ableitungen bei Signalen mit hohem SNR problematisch sein, da sie das Rauschen verstärken.)

Die schwarzen Linien zeigen dann an, bei welchen Schwellenwerten wir das Bild "segmentieren" sollten, so dass wir das gesamte Grundrauschen vermeiden. Wenn wir dann auf unsere ursprünglichen Signale anwenden, erhalten wir die folgenden Diagramme, wobei die schwarzen Linien den Beginn der Energie unserer Signale angeben:

Bildbeschreibung hier eingeben

δt=241

Ich hoffe das hat geholfen.

Spacey
quelle
Beeindruckend. Ich danke dir sehr. Für mich sind das alles neue Techniken, mit deren Erforschung ich beginnen werde. Kann ich mir den Matlab-Code ansehen, den Sie verwendet haben?
Nick Sinas
Ich habe also die Schritte 1 und 2 in Matlab ausgeführt und meine Ergebnisse stimmen mit Ihren überein, aber ich habe Probleme mit Schritt 3. Welche Funktionen haben Sie genutzt?
Nick Sinas
@NickS. Fragen Sie, .. und Sie erhalten, schreiben Sie mir eine E-Mail und ich kann Ihnen den von mir verwendeten Code senden.
Spacey
@ Mohammed Könnten Sie bitte Ihren Code posten, um die Zeitverzögerung zu schätzen. Ich habe Ihnen diesbezüglich eine E-Mail gesendet. Bitte helfen Sie
6

Bei der Autokorrelation treten einige Probleme auf

  1. Enormer DC-Offset (bereits behoben)
  2. Zeitfenster: Matlabs xcorr () hat die lästige Konvention, das Signal an beiden Enden im Wesentlichen auf Null zu setzen, wenn Sie die Zeitverzögerung verschieben. Dh das Datenfenster ist eine Funktion der Zeitverzögerung. Dadurch wird für alle stationären Signale (einschließlich Sinuswellen) eine dreieckige Form erstellt. Besser ist es, ein Korrelationsfenster so zu wählen, dass die Fenstergröße und die maximale Zeitverzögerung in Ihr gesamtes Datenfenster passen, oder die Kreuzkorrelation durch die Anzahl der nicht aufgefüllten Stichproben zu normalisieren.
  3. Die beiden Signale sehen für mich nicht besonders korreliert aus. Die Form ist etwas ähnlich, aber der spezifische Abstand von Peaks und Dips ist sehr unterschiedlich, so dass ich bezweifle, dass selbst eine korrekte Autokorrelation hier viel Einsicht bringen würde.

Ein viel einfacherer Ansatz wäre die Verwendung eines Schwellenwertdetektors, um die Startpunkte zu finden und einfach die Differenz zwischen diesen Punkten als Verzögerung zu verwenden.

Hilmar
quelle
4

In diesem Fall zeigt eine Spitze in der Mitte der Ausgabe, wie Pichenetten gezeigt haben, eine Verzögerung von 0 an. Der Abstand des Peaks vom Mittelpunkt ist Ihre Zeitverzögerung.

EDIT: Es geht mir darum, dass die Korrelation so fast ein perfektes Dreieck ist. Das zeigt mir, dass die Kreuzkorrelation keine Leistungsnormalisierung durchführt. Dies gibt kleineren Verzögerungen eine unfaire Tendenz gegenüber größeren Verzögerungen. Ich würde Ihren xcorr-Aufruf zu "cc = xcorr (x1, x2, 'unvoreingenommen');" ändern.

Dies ist jedoch keine perfekte Lösung, da die Ergebnisse mit großen Verzögerungen jetzt instabiler sind als die Ergebnisse mit geringen Verzögerungen, da sie auf weniger Daten basieren. Ein großer Peak an den Extremitäten kann aus dem gleichen Grund falsch sein, dass Sie bei nur wenigen Münzwürfen 100% Köpfe und keine Schwänze erhalten können, während es bei vielen Würfen äußerst unwahrscheinlich ist, dass dies geschieht.

Jim Clay
quelle
Das heißt, dass die Signale nicht verzögert werden?
Nick Sinas
Ich bin nicht sicher, wo ist der Höhepunkt? Ich kann sehen, dass es nahe der Mitte ist, aber es ist nicht klar, dass es tatsächlich in der Mitte ist. Außerdem gibt es ein Problem mit der Energienormalisierung, auf das ich in einer Bearbeitung meiner Antwort eingehen werde.
Jim Clay
'unvoreingenommene' Parameter lässt es definitiv besser aussehen. mehr von dem, was ich erwarten würde. Ich werde das weiter untersuchen. Vielen Dank.
Nick Sinas
@ JimClay Vielleicht korreliert Nick S die Hüllkurven seiner Signale und nicht die tatsächlichen Signale (Nick ist das wahr?). Dies würde (ungefähr) diese dreieckige Form ergeben, die ich mir vorstelle.
Spacey
2
@NickS. Durch Mohammeds Kommentar wurde mir klar, dass Sie einen riesigen Gleichstromversatz haben, der Ihre Ergebnisse durcheinander bringt. Subtrahieren Sie den Mittelwert von beiden Signalen und führen Sie dann xcorr für sie aus. Ich würde es zuerst ohne die "unvoreingenommene" Option versuchen.
Jim Clay
4

Wie die anderen bereits betont haben und Sie anscheinend aufgrund Ihrer letzten Bearbeitung der Frage festgestellt haben, lässt sich die Zeitverzögerung für die angezeigten Datensätze nicht gut durch Kreuzkorrelation abschätzen. Die Korrelation misst die Ähnlichkeit in der Form zwischen zwei Zeitreihen, indem sie für eine Reihe von Zeitverzögerungen übereinander gleitet und bei jeder Verzögerung ein inneres Produkt zwischen den beiden Reihen berechnet. Das Ergebnis wird eine große Größe haben, wenn die beiden Reihen qualitativ ähnlich sind oder miteinander "korreliert" sind. Dies ist vergleichbar damit, wie ein inneres Produkt zweier Vektoren am größten ist, wenn die beiden Vektoren in dieselbe Richtung weisen.

Das Problem mit den Daten, die Sie gezeigt haben, ist, dass (zumindest für die Schnipsel, die wir sehen können) die Form nicht sehr ähnlich zu sein scheint. Es gibt keine Verzögerung, die Sie auf eines der Signale anwenden können, damit es wie das andere aussieht. Genau das tun Sie, indem Sie deren Kreuzkorrelation berechnen.

Es gibt jedoch Fälle, in denen eine Kreuzkorrelation sinnvoll ist. Angenommen, Ihr zweites Signal war wirklich eine zeitversetzte Version des Originals, auch wenn ein zusätzliches Rauschen hinzugefügt wurde:

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)

Bildbeschreibung hier eingeben

Es ist jetzt nicht sofort klar, dass die beiden Signale zeitverzögert zusammenhängen. Wenn wir jedoch die Kreuzkorrelation nehmen, erhalten wir:

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

Bildbeschreibung hier eingeben

Dies zeigt einen Peak mit der korrekten Verzögerung von 200 Samples. Korrelation kann ein nützliches Werkzeug zur Bestimmung der Zeitverzögerung sein, wenn sie auf Datensätze angewendet wird, die den richtigen Ähnlichkeitstyp enthalten.

Jason R
quelle
Irgendwelche Ideen, was ich sonst noch tun könnte? Vielleicht eine andere Technik als Kreuzkorrelation oder vielleicht eine Art Filter? Vielen Dank.
Nick Sinas
@NickS. Ich habe es mir auch angeschaut und es handelt sich nicht um verspätete Kopien voneinander. Möchten Sie die Verzögerung der Energie einschätzen ? Ich denke das wäre in diesem Fall sinnvoller, VS Verzögerung der Signale . Wenn Sie uns mehr über den zugrunde liegenden Kanal / das zugrunde liegende Experiment erzählen, können wir Ihnen mehr über mögliche Pfade erzählen.
Spacey
@ Mohammad Danke. Der darunter liegende Kanal ist Stahl. Irgendeine Idee, wie man die Verzögerung der Energie abschätzt?
Nick Sinas
@Mohammad glaubst du, dass die Verzerrung der Signale eine Art Nachhall sein könnte, der durch Filtern beseitigt werden könnte?
Nick Sinas
@NickS. Möglicherweise gibt es einige Tricks zur Reinigung des Nachhalls (ich weiß nicht, wie diese ausgeführt werden sollen), aber ich habe etwas Einfaches zusammengestellt, das eine Energiekalkulation darstellt, wenn Sie einen Blick darauf werfen möchten.
Spacey
0

Auf Mohammeds Vorschlag hin habe ich versucht, ein Matlab-Skript zu erstellen. Ich kann jedoch nicht ableiten, ob er eine Gaußsche Verteilung auf der Grundlage der Varianzen erstellt und dann eine KDE-Schätzung vornimmt oder eine KDE-Schätzung mit Gaußscher Annahme vornimmt.

Es ist auch schwierig abzuleiten, wie er die KDE-Versatzzeit in den Zeitbereich übersetzt. Hier ist mein Versuch dazu. Jeder Benutzer, der an der Verwendung des Skripts interessiert ist, kann die verbesserte Version nach Möglichkeit kostenlos aktualisieren.

%% 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
quelle