Zeitverzögerungsschätzung von Oszilloskopsignalen unter Verwendung von Kreuzkorrelation

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

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:

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:

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

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

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:

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:

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:

$\delta{t} = 241$

Ich hoffe das hat geholfen.

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.

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.

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)

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

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.

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