Friday, June 27, 2008

Interpolating for Zero-Crossing Detection in MATLAB

Today someone who read my post on sinc interpolation in MATLAB e-mailed me to ask about how to interpolate between data points in MATLAB to do accurate zero-crossing detection. Here, I summarize my response.

Detecting zero crossings of (nearly) square waves that are sampled is difficult to do accurately because it pushes the limits of the "band-limited" signal approximation we use when sampling. That being said, if you'd like to use sinc interpolation, you should decide *HOW* close you need to be to the zero crossing. For example, your original signal might have 50 picoseconds between samples, and so without interpolation, you can detect zero crossings accurate within 50 picoseconds (ps). Maybe you need a closer result. How close do you need to be? 10ps? 5ps? 1ps? Less than that?

Figure out how close you need to be and then generate the appropriate upsampled time vector:
uptime = 0:step_size:final_time;
where step_size is your tolerance (e.g., 5 ps) and final_time is the maximum value of your sampled time vector. It is a good idea to engineer your step_size and final_time so that length(uptime)/length(old_time) is an INTEGER. It's also a good idea to make sure that all of your old_time points (except maybe your last point) exist in your uptime vector (i.e., your uptime vector just has time ADDED between old_time points).

Then you can use any interpolation mechanism you'd like. If you want to use the sinc_interp function that I wrote, do something like:
new_values = sinc_interp( old_values, old_time, uptime );
If you use sinc_interp, make sure you use the VECTORIZED version, which is MUCH faster. Actually, MATLAB provides a MUCH faster function that does NEARLY the exact same thing — the interpft command. In that case, you do
new_values = interpft( old_values, length(uptime) )
Otherwise, you might have luck using the resample command in MATLAB (it does NOT do pure sinc interpolation; it's more sophisticated, and it may not match at every data point):
new_values = resample( old_values, P, Q );
where P/Q is the INTEGER representing length(new_time)/length(old_time) (again, engineer your step_size to make this possible). ALTERNATIVELY, interp works OK too:
new_values = interp( old_values, P/Q );
Try all three and see which one works best for you.

Once you have this new interpolated vector, use some mechanism to detect zero crossings. Off the top of my head, here's the fastest way I can think of to do the zero-detection. Assume that we have
  • y (column vector of data -- e.g., new_values)
  • t (column vector of time -- e.g., uptime)
Also assume that t is already sorted. Then... (note that union sorts too)
i = union( find(y==0), find(conv(sign(y),[1 1]) == 0) );
should give you the MATLAB indexes where y crosses zero.

To test, try...
plot( t, y );
hold on;
stem( t(i), max(y)*ones(size(i)) );
hold off;
The result should be a plot of your data with spikes added at each zero crossing. That example also demonstrates how you recover the particular zero-crossing times with
t(i)
For example, when I apply this algorithm to a sin(t) waveform, t(i) gives me 0, pi, 2*pi, and 3*pi out, as expected.

No comments: