Szacowanie opóźnienia sygnałów oscyloskopowych z wykorzystaniem korelacji krzyżowej

Nagrałem 2 sygnały z oscopa. Wyglądają tak:

Chcę zmierzyć opóźnienie między nimi w Matlabie. Każdy sygnał ma 2000 próbek o częstotliwości próbkowania 2001000.5.

Dane znajdują się w pliku csv. To właśnie mam.
Usunąłem dane czasu z pliku csv, tak że tylko poziomy napięcia znajdują się w pliku csv.

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

To daje ten wynik:

Z tego, co przeczytałem, muszę wziąć korelację krzyżową tych sygnałów, co powinno dać mi szczyt związany z opóźnieniem czasowym. Jednak gdy wezmę korelację krzyżową tych sygnałów, otrzymam pik przy 2000, który, jak wiem, jest nieprawidłowy. Co powinienem zrobić z tymi sygnałami, zanim je skoreluję? Po prostu szukam jakiegoś kierunku.

EDYCJA: po usunięciu przesunięcia DC jest to wynik, który teraz otrzymuję:

Czy istnieje sposób, aby to wyczyścić, aby uzyskać bardziej określone opóźnienie czasowe?

EDYCJA 2: Oto pliki:
http://dl.dropbox.com/u/10147354/scope1col.csv
http://dl.dropbox.com/u/10147354/scope2col.csv

@NickS

Ponieważ nie jest pewne, czy drugi sygnał na wykresach jest w rzeczywistości jedynie opóźnioną wersją pierwszego, należy spróbować innych metod oprócz klasycznej korelacji krzyżowej. Wynika to z faktu, że korelacja krzyżowa (CC) jest jedynie estymatorem maksymalnego prawdopodobieństwa, jeśli twój sygnał (sygnały) są wzajemnie opóźnionymi wersjami. W tym przypadku najwyraźniej tak nie jest, nie mówiąc już o ich niestacjonarności.

W tym przypadku uważam, że to, co może zadziałać, to oszacowanie w czasie znaczącej energii sygnałów. To prawda, że ​​„znaczący” może lub nie może być nieco subiektywny, ale wierzę, że patrząc na twoje sygnały ze statystycznego punktu widzenia, będziemy w stanie określić ilościowo „znaczący” i stamtąd.

W tym celu wykonałem następujące czynności:

KROK 1: Oblicz obwiednie sygnału:

Ten krok jest prosty, ponieważ obliczana jest bezwzględna wartość wyjściowa transformaty Hilberta każdego z sygnałów. Istnieją inne metody obliczania kopert, ale jest to całkiem proste. Ta metoda zasadniczo oblicza analityczną formę twojego sygnału, innymi słowy reprezentację fazora. Kiedy bierzesz wartość bezwzględną, niszczysz fazę i tylko po energii.

Ponadto, ponieważ staramy się oszacować opóźnienie energii twoich sygnałów, takie podejście jest uzasadnione.

KROK 2: Usuwanie szumów za pomocą nieliniowych filtrów przyśrodkowych chroniących krawędzie:

To ważny krok. Celem jest wygładzenie kopert energetycznych, ale bez zniszczenia lub wygładzenia krawędzi i szybkiego wzrostu. W rzeczywistości poświęcono temu całe pole, ale do naszych celów możemy po prostu użyć łatwego do wdrożenia nieliniowego filtra Medial . (Filtrowanie mediany). Jest to potężna technika, ponieważ w przeciwieństwie do średniego filtrowania, filtrowanie przyśrodkowe nie niweluje twoich krawędzi, ale jednocześnie „wygładza” twój sygnał bez znaczącej degradacji ważnych krawędzi, ponieważ w żadnym momencie nie jest wykonywana żadna arytmetyka na twoim sygnale (pod warunkiem, że długość okna jest nieparzysta). W naszym przypadku wybrałem filtr przyśrodkowy o rozmiarze okna 25 próbek:

KROK 3: Usuń czas: Skonstruuj funkcje szacowania gęstości jądra Gaussa:

Co by się stało, gdybyś spojrzał na powyższą działkę bokiem zamiast normalnej drogi? Z matematycznego punktu widzenia oznacza to, co byś otrzymał, gdybyś rzucił każdą próbkę naszych denominowanych sygnałów na oś amplitudy y? W ten sposób uda nam się usunąć czas, że tak powiem, i będziemy w stanie studiować wyłącznie statystyki sygnałów.

Intuicyjnie co wyskakuje z powyższego rysunku? Chociaż energia hałasu jest niska, ma tę zaletę, że jest bardziej „popularna”. W przeciwieństwie do tego, podczas gdy obwiednia sygnału, która ma energię, jest bardziej energetyczna niż hałas, jest podzielona na progi. Co jeśli uznamy „popularność” za miarę energii? To właśnie zrobimy z (moją prymitywną) implementacją funkcji gęstości jądra (KDE) z jądrem gaussowskim.

Aby to zrobić, każda próbka jest pobierana, a funkcja gaussowska jest konstruowana z wykorzystaniem jej wartości jako średniej, a wstępnie ustawiona szerokość pasma (wariancja) jest wybierana a priori. Ustawienie wariancji gaussa jest ważnym parametrem, ale można go ustawić na podstawie statystyk hałasu na podstawie aplikacji i typowych sygnałów. (Mam tylko 2 twoje pliki do włączenia). Jeśli następnie zbudujemy oszacowanie KDE, otrzymamy następujący wykres:

Można myśleć o KDE jako ciągłej formie histogramu, a wariancję jako o szerokości bin. Ma jednak tę zaletę, że gwarantuje płynny plik PDF, na którym możemy następnie wykonać pierwszy i drugi rachunek różniczkowy. Teraz, gdy mamy gaussowskie KDE, możemy zobaczyć, gdzie próbki szumu osiągają największą popularność. Pamiętaj, że oś X reprezentuje tutaj rzuty naszych danych na przestrzeń amplitudy. W ten sposób możemy zobaczyć, w których progach hałas jest najbardziej „energetyczny”, a te mówią nam, jakich progów należy unikać.

Na drugim wykresie pobierana jest pierwsza pochodna gaussowskiego KDE i wybieramy odciętą pierwszą próbkę po pierwszej pochodnej po szczycie mieszaniny Gaussa, aby osiągnąć pewną wartość bliską zeru. (Lub pierwsze przejście przez zero). Możemy skorzystać z tej metody i być „bezpiecznym”, ponieważ nasze KDE zostało zbudowane z gładkich Gaussów o rozsądnej przepustowości, i wzięto pierwszą pochodną tej gładkiej i pozbawionej szumów funkcji. (Zwykle pierwsze pochodne mogą być problematyczne w niczym innym, jak w przypadku wysokich sygnałów SNR, ponieważ zwiększają szum).

Czarne linie pokazują wtedy, przy jakich progach rozsądnie byłoby „segmentować” obraz, aby uniknąć całej podłogi szumu. Jeśli następnie zastosujemy się do naszych oryginalnych sygnałów, uzyskamy następujące wykresy, z czarnymi liniami wskazującymi początek energii naszych sygnałów:

To daje zatem próbek.$\delta{t} = 241$

Mam nadzieję, że to pomogło.

Jest kilka problemów z robieniem tego z autokorelacją

1. Ogromne przesunięcie DC (już naprawione)
2. Okno czasu: xcorr () Matlaba ma irytującą konwencję polegającą na „zerowaniu padu” sygnału na obu końcach podczas przesuwania opóźnienia czasowego. To znaczy, okno danych jest funkcją opóźnienia czasowego. Stworzy to trójkątny kształt dla dowolnego sygnału stacjonarnego (w tym fal sinusoidalnych). Lepszym wyborem jest wybranie okna korelacji, tak aby rozmiar okna plus maksymalne opóźnienie pasowały do ​​całego okna danych, lub znormalizowanie korelacji krzyżowej na podstawie liczby próbek nie wypełnionych.
3. Te dwa sygnały nie wydają mi się szczególnie skorelowane. Kształt jest nieco podobny, ale specyficzne odstępy między pikami i spadkami są zupełnie inne, więc wątpię, aby nawet właściwa autokorelacja dała tu wiele wglądu.

Dużo prostszym podejściem byłoby użycie detektora progowego do znalezienia punktów początkowych i po prostu wykorzystanie różnicy między tymi punktami jako opóźnienia.

Jak wskazano w pikenetach, w tym przypadku pik na środku wyjścia wskazuje 0 opóźnienie. Przesunięcie piku od punktu środkowego to opóźnienie czasowe.

EDYCJA: Martwi mnie to, że korelacja jest tak idealnym trójkątem. To wskazuje mi, że korelacja krzyżowa nie powoduje normalizacji mocy. Daje to niesprawiedliwe nastawienie na mniejsze opóźnienia w stosunku do większych opóźnień. Zmodyfikowałbym twoje wywołanie xcorr na „cc = xcorr (x1, x2,„ bezstronne ”);”.

To nie jest idealne rozwiązanie, ponieważ wyniki dużego opóźnienia są teraz bardziej niestabilne niż wyniki niskiego opóźnienia, ponieważ są oparte na mniejszej ilości danych. Duży szczyt na końcach może być fałszywy z tego samego powodu, że możesz uzyskać 100% główek i żadnych ogonów na kilku rzutach monetą, podczas gdy jest bardzo mało prawdopodobne, aby zdarzyło się to przy wielu rzutach.

Jak zauważyli inni i wydaje się, że uświadomiłeś sobie na podstawie ostatniej edycji pytania, nie wydaje się, że korelacja krzyżowa da ci dobre oszacowanie opóźnienia czasowego dla pokazanych zestawów danych. Korelacja mierzy podobieństwo kształtu między dwoma szeregami czasowymi, przesuwając jeden po drugim dla szeregu opóźnień i obliczając iloczyn wewnętrzny między dwoma szeregami dla każdego opóźnienia. Wynik będzie miał dużą wielkość, gdy dwie serie będą jakościowo podobne lub „skorelowane” ze sobą. Jest to podobne do tego, jak iloczyn wewnętrzny dwóch wektorów jest największy, gdy dwa wektory są skierowane w tym samym kierunku.

Problem z danymi, które pokazałeś, polega na tym, że (przynajmniej w przypadku fragmentów, które widzimy) wydaje się, że nie ma bardzo podobnego kształtu. Nie ma żadnego opóźnienia, które można zastosować do jednego z sygnałów, aby wyglądał jak drugi, co dokładnie robisz, obliczając ich korelację krzyżową.

Są jednak przypadki, w których przydatna jest korelacja krzyżowa. Powiedz, że twój drugi sygnał był naprawdę przesuniętą w czasie wersją oryginału, nawet z dodanym dodatkowym szumem:

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)

Teraz nie jest od razu jasne, że oba sygnały są powiązane opóźnieniem czasowym. Jeśli jednak weźmiemy korelację krzyżową, otrzymamy:

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

który pokazuje pik przy prawidłowym opóźnieniu 200 próbek. Korelacja może być użytecznym narzędziem do określania opóźnienia czasowego w przypadku zestawów danych, które zawierają odpowiedni rodzaj podobieństwa.

W oparciu o sugestie Mahometa próbowałem napisać skrypt Matlab. Nie jestem jednak w stanie wywnioskować, czy konstruuje rozkład Gaussa na podstawie wariancji, a następnie przyjmuje oszacowanie KDE lub wykonuje oszacowanie KDE przy założeniu Gaussa.

Trudno też wywnioskować, w jaki sposób tłumaczy on czas przesunięcia KDE na dziedzinę czasu. Oto moja próba. Każdy użytkownik zainteresowany korzystaniem ze skryptu może bezpłatnie aktualizować ulepszoną wersję, jeśli to możliwe.

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