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
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);
s2ana = hilbert(sig2_with_noise);
angle1 = unwrap(angle(s1ana));
angle2 = unwrap(angle(s2ana));
anglerel = unwrap(angle(s2ana./s1ana));