How to mark curve intersections on a graph

42 次查看（过去 30 天）
Kristina
Kristina 2023-2-15
编辑： William Rose 2023-2-15
I've tried to look up the various answers in the community, but it doesn't seem to be getting me anywhere. I know there are file exchanges available such as intersections() and interx() but I wanted to be able to (hopefully) write a couple of extra lines to mark these intersections of the curves. I tried
i=find(xr==xt)
x=t(i)
y=xt(i) --- or even xr(i)
...
plot(x,y,'go','MarkerSize',8)
as well and it did not plot all of the intersections between the curves. Please let me know if there is a way to plot the intersections on the graph.
Current code
t=0:0.01:2;
n=0:0.01:2;
xt=cos(2*pi*t)+2*cos(12*pi*t);
xr=cos(2*pi*n);
hold on
plot(t,xt,'r','LineStyle','--')
plot(n,xr,'b')
hold off
xlabel('t (sec)')
title('fs = 5kHz')
legend('Original','Reconstructed')

请先登录，再回答此问题。

采纳的回答

Emma Smith Zbarsky
Emma Smith Zbarsky 2023-2-15
Ran in:
The challenge here is that the points of intersection are not exactly at points where you are evaluating your functions.
A quick solution that does just use the known values of your data might look like this:
t=0:0.01:2;
n=0:0.01:2;
xt=cos(2*pi*t)+2*cos(12*pi*t);
xr=cos(2*pi*n);
hold on
plot(t,xt,'r','LineStyle','--')
plot(n,xr,'b')
hold off
xlabel('t (sec)')
title('fs = 5kHz')
% New Code
idxNeg = xr<xt; % Find all the points where xr is below xt
findChangePts = diff(idxNeg); % Check for values where this changes
myIdx = ~(findChangePts==0); % Identify the indices of the values where the sign changes
myIdx = logical([0 myIdx]) | xr==xt ; % add the first point back in as well as any actual exact zero crossings
hold on
plot(t(myIdx),xr(myIdx),"y*")
hold off
legend('Original','Reconstructed','Intersections')
If you care about more precision and you are working from data you would need to add some interpolation (or in this case just make your step size smaller when defining t).

更多回答（1 个）

William Rose
William Rose 2023-2-15
编辑：William Rose 2023-2-15
Ran in:
@Kristina,
[edit: Correcting my spelling mistakes. No changes to the code.]
The reason that this is not as easy a problem as you might think is that the vectors xr and xt are never equal to one another at the times sampled. Therefore if you look for the times when xr=xt, you won't find any. Of course we know they must be equal at some instant, since the plot shows that the signals cross.
Strategy 1: March along the curves, looking for times where the plots cross. When you find a crossing, estimate the intersection as midway between the points before and after the crossing. By midway, I mean midway in time, and in the middle of the xr() and xt() before the crossing and xr() and xt() after the crossing. Not very sophisticated, but let's check out what this crude strategy does:
t=0:0.01:2;
xt=cos(2*pi*t)+2*cos(12*pi*t);
xr=cos(2*pi*t);
numint=0;
for i=1:length(xt)-1
if (xt(i)>xr(i) && xt(i+1)<=xr(i+1)) || (xt(i)<xr(i) && xt(i+1)>=xr(i+1))
numint=numint+1;
tint(numint)=(t(i)+t(i+1))/2; %estimated intersection time
xint(numint)=(xt(i)+xr(i)+xt(i+1)+xr(i+1))/4; %estimated x-value at intersection
end
end
plot(t,xt,'-b',t,xr,'-r',tint,xint,'k*')
legend('xt','xr','Int')
That is better than nothing, but not great. How can we do better?
Strategy 2: As in strategy 1, work along the curves one point at a time. As in strategy 1, find successive times when the lines cross. Unlike strategy 1, approximate each cirve with a straight line during the crossing interval. Find the time and x-value where the two lines intersect.
t=0:0.01:2;
xt=cos(2*pi*t)+2*cos(12*pi*t);
xr=cos(2*pi*t);
numint=0;
for i=1:length(xt)-1
if (xt(i)>xr(i) && xt(i+1)<=xr(i+1)) || (xt(i)<xr(i) && xt(i+1)>=xr(i+1))
numint=numint+1;
b1=xt(i); m1=xt(i+1)-xt(i); %slope, intercept for xt line
b2=xr(i); m2=xr(i+1)-xr(i); %slope, intercept for xr line
tint(numint)=t(i)+(t(i+1)-t(i))*(b2-b1)/(m1-m2); %estimated intersection time
xint(numint)=b1+m1*(b2-b1)/(m1-m2); %estimated x-value at intersection
end
end
plot(t,xt,'-b',t,xr,'-r',tint,xint,'k*')
legend('xt','xr','Int')
That's not so bad!
Good luck.

类别

在 Help Center 和 File Exchange 中查找有关 MATLAB 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by