Both of them have the same dominant frequency in the FFT (and same sampling frequency), but the result I get is ϕ= 5.54164
Which is way too high, the relative phase should be much less than π/2 as can be seen by the curves.
Can anyone shed some light on what am I missing here?

EDIT2: Just to better understand what the phase of the FFTW means, I’m trying to reconstruct the signal using:

function sine_build(time, signal, fs)
F = rfft(signal)
idx = argmax(abs.(F))
freqs = rfftfreq(length(time), fs)
f = freqs[idx]
A = (maximum(signal) - minimum(signal)) / 2
ϕ = angle(F[idx])
sine = @. A * sin(2 * π * f * time - ϕ)
return sine
end

But for some reason, when I compare the result to the original series, the frequency and amplitude are fine, but the phase is off. I’m missing some crucial understanding about the phase of the FFT.

My first guess would be that y2 signal is not completely symmetric relative to zero and, as a result, has some DC component.
If you provide the signals data, it will be easier to investigate.

The figure actually looks fine to me, a phase difference of 5.54164 is equivalent to a phase difference of 2pi - 5.54164 i.e. 43 degrees, depending on which signal you take as reference. This is about what we see here:

using GLMakie
lines(0..30, sin)
lines!(0..30, x->sin(x+5.54))

I thought about it, but then why would I need to subtract the value from 2π ? Is there a way to extract the right phase difference universally? (I want to calculate the phase for a bunch of different pairs).

Note that 5.54 radians is just 0.882(2\pi) radians or +318 degrees, i.e. the FFT predicts that one signal is advanced by that phase angle relative to the other. Since the Fourier transform assumes pure (CW) sinusoids, so this is equivalent to saying that there’s approximately a 42 degree magnitude phase difference between the signals (360-318 or |318-360|).

Cool thank you, I didn’t think about it. That’s helpful.
Another interesting thing is that I find it difficult to reconstruct the signal using the phase from the FFT, for some reason the frequency is fine but the phase is off.

You determine the frequencies of both sine waves separately. But if that frequency is different between them, then the phase difference between them will rotate with the difference frequency, i.e. will not be fixed. If you know the frequency f of both waves, you could instead of using an FFT simply multiply both signals with phasor = cispi.(2*f*eachindex(sig1)/fs) of that same frequency f, and low-pass filter. The result will be a complex number the angle of which represents the phase of each sine wave. The angle of the quotient between these numbers for both signals will be the phase difference. Basically I/Q downconversion.

Also: if you use the FFT for frequency estimation, always multiply with a window function first.

Your calculation for magnitude A makes a lot of implicit assumptions, including that the signal is purely composed of a single frequency. This may or may not be a safe assumption, depending on the source of these signals. You could also estimate this magnitude directly from the magnitude of the FFT results.

The typical way to reconstruct a signal from FFT results is using the inverse FFT. If you only want the signal contribution at a single frequency, you could simply zero all of the other vector entries and then irfft the rfft result.

If you’re really set on implementing your own symbolic reconstruction, you might need to test whether you’re using the appropriate basis function. The FFT is fundamentally a function \mathbb{C}\mapsto\mathbb{C} on the exponential \exp(j\omega t)=\cos(\omega t) + j \sin (\omega t). When the time-domain signal is purely real-valued, you get equal contributions on two output frequencies \pm f_0 so that the imaginary-valued \sin contribution goes to zero. An rfft speeds up the math for real-valued inputs but produces the same result, so to reconstruct the real-valued signal directly you’d probably want to truncate the exponential above to just the real-valued component \cos(j\omega t).

using DSP, DataFrames
t1, w = 1:600, 2π/12
X = [sin.(w*t1) sin.(w*t1.+0.743)]
csp = mt_cross_power_spectra(permutedims(X))
pwr = power(csp)
D = DataFrame(freq=freq(csp), csp1=pwr[1,2,:], csp2=pwr[2,1,:])
D.power = abs.(D.csp1)
# Change to csp1 or csp2 as you prefer:
D.Δphase = angle.(D.csp2)
D[argmax(D.power),:]