Weird phase result between two simple sine curves

Specific Domains Signal and Image Processing
1

Hi all,
I’m trying to get the relative phase between two simple sine series:
sine

I’m trying to use FFTW to get their angle from the dominant frequency and calculate the relative phase:

function FFT_phase_rel(sig1, sig2, t, fs)
	F1 = rfft(sig1)
	F2 = rfft(sig2)
	idx = argmax(abs.(F1))
	idx2 = argmax(abs.(F2))
	ϕ1 = angle(F1[idx])
	ϕ2 = angle(F2[idx2])
	ϕ = abs(ϕ1 - ϕ2)
	return ϕ
end

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?

EDIT: the data is here (jld2 files): https://file.io/4UCDYj17qcdZ

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.

Thank you!

2

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.

3

Thank you, I added the data in an edit

4

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))

sin
sin1720×1339 243 KB

2 Likes
5

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).

6

Maybe you want something like this:

mod2pi(ϕ1 - ϕ2 + π) - π

to always have a value in [-\pi,\pi). Then take the absolute value if you want.

1 Like
7

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|).

8

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.

9

What code are you using to reconstruct the signal?

10

I added it to the original post above, thank you!

11

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.

12

Couple of notes:

  • 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).
13

There’s also the cross-spectrum:

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),:]
Row freq csp1 csp2 power Δphase
88 0.0849609 30.628-28.1375im 30.628+28.1375im 41.5908 0.743043
1 Like