Fast location of zero crossings with interpolation

Hi all, I've got data that looks like this:

clean = sin(1:0.2:2e5);
noise = -0.03 + (0.06).*rand(size(clean));
data = clean + noise;

I need to locate the interpolated zero crossing values of this data (the predicted/interpolated X values, not just the nearest index). Right now I'm locating the nearest index to each zero crossing using sign change, which is easy/fast, and then calling interp1 with two points on each side of the nearest index (five points total). This works fine however it means that I'm calling interp1 about a million times, which adds up in terms of time :)

I'm having trouble coming up with a way to vectorize/speed this up. I've switched over to using interp1q which helps somewhat, but I still need to speed this operation up by quite a lot. Any ideas would be greatly appreciated!

Thanks, Eric

Answer 1

[Star Strider]

This tells me it takes about 2.7 seconds to calculate 63661 intercepts. (I didn’t run a comparison with interp1, so I can’t compare the times.) It seems to me to produce the result you want:

an = 1:0.2:2e5;
clean = sin(an);
noise = -0.03 + (0.06).*rand(size(clean));
data = clean + noise;

% Find zero crossings:
T0 = clock;
cnv = data .* circshift(data,[0 1]);
zx = find(cnv < 0);
% Interpolate zeros:
for k1 = 2:size(zx,2)-1
    ixrng = zx(k1)-2:zx(k1)+2;
    X = an(ixrng);
    Y = data(ixrng);
    b = [ones(size(X)); X]'\Y';
    Xi(k1) = -b(1)/b(2);
end
T1 = clock;
DT = etime(T1,T0)

figure(1)
plot(an, data)
hold on
% plot(an(zx), zeros(size(zx)), 'b*')
plot(Xi, zeros(size(Xi)), 'pr')
hold off
axis([0  25    ylim])
grid

It implements a simple linear least-squares regression, then solves for the X-intercepts.

[Eric Sampson]

Thanks! That is faster. I also found that I could use Doug Schwarz's INTERSECTIONS File Exchange function to do this same operation, with a horizontal line at Y=0 as the second 'curve', and it seemed even faster than the code above.

[Star Strider]

My pleasure!

I’ll take a look at Doug’s routine, since I like to see how people with more experience than I have with MATLAB solve specific problems. (I put my code together in about a half hour.)

Answer 2

[Shashank]

Hello Star Strider

Could you please explain bit about your code, specifically these two steps.

    b = [ones(size(X)); X]'\Y';
    Xi(k1) = -b(1)/b(2);

what are you trying to calculate ?

Thank you Shasha

[Star Strider]

Those two lines calculate the linear regression.

This line:

b = [ones(size(X)); X]'\Y';

calculates the parameters for the regression, Y = b(1) + b(2)*X. The ones vector creates the y-intercept, b(1). The transpose (') operator is necessary here because the data are in a row vector.

This line:

Xi(k1) = -b(1)/b(2);

calculates the X-intercept. Set Y = 0 and solve for X. I called that array Xi for ‘X-intercept’.

See Systems of Linear Equations in the MATLAB documentation for details.