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

12

Nagrałem 2 sygnały z oscopa. Wyglądają tak: wprowadź opis zdjęcia tutaj

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: wprowadź opis zdjęcia tutaj

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ę:
wprowadź opis zdjęcia tutaj

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

Nick Sinas
źródło
Jak dokładnie wykonujesz korelację krzyżową? W odpowiedzi na twoje bezpośrednie pytanie nie powinieneś nic robić z sygnałami przed korelacją krzyżową, chociaż w niektórych przypadkach filtrowanie najpierw pomaga pozbyć się szumu, który może zniekształcić wyniki.
Jim Clay,
1
Proszę zamieścić użyty kod, a co ważniejsze, wykres sygnału korelacji krzyżowej. Niektóre narzędzia / biblioteki umieszczają wynik (lag = 0) na środku wykresu; Nie pamiętam, czy Matlab to robi.
pikenety
@pichenettes: zaktualizowany post
Nick Sinas,
@JimClay: zaktualizowany post
Nick Sinas,
@NickS. Jeśli twoje sygnały są idealnie wyrównane, uzyskasz pik pośrodku wykresu CC. Zatem szczyt przy 2000 oznacza brak opóźnienia. Powiedzmy teraz, że masz opóźnienie o 10 próbek, co oznacza, że ​​sygnał2 jest o 10 próbek wyłączony od sygnału1. To po prostu przeniesie twój szczyt w DW od 2000 do 2010 (lub 1990). Zatem opóźnienie czasowe odpowiada rzeczywistej lokalizacji szczytowej MINUS 2000.
Spacey

Odpowiedzi:

11

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

wprowadź opis zdjęcia tutaj

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:

wprowadź opis zdjęcia tutaj

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:

wprowadź opis zdjęcia tutaj

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:

wprowadź opis zdjęcia tutaj

To daje zatem próbek.δt=241

Mam nadzieję, że to pomogło.

Spacey
źródło
łał. Dziękuję bardzo. To są dla mnie nowe techniki, które zacznę badać. Czy jest jakiś sposób, by spojrzeć na użyty kod Matlab?
Nick Sinas,
Mam więc kroki 1 i 2 wykonane w Matlabie, a moje wyniki są zgodne z Twoimi, ale mam problemy z krokiem 3. Z jakich funkcji korzystałeś?
Nick Sinas,
@NickS. Zapytaj ... a otrzymasz, strzel do mnie e-mail, a ja mogę przesłać Ci kod, którego użyłem.
Spacey
@Mohammed Czy możesz podać swój kod, aby oszacować opóźnienie. Wysłałem ci e-mail w tej sprawie, więc proszę o pomoc
6

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.

Hilmar
źródło
4

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.

Jim Clay
źródło
Czyli sygnały nie są opóźnione?
Nick Sinas,
Nie jestem pewien - gdzie jest szczyt? Widzę, że jest blisko środka, ale nie jest jasne, że tak naprawdę jest na środku. Ponadto istnieje problem normalizacji mocy, który rozwiążę w edycji mojej odpowiedzi.
Jim Clay,
„Bezstronny” parametr zdecydowanie poprawia wygląd. więcej tego, czego bym się spodziewał. Będę się nad tym zastanawiał. Dzięki.
Nick Sinas,
@JimClay Być może Nick S koreluje obwiednie swoich sygnałów, a nie rzeczywistych sygnałów (Nick, czy to prawda?). To dałoby (z grubsza) ten trójkątny kształt, jaki sobie wyobrażam.
Spacey,
2
@NickS. Komentarz Mohammada sprawił, że spojrzałem i zdałem sobie sprawę, że masz ogromne przesunięcie DC, które psuje twoje wyniki. Odejmij średnią z obu sygnałów, a następnie uruchom na nich xcorr. Najpierw spróbowałbym tego bez opcji „obiektywnej”.
Jim Clay,
4

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)

wprowadź opis zdjęcia tutaj

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

wprowadź opis zdjęcia tutaj

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.

Jason R.
źródło
Jakieś pomysły na to, co jeszcze mogę zrobić? Może inna technika niż korelacja krzyżowa, a może jakiś filtr? Dzięki.
Nick Sinas,
@NickS. Ja też na to spojrzałem i nie są one opóźnionymi kopiami. Biorąc to pod uwagę, czy chcesz oszacować opóźnienie energii ? Myślę, że w tym przypadku byłoby to bardziej sensowne, opóźnienie VS sygnałów . Jeśli powiesz nam więcej o prowadzonym kanale / eksperymencie, możemy powiedzieć więcej o możliwych ścieżkach.
Spacey,
@Mohammad Thanks. Podstawowym kanałem jest stal. Wiesz, jak oszacować opóźnienie energii?
Nick Sinas,
@Mohammad, czy uważasz, że zniekształcenie sygnałów może być rodzajem pogłosu, który można by usunąć za pomocą filtrowania?
Nick Sinas,
@NickS. Mogą istnieć pewne sztuczki czyszczenia pogłosu (nie wiem, jak by się to udało), ale zebrałem razem coś prostego, co będzie estymatorem energii, jeśli chcesz rzucić okiem.
Spacey
0

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')
AshG
źródło