Find intersections of two sin wave function

8 views (last 30 days)
Zihao Hu on 25 Apr 2024
Commented: Mathieu NOE on 28 May 2024 at 17:42
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 on 25 Apr 2024
There is OptimalityTolerance and StepTolerance but not Tolerance

Mathieu NOE on 25 Apr 2024
hello
why not using this function available from Fex :
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 on 25 Apr 2024
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 on 28 May 2024 at 17:42
hello again
problem solved ?

David Goodmanson on 25 Apr 2024
Edited: David Goodmanson on 25 Apr 2024
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')
.

Categories

Find more on Fourier Analysis and Filtering in Help Center and File Exchange

R2024a

Community Treasure Hunt

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

Start Hunting!