MATLAB Answers

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

1 view (last 30 days)
Jay Vaidya
Jay Vaidya on 17 Dec 2020
Edited: David Goodmanson on 19 Dec 2020
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:
  3 Comments
Matt Gaidica
Matt Gaidica on 17 Dec 2020
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?

Sign in to comment.

Answers (2)

VBBV
VBBV on 17 Dec 2020
%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.
  4 Comments
VBBV
VBBV on 17 Dec 2020
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

Sign in to comment.


David Goodmanson
David Goodmanson on 18 Dec 2020
Edited: David Goodmanson on 19 Dec 2020
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)

Tags

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!