matlab - Time delay estimation of oscilloscope signals using cross correlation

I have recorded 2 signals from an oscope. They look like this:

I want to measure the time delay between them in Matlab. Each signal has 2000 samples with a sampling frequency of 2001000.5.

The data is in a csv file. This is what I have so far.
I erased the time data out of the csv file so that only the voltage levels are in the csv file.

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

This gives this result:

From what I've read I need to take the cross correlation of these signals and this should give me a peak relating to the time delay. However when I take the cross correlation of these signals I get a peak at 2000 which I know is not correct. What should I do to these signals before I cross correlate them? Just looking for some direction.

EDIT: after removing the DC offset this is the result I am now getting:

Is there a way to clean this up to get a more defined time delay?

EDIT 2: Here are the files:
http://dl.dropbox.com/u/10147354/scope1col.csv
http://dl.dropbox.com/u/10147354/scope2col.csv

Answer

@NickS

Since it is far from certain that the second signal in the plots is in fact a solely delayed version of the first, other methods besides the classical cross-correlation have to be attempted. This is because the cross-correlation (CC) is merely a maximum likelihood estimator if your signal(s) are delayed versions of each other. In this case, they are clearly not, to say nothing about the non-stationarity of them either.

In this case, I believe what may work is a time-estimation of the significant energy of the signals. Granted, 'significant' can or cant be somewhat subjective, but I believe that by looking at your signals from a statistical point of view, we will be able to quantify 'significant' and go from there.

To this end, I did the following:

STEP 1: Compute the signal envelopes:

This step is simple, as the absolute value of output of the Hilbert-Transform of each of your signals is computed. There are other methods to compute envelopes, but this is pretty straight forward. This method essentially computes the analytical form of your signal, in other words, the phasor representation. When you take the absolute value, you are destroying phase and only after energy.

Furthermore since we are pursuing a time-delay-estimate of the energy of your signals, this approach is warranted.

STEP 2: De-noise with edge-preserving non-linear Medial Filters:

This is an important step. The objective here is to smooth out your energy envelopes, but without destruction or smoothing out your edges and fast rise times. There is actually an entire field devoted to this, but for our purposes here, we can simply use an easy to implement non-linear Medial filter. (Median Filtering). This is a powerful technique because unlike mean filtering, medial filtering will not null out your edges, but at the same time 'smooth' out your signal without significant degradation of the important edges, since at no time is any arithmetic being performed on your signal (provided the window length is odd). For our case here, I selected a medial filter of window size 25 samples:

STEP 3: Remove Time: Construct Gaussian Kernel Density Estimation Functions:

What would occur if you looked at the above plot sideways instead of the normal way? Mathematically speaking, that means, what would you get if you projected every sample of our denoised signals onto the y-amplitude-axis? In doing this we will manage to remove time so to speak, and be able to study the signal statistics solely.

Intuitively what pops out of the figure above? While noise energy is low, it has the advantage that it is more 'popular'. In contrast, while the signal envelope that has energy is more energetic than the noise, it is fragmented across thresholds. What if we considered the 'popularity' as a measure of energy? This is what we will do with (my crude) implementation of a Kernel Density Function, (KDE), with a Gaussian Kernel.

To do this, every sample is taken and a gaussian function constructed using its value as a mean, and a pre-set bandwidth (variance) selected a-priori. Setting the variance of your gaussian is an important parameter, but you can set it based on noise statistics based on your application and typical signals. (I only have your 2 files to go off on). If we then construct the KDE Estimation, we get the following plot:

You can think of the KDE as a continuous form of a histogram so to speak, and the variance as your bin-width. However it has the advantage of guaranteeing a smooth PDF that we can then perform first and second derivitave calculus on. Now that we have the Gaussian KDEs, we can see where the noise samples peak in popularity. Remember that the x-axis here represents the projections of our data onto the amplitude space. Thus, we can see which thresholds the noise is the most 'energetic' in, and those tell us which thresholds to avoid.

In the second plot, the first-derivative of the Gaussian KDEs is taken, and we pick the abscissa of the first sample after first-derivative after the peak of the mixture of Gaussians to attain a certain value close to zero. (Or first zero-crossing). We can use this method and be 'safe' because our KDE was constructed of smooth Gaussians of reasonable bandwidth, and the first derivative of this smooth and noise-less function was taken. (Typically first-derivatives can be problematic in anything but high SNR signals since they magnify noise).

The black lines show then at what thresholds we would be wise to 'segment' the image at, such that we avoid the entire noise floor. If we then apply to our original signals, we attain the following plots, with the black lines indicating the start of energy of our signals:

This thus yields a $\delta{t} = 241$ samples.

I hope this helped.