Estimasi waktu tunda sinyal osiloskop menggunakan korelasi silang

Saya telah merekam 2 sinyal dari oscope. Mereka terlihat seperti ini:

Saya ingin mengukur waktu tunda di antara mereka di Matlab. Setiap sinyal memiliki 2000 sampel dengan frekuensi sampling 2001000.5.

Data dalam file csv. Inilah yang saya miliki sejauh ini.
Saya menghapus data waktu dari file csv sehingga hanya level tegangan yang ada di file csv.

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

Ini memberikan hasil ini:

Dari apa yang saya baca saya perlu mengambil korelasi silang dari sinyal-sinyal ini dan ini akan memberi saya puncak terkait dengan penundaan waktu. Namun ketika saya mengambil korelasi silang dari sinyal-sinyal ini saya mendapatkan puncak pada 2000 yang saya tahu tidak benar. Apa yang harus saya lakukan terhadap sinyal-sinyal ini sebelum saya mengkorelasikannya? Hanya mencari arah.

EDIT: setelah menghapus DC offset, inilah hasil yang sekarang saya dapatkan:

Apakah ada cara untuk membersihkan ini untuk mendapatkan waktu tunda yang lebih jelas?

EDIT 2: Berikut adalah file-nya:
http://dl.dropbox.com/u/10147354/scope1col.csv
http://dl.dropbox.com/u/10147354/scope2col.csv

@NickS

Karena jauh dari pasti bahwa sinyal kedua dalam plot sebenarnya adalah semata-mata versi tertunda dari yang pertama, metode lain selain klasik korelasi silang harus dicoba. Ini karena korelasi silang (CC) hanyalah penaksir kemungkinan maksimum jika sinyal Anda adalah versi yang tertunda satu sama lain. Dalam hal ini, mereka jelas tidak mengatakan apa-apa tentang ketidakstabilan mereka juga.

Dalam hal ini, saya percaya apa yang mungkin berhasil adalah estimasi waktu dari energi signifikan dari sinyal. Memang, 'signifikan' dapat atau tidak bisa agak subyektif, tetapi saya percaya bahwa dengan melihat sinyal Anda dari sudut pandang statistik, kami akan dapat mengukur 'signifikan' dan pergi dari sana.

Untuk tujuan ini, saya melakukan yang berikut:

LANGKAH 1: Hitung amplop sinyal:

Langkah ini sederhana, karena nilai absolut output Hilbert-Transform dari masing-masing sinyal Anda dihitung. Ada metode lain untuk menghitung amplop, tetapi ini cukup mudah. Metode ini pada dasarnya menghitung bentuk analitik sinyal Anda, dengan kata lain, representasi fasor. Ketika Anda mengambil nilai absolut, Anda menghancurkan fase dan hanya setelah energi.

Selain itu karena kami mengejar estimasi waktu tunda energi sinyal Anda, pendekatan ini diperlukan.

LANGKAH 2: De-noise dengan Filter Medial non-linear yang mempertahankan tepi:

Ini merupakan langkah penting. Tujuannya di sini adalah untuk menghaluskan amplop energi Anda, tetapi tanpa merusak atau menghaluskan tepi dan waktu kenaikan cepat Anda. Sebenarnya ada seluruh bidang yang dikhususkan untuk ini, tetapi untuk tujuan kita di sini, kita cukup menggunakan filter Medial non-linear yang mudah diterapkan . (Penyaringan Median). Ini adalah teknik yang ampuh karena tidak seperti pemfilteran rata-rata , pemfilteran medial tidak akan menghilangkan tepi Anda, tetapi pada saat yang sama 'menghaluskan' sinyal Anda tanpa degradasi signifikan pada tepi-tepi penting, karena pada waktu tidak ada aritmatika yang dilakukan pada sinyal Anda (asalkan panjang jendela aneh). Untuk kasus kami di sini, saya memilih filter medial ukuran jendela 25 sampel:

LANGKAH 3: Hapus Waktu: Bangun Fungsi Estimasi Kerapatan Kernel Gaussian:

Apa yang akan terjadi jika Anda melihat plot di atas daripada cara yang biasa? Secara matematis, itu berarti, apa yang akan Anda dapatkan jika Anda memproyeksikan setiap sampel sinyal denoised kami ke sumbu y-amplitudo? Dalam melakukan ini, kami akan mengatur untuk menghilangkan waktu untuk berbicara, dan dapat mempelajari statistik sinyal semata-mata.

Secara intuitif apa yang muncul dari gambar di atas? Walaupun energi kebisingan rendah, ia memiliki keuntungan karena lebih 'populer'. Sebaliknya, sementara amplop sinyal yang memiliki energi lebih energik daripada noise, itu terfragmentasi melintasi ambang batas. Bagaimana jika kita menganggap 'popularitas' sebagai ukuran energi? Inilah yang akan kita lakukan dengan implementasi (Kasar saya) dari Fungsi Kernel Density , (KDE), dengan Kernel Gaussian.

Untuk melakukan ini, setiap sampel diambil dan fungsi gaussian dibangun menggunakan nilainya sebagai rata-rata, dan bandwidth yang ditetapkan sebelumnya (varians) dipilih a-priori. Mengatur varian gaussian Anda adalah parameter penting, tetapi Anda dapat mengaturnya berdasarkan statistik kebisingan berdasarkan aplikasi Anda dan sinyal khas. (Saya hanya ingin mengaktifkan 2 file Anda). Jika kita membangun Estimasi KDE, kita mendapatkan plot berikut:

Anda dapat menganggap KDE sebagai bentuk histogram terus-menerus sehingga bisa dikatakan, dan varians sebagai lebar bin Anda. Namun memiliki keuntungan menjamin kelancaran PDF sehingga kita dapat melakukan kalkulus derivitave pertama dan kedua. Sekarang kita memiliki Gaussian KDE, kita dapat melihat di mana sampel suara memuncak dalam popularitas. Ingatlah bahwa sumbu x di sini mewakili proyeksi data kami ke ruang amplitudo. Dengan demikian, kita dapat melihat ambang mana kebisingan yang paling 'energik' dalam, dan mereka memberitahu kita ambang yang harus dihindari.

Dalam plot kedua, turunan pertama dari KDE Gaussian diambil, dan kami mengambil absis dari sampel pertama setelah turunan pertama setelah puncak campuran Gaussians untuk mencapai nilai tertentu mendekati nol. (Atau zero-crossing pertama). Kita dapat menggunakan metode ini dan menjadi 'aman' karena KDE kami dibuat dari Gaussians dengan bandwidth yang masuk akal, dan turunan pertama dari fungsi halus dan tanpa noise ini diambil. (Biasanya turunan pertama bisa bermasalah dalam apa pun kecuali sinyal SNR tinggi karena mereka memperbesar kebisingan).

Garis hitam menunjukkan kemudian pada ambang berapa kita akan bijaksana untuk 'segmen' gambar di, sehingga kita menghindari seluruh lantai kebisingan. Jika kemudian kita menerapkan sinyal asli kita, kita mendapatkan plot berikut, dengan garis hitam menunjukkan awal energi sinyal kita:

Ini dengan demikian menghasilkan sampel.$\delta{t} = 241$

Saya harap ini membantu.

Ada beberapa masalah dalam melakukan ini dengan autokorelasi

1. Offset DC besar (sudah diperbaiki)
2. Jendela waktu: Matlab's xcorr () memiliki konvensi yang mengganggu pada dasarnya "nol pad" sinyal di kedua ujung saat Anda menggeser jeda waktu. Yaitu jendela data adalah fungsi dari jeda waktu. Ini akan membuat bentuk segitiga untuk sinyal diam apa pun (termasuk gelombang sinus). Pilihan yang lebih baik adalah memilih jendela korelasi sehingga ukuran jendela plus jeda waktu maksimum masuk ke dalam jendela data total Anda, atau untuk menormalkan korelasi silang dengan jumlah sampel yang tidak berlapis.
3. Kedua sinyal tidak terlihat berkorelasi khusus dengan saya. Bentuknya agak mirip tetapi jarak spesifik puncak dan kemiringan sangat berbeda, jadi saya ragu bahwa bahkan korelasi otomatis yang baik akan menghasilkan banyak wawasan di sini.

Pendekatan yang jauh lebih sederhana adalah dengan menggunakan detektor ambang untuk menemukan titik awal dan hanya menggunakan perbedaan antara titik-titik ini sebagai penundaan.

Seperti yang ditunjukkan pichenettes, dalam hal ini puncak di tengah output menunjukkan 0 lag. Offset puncak dari titik tengah adalah jeda waktu Anda.

EDIT: Ini mengkhawatirkan saya bahwa korelasinya hampir merupakan segitiga sempurna. Itu menunjukkan kepada saya bahwa korelasi silang tidak melakukan normalisasi daya. Itu memberikan bias tidak adil untuk kelambatan yang lebih kecil dari kelambatan yang lebih besar. Saya akan memodifikasi panggilan xcorr Anda menjadi "cc = xcorr (x1, x2, 'tidak bias');".

Itu bukan, pikiran Anda, solusi sempurna karena hasil jeda besar sekarang lebih tidak stabil daripada hasil jeda rendah karena mereka didasarkan pada data yang lebih sedikit. Puncak besar di ekstremitas bisa palsu karena alasan yang sama bahwa Anda bisa mendapatkan 100% kepala dan tidak ada ekor hanya dengan beberapa kali lemparan koin, sementara itu sangat tidak mungkin terjadi pada banyak lemparan.

Seperti yang telah ditunjukkan oleh yang lain, dan tampaknya Anda telah menyadari berdasarkan edit terakhir Anda pada pertanyaan, tampaknya tidak ada korelasi silang yang akan memberi Anda perkiraan waktu tunda yang baik untuk dataset yang ditampilkan. Korelasi mengukur kesamaan bentuk antara dua deret waktu dengan menggeser yang satu ke yang lainnya untuk rentang jeda waktu dan menghitung produk dalam antara kedua seri pada setiap jeda. Hasilnya akan memiliki besaran yang besar ketika kedua seri secara kualitatif serupa, atau mereka "berkorelasi" satu sama lain. Ini mirip dengan bagaimana produk dalam dua vektor terbesar ketika dua vektor menunjuk ke arah yang sama.

Masalah dengan data yang Anda tunjukkan adalah (setidaknya untuk cuplikan yang dapat kita lihat) tampaknya tidak ada banyak kesamaan bentuk. Tidak ada penundaan yang dapat Anda terapkan pada salah satu sinyal untuk membuatnya terlihat seperti yang lain, yang persis sama dengan apa yang Anda lakukan dengan menghitung korelasi-silangnya.

Namun, ada beberapa contoh di mana korelasi silang bermanfaat. Katakan bahwa sinyal kedua Anda benar-benar merupakan versi asli dari time-shifted, bahkan dengan tambahan noise:

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)

Sekarang tidak segera jelas bahwa kedua sinyal terkait oleh penundaan waktu. Namun, jika kita mengambil korelasi silang, kita mendapatkan:

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

yang menunjukkan puncak pada lag yang benar dari 200 sampel. Korelasi dapat menjadi alat yang berguna untuk menentukan waktu tunda, ketika diterapkan pada kumpulan data yang berisi jenis kesamaan yang tepat.

Berdasarkan saran Muhammad, saya mencoba membuat naskah Matlab. Namun, saya tidak dapat menyimpulkan jika dia membangun distribusi Gaussian berdasarkan varians dan kemudian mengambil estimasi KDE atau dia melakukan estimasi KDE dengan asumsi Gaussian.

Juga sulit untuk menyimpulkan bagaimana ia menerjemahkan waktu offset KDE ke domain waktu. Ini usaha saya untuk itu. Setiap pengguna yang tertarik menggunakan skrip bebas untuk dan memperbarui versi yang disempurnakan jika memungkinkan.

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