Find intersections of two sin wave function
4 views (last 30 days)
Show older comments
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 = 0.079
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
T=25;
Fs = 44100;
t = 0:1/44100:T-1/44100
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
disp("Intersection points (refined):")
disp(refined_intersections)
1 Comment
Answers (2)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1678896/image.png)
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 ' );
2 Comments
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
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
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')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1679016/image.png)
.
0 Comments
See Also
Categories
Find more on Fourier Analysis and Filtering in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!