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.
