How to find the delta cycles (change in the number of cycles) in two sine signals with nearly identical frequencies?

I have these two signals generated from the following code;

clc
clear all
close all

%% 
f = 4000;
t = 0:1e-5:50e-3;
sig1 = sin(2*pi*(f)*t);
sig2 = sin(2*pi*(f+50)*t);

noiseAmplitude = 0.05;
sig1_with_noise = sig1 + noiseAmplitude * randn(1, length(sig1));
sig2_with_noise = sig2 + noiseAmplitude * randn(1, length(sig2));

%%
figure
plot(t,sig1_with_noise)
hold on
plot(t,sig2_with_noise)
ylim([-1.1 1.1]);

These are two sinusoids (with a small noise) with nearly identical frequencies. Their frequencies differ by 50 Hz as you can see in the above code. But since these are not perfectly identical, one of the waveforms overtakes the other as a function of time, meaning the number of cycles over time from one signal will be greater than the other one. I want to find the delta_cycles Vs time plot keeping the sig_1 waveform as the reference. Currently, I can just tell by looking that the delta_cycles will be ~ 2.5 cycles over the time duration that showed in the top figure. But I want to plot the delta_cycles Vs time. Should I just take the difference? Or should I use the derivative to see the change? Thanks in advance.

Plot:

[Matt Gaidica]

It's unclear: will these be signals you don't have any information about?

[Jay Vaidya]

I will have only the time domain data of those signals. And about the information: they are not sinusoidal signals, but they can be any periodic signal.

[Matt Gaidica]

I feel like this might be a time when you would want to use a FFT to identify the frequency of each signal? It's hard to say without knowing you application. With very noisy data you might miss peaks and valleys of data, is that a fair assumption?

Answer 1

[VBBV]

%ue
plot(t,(sig1_with_noise-sig2_with_noise)./(2*pi))

As you said, you can take the difference between two signals and divide by 2*pi to get the delta number of cycles and plot it.

[Jay Vaidya]

Yes, I did that. But the difference in the number of cycles should continuously increase with time, monotonically. This plot gives just where the two signals cross each other. In other words, delta_cycles needs to keep increasing. (assuming that the sig_1 freq is larger than sig_2).

[VBBV]

The delta cycle needs to be constant ignoring noise and show as straightline. It cannot increase continuously with time as the difference between the two signals is always 50 at any instant

[Jay Vaidya]

No, delta_cycles is going to be an increasing with time because delta_cycles is not the difference in the time periods. What you are referring to is the difference in the time periods/frequency. But what I am intending to do is:

Let's say in a sample time of 50ms, because the frequency of the 2nd signal is slighlty higher than the first one, so within the same time (50ms) sig_2 will undergo more oscillations compared to the first one. Do you agree?

Which means, that the

delta_cycles  = (number_of_cycles_of_sig_1) - (number_of_cycles_of_sig_2)

is not constant but a increasing function as the

number_of_cycles_of_sig_1>number_of_cycles_of_sig_2

in 50ms.

Also,

number_of_cycles_of_sig_1>number_of_cycles_of_sig_2

because the

freq_sig_1 < freq_sig_2;

[VBBV]

I agree with you that 2nd signal will have more oscillations than first as it's frequency is higher than first BUT as both signals periodic in nature, at some some point in time they must cross each other. Do you agree? So that way the plot which you get should be correct since you can see the delta cycles increase and then decrease periodically rather than a linear increase in time

Answer 2

[David Goodmanson]

Hi Jay,

Using the hilbert transform on the signal gives the so-called analytic signal. The transform creates an imaginary part and adds it to the original signal (the real part) to gilve a complex signal of the form

const*exp(2*pi*i*f*t)

as in fig(1).

[NOTE that the Matlab definition of the hilbert transform is NONstandard. The actual hilbert transform gives just the red waveform in fig 1. Mathworks adds in the original signal (blue) as well, to create the full analytic signal ]

Finding the angle of the analytic signal, and using unwrap, shows an angle that increases linearly with with t as in fig(2). As you can see the higher freq signal has the steeper slope. Taking the ratio of the two analytic signals gives the relative phase as in fig(3). Since the difference frequency is 50 Hz and the total time is 50 msec, you would expect the accumulated phase difference to be

2*pi*50*50e-3 = 15.7 radians

which is what happens.

The code works with f = 4000, but I dropped f to 1000 to make things easier to see.

f = 1000;
t = 0:1e-5:50e-3;
sig1 = sin(2*pi*(f)*t);
sig2 = sin(2*pi*(f+50)*t);
noiseAmplitude = 0.0;
sig1_with_noise = sig1 + noiseAmplitude * randn(1, length(sig1));
sig2_with_noise = sig2 + noiseAmplitude * randn(1, length(sig2));
s1ana = hilbert(sig1_with_noise);   % analytic signal
s2ana = hilbert(sig2_with_noise);   % analytic signal   
figure(1)
plot(t,real(s1ana),t,imag(s1ana))
xlim([0 10e-3])
grid on
angle1 = unwrap(angle(s1ana));
angle2 = unwrap(angle(s2ana));
figure(2)
plot(t,angle1,t,angle2)
legend(['sigl'; 'sig2'],'location','east')
anglerel = unwrap(angle(s2ana./s1ana));
figure(3)
plot(t,anglerel)