Find intersections of two sin wave function

Hi guys,

I am trying to find the intersection points of these two sin waves. I dont know what I did wrong.

f_1=120

f_1 = 120

f = 0.079

f = 0.0790

f1 = @(t) 35 * sin(2*pi*f_1*t + 0.08*pi).^2;
f2 = @(t) 35 * sin(2*pi*f*t + 0.08*pi).^2;
tolerance = 1e-5

tolerance = 1.0000e-05

T=25;
Fs = 44100;
t = 0:1/44100:T-1/44100

t = 1x1102500

1.0e+00 * 0 0.0000 0.0000 0.0001 0.0001 0.0001 0.0001 0.0002 0.0002 0.0002 0.0002 0.0002 0.0003 0.0003 0.0003 0.0003 0.0004 0.0004 0.0004 0.0004 0.0005 0.0005 0.0005 0.0005 0.0005 0.0006 0.0006 0.0006 0.0006 0.0007

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

plot(t, f1(t), t, f2(t));
xlabel('Time (s)');
ylabel('Amplitude');
title('Potential Intersection Regions');
legend('f1(t)', 'f2(t)');

difference = abs(f1(t) - f2(t));
potential_intersections = t(difference < tolerance);

refined_intersections = [];
for i = 1:length(potential_intersections)
    p_intersection = potential_intersections(i);
    try
        intersection_func = @(t) f1(t) - f2(t);
        options = optimoptions('fsolve'); 
        options.Tolerance = tolerance; 
        refined_intersection = fsolve(intersection_func, p_intersection, options);
        refined_intersections.append(refined_intersection)
    catch ME
        if strcmp(ME.identifier, 'MATLAB:fsolve:NotVectorValues')
            warning('fsolve did not converge for guess: %f', p_intersection);
        else
            rethrow(ME); 
        end
    end
end

Unrecognized property 'Tolerance' for class 'optim.options.Fsolve'.

disp("Intersection points (refined):")
disp(refined_intersections)

[Walter Roberson]

There is OptimalityTolerance and StepTolerance but not Tolerance

Answer 1

[Mathieu NOE]

hello

why not using this function available from Fex :

Fast and Robust Curve Intersections - File Exchange - MATLAB Central (mathworks.com)

I modified a bit the signal frequencies and duration to make the result easier to see here

if the function is getting too slow or creating out of memory issues on long signals, we could split it into smaller segments and concatenate the results

f_1 = 12;
f = 0.79;
f1 = @(t) 35 * sin(2*pi*f_1*t + 0.08*pi).^2;
f2 = @(t) 35 * sin(2*pi*f*t + 0.08*pi).^2;
T=1;
Fs = 44100;
t = 0:1/44100:T-1/44100;

[X0,Y0,I,J] = intersections(t, f1(t), t, f2(t),1);

plot(t, f1(t), t, f2(t),X0,Y0,'dk');
xlabel('Time (s)');
ylabel('Amplitude');
title('Potential Intersection Regions');
legend('f1(t)', 'f2(t)' , ' intersections ' );

[Mathieu NOE]

a home made solution (with linear interpolation) - much faster

beside that , I don't understand the need to sample at 44100 Hz for sine waves with max freq = 120 Hz

Fs around 2 kHz would largely suffice , especially here where we use interpolation to find the "true" intersection point

clc
clearvars

f_1 = 12;
f = 0.79;
f1 = @(t) 35 * sin(2*pi*f_1*t + 0.08*pi).^2;
f2 = @(t) 35 * sin(2*pi*f*t + 0.08*pi).^2;
T=1;
Fs = 44100;
t = 0:1/Fs:T-1/Fs;

% % first solution 
% tic 
% [X0,Y0,I,J] = intersections(t, f1(t), t, f2(t),1);
% toc
% 
% figure
% plot(t, f1(t), t, f2(t),X0,Y0,'dk');
% xlabel('Time (s)');
% ylabel('Amplitude');
% title('Potential Intersection Regions');
% legend('f1(t)', 'f2(t)' , ' intersections ' );

% second solution 
tic 
diff = f1(t) - f2(t);
[ZxP,ZxN] = find_zc(t,diff,0); % "zero crossing points with linear interpolation)
Xi = sort([ZxP;ZxN]);
Yi = interp1(t,f2(t),Xi);
toc

Elapsed time is 0.015296 seconds.

figure
plot(t, f1(t), t, f2(t),Xi,Yi,'dk');
xlabel('Time (s)');
ylabel('Amplitude');
title('Potential Intersection Regions');
legend('f1(t)', 'f2(t)' , ' intersections ' );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)

% put data in rows
x = x(:);
y = y(:);

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

% negative slope "zero" crossing detection, using linear interpolation
zci = @(data) find(diff(sign(data))<0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

[Mathieu NOE]

hello again

problem solved ?

Answer 2

[David Goodmanson]

Hello ZH

The symbol f is a bit overused in your example, but if you have two functions

g1 = A*sin(2*pi*f1*t +phi).^2 
g2 = A*sin(2*pi*f2*t +phi).^2 

then the points of intersection are the union of points (assume f1 < f2 )

tn = (n/2)/(f2-f1)    and        tm  = ((m/2)-phi/pi)/(f2+f1)

where both n and m are sets of consecutive integers

na<=n<=nb  and ma<=m<=mb

***** This is basically a beat frequency situation, hence the equal spacing of tn and of tm. *****

The values of n and m need to cover all the possible intersections of g1 and g2 for the given domain of t. One way to establish n and m is shown in the code below, although there may be an off-by-one problem at the ends of the time domain.

f1 = .877;
f2 = 10;
phi = pi/7;
t = 0:.001:1;
g1 = sin(2*pi*f1*t +phi).^2;
g2 = sin(2*pi*f2*t +phi).^2;
nb = round(max(t)*2*(f2-f1));
na = round(min(t)*2*(f2-f1));
tn = (na:nb)*(1/2)/(f2-f1);                 % crossing times
mb = round(2*(max(t)*(f2+f1)+phi/pi));
ma = round(2*(min(t)*(f2+f1)+phi/pi));
tm = ((ma:mb)*(1/2)-phi/pi)/(f2+f1);        % crossing times
gn = sin(2*pi*f2*tn +phi).^2;
gm = sin(2*pi*f2*tm +phi).^2;
fig(1)
plot(t,g2,t,g1,tn,gn,'ok',tm,gm,'ok')

.