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.

sinc interpolation in MATLAB

UPDATE: The examples given here are meant to give mathematical insight into how sinc interpolation works by using a finite-time APPROXIMATION. Sinc interpolation necessarily involves a circular convolution, which is not a finite computation in the time domain. If you actually need to do sinc interpolation, use the interpft function. It does an FFT, pads the FFT with zeros, and does an IFFT. Consequently, it is VERY FAST. Moreover, using an FFT (or DFT, in general) is the only way to use a finite computation to do a sinc interpolation (or circular convolution in general). Just make sure that the resampled time vector you use has a length that is an integer multiple of your original time vector. Also make sure that it lands on the same points as your original time vector (i.e., it should only add new points between old points).

I've also created a interpftw that does the same job as interpft but allows you to reconstruct samples from different aliasing windows. Whereas interpft pads the end of the FFT with zeros, interpftw pads both the beginning and the end. In other words, the former is a low-pass filter and the latter is a bandpass filter.

If you search Google for sinc interpolation in MATLAB, many pages will reference the sinc_interp example from John Loomis. Unfortunately, I've found few sites that recognize that the function is not meant to do general purpose sinc interpolation. That is, it makes a few assumptions about the sampling rates that may not be evident to the average user.

So, I'm giving some of my students this example:

% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u)
    % Interpolates x sampled sampled at "s" instants
    % Output y is sampled at "u" instants ("u" for "upsampled")

    % Find the sampling period of the undersampled signal
    T = s(2)-s(1);

    for i=1:length(u)
        y( i ) = sum( x .* sinc( (1/T)*(u(i) - s) ) );
    end

    % Make sure y is same shape as u (row->row, col->col)
    y = reshape(y, size(u));
end

Here's a vectorized (i.e., MUCH FASTER) version:

% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u)
    % Interpolates x sampled at "s" uniformly spaced instants
    % Output y is sampled at "u" uniformly spaced instants
    % ("s" for "sampled" and "u" for "upsampled")
    % (consequently, length(x)=length(s))

    % Find the period of the undersampled signal
    T = s(2)-s(1);

    % The entries of this matrix are each u-s permutation.
    % It will be used to generate the sinc transform that will
    % be convolved below with the input signal to do the
    % interpolation.
    %
    % (recall that u(:) will be a column vector regardless
    % of the row-ness of u. So u(:) is a row, and s(:) is a
    % column)
    sincM = repmat( u(:), 1, length(s) ) ...
           - repmat( s(:)', length(u), 1 );

    % * Sinc is the inverse Fourier transform of the boxcar in
    % the frequency domain that was used to filter out the
    % ambiguous copies of the signal generated from sampling.
    % * That sinc, which is now sampled at length(u) instants,
    % is convolved with the input signal becuse the boxcar was
    % multipled with its Fourier transform.
    % So this multiplication (which is a matrix transformation
    % of the input vector x) is an implementation of a
    % convolution.
    % (reshape is used to ensure y has same shape as upsampled u)
    y = reshape( sinc( sincM/T )*x(:) , size(u) );
end

My function sinc_interp resamples the data in the x vector. The original time vector is given by the s vector and the new time vector is given by the u vector.

Alternatively, you can try this more advanced one that lets you pick which aliasing window you want to use to reconstruct (i.e., it's an ideal bandpass filter rather than just a low-pass filter).

% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u,N)
    % Interpolates x sampled sampled at "s" instants
    % Output y is sampled at "u" instants ("u" for "upsampled")
    % Optionally, uses the Nth sampling window where N=0 is DC
    %     (so non-baseband signals have N = 1,2,3,...)

    if nargin < 4
        N = 0;
    end

    % Find the period of the undersampled signal
    T = s(2)-s(1);

    for i=1:length(u)
        y( i ) = ...
            sum( x .* ...
                ( (N+1)*sinc( ((N+1)/T)*(u(i) - s) ) - ...
                  N*sinc( (N/T)*(u(i) - s) ) ) );
    end
end

Here's a vectorized (i.e., MUCH FASTER) version:

% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u,N)
    % Interpolates x sampled sampled at "s" instants
    % Output y is sampled at "u" instants ("u" for "upsampled")
    % Optionally, uses the Nth sampling window where N=0 is DC
    % (so non-baseband signals have N = 1,2,3,...)

    if nargin < 4
        N = 0;
    end

    % Find the period of the undersampled signal
    T = s(2)-s(1);

    % When generating this matrix, remember that "s" and "u" are
    % passed as ROW vectors and "y" is expected to also be a ROW
    % vector. If everything were column vectors, we'd do.
    %
    % sincM = repmat( u, 1, length(s) ) - repmat( s', length(u), 1 );
    %
    % So that the matrix would be longer than it is wide.
    % Here, we generate the transpose of that matrix.
    sincM = repmat( u, length(s), 1 ) - repmat( s', 1, length(u) );

    % Equivalent to column vector math:
    % y = sinc( sincM'(N+1)/T )*x';
    y = x*( (N+1)*sinc( sincM*(N+1)/T ) - N*sinc( sincM*N/T ) );
end