Fast location of zero crossings with interpolation

27 次查看(过去 30 天)
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

采纳的回答

Star Strider
Star Strider 2014-2-25
编辑:Star Strider 2014-2-26
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.
  2 个评论
Eric Sampson
Eric Sampson 2014-3-3
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
Star Strider 2014-3-3
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.)

请先登录,再进行评论。

更多回答(1 个)

Shashank
Shashank 2014-3-6
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
  1 个评论
Star Strider
Star Strider 2014-3-6
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.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 MATLAB 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by