# Discrete Fourier Transform to find phase shift - Mathematica

• Mathematica
If I use the following code in Mathematica

Code:
f1[t_] := Cos[w t + d1]; f2[t_] := Cos[w t + d2];
data1 = Table[f1[t], {t,1,10000}]; data2 = Table[f2[t], {t,1,10000}];
ft1 = Fourier[data1]; ft2 = Fourier[data2];

To take the fourier transform of two data sets, how can I use the resulting data to determine the difference d2-d1?

berkeman
Mentor
Thread moved from Calculus to EE for a bit.

I'm not sure I understand what you are wanting to do. Are you trying to get some information about the phase shifts in the initial datasets after you take a Fourier transform of the datasets? Seems like a standard FT would not preserve that information...

That's exactly what I'm trying to do. I'm just using the test functions for now, once I figure it out, I will have two data sets, which (should) have the same dominant frequencies, but I need to find out the phase shift between them.

berkeman
Mentor
Does the FT in Mathematica give you complex data? If there is phase shift information, it seems like it would have to be in the complex domain of the FT output.... (I'm no expert in FT stuff, BTW)

Yes the output is complex, but it produces a range of frequencies, not just the dominant one. Should I just take the maximum (magnitude) of each transform, and divide them to find the phase shift?

berkeman
Mentor
Yes the output is complex, but it produces a range of frequencies, not just the dominant one. Should I just take the maximum (magnitude) of each transform, and divide them to find the phase shift?

Magnitudes don't help with phase information. Look at the Re and Im components at the dominant frequency, and compare those to the phase shift in your original data. Do the Re and Im components correspond to the original phase shift? Like, if you put in a 45 degree phase shift (PI/4) in the original sine wave data, do you get an FT where the complex components at that frequency show the PI/4 phase shift from Re to Im (or whichever direction it rotates)?

i wouldn't even bother with the Fourier Transform. but if neither f1[.] or f2[.] are a known sinusoid (both are measured), then i might bother with the Hilbert Transform. but let's say that the parameters of f1[.] are known. if you can generate (from the paramenters):

$$f_1(t) = \cos(\omega t + d_1)$$

(why not set d1 to 0?) then you can also generate (from the same parameters)

$$g_1(t) = \sin(\omega t + d_1)$$

now, if you cannot generate f1(.) (and thus not g1(.)), and are given f1 as data, then you need to Hilbert transform f1 to get g1. from f1 and g1, you cross-correlate these against f2 and you will get relative amplitudes of the components of f2 that is in-phase with f1 and g1. from those two relative amplitudes, apply the complex arg and you have your phase difference. explicitly

$$A = \sum_n f_1[n] f_2[n] w[n]$$

$$B = \sum_n g_1[n] f_2[n] w[n]$$

where w[n] is a suitable window function of your choice.

$$d_2-d_1 = \arg\{A+iB\} = \arctan\left( \frac{B}{A} \right)$$

the latter equality is true only if A>0. you need to use the atan2(B,A) function to put it in the correct quadrant.

Dale
Mentor
2021 Award
I think you will be better off using cross correlation than the FFT. For example:

ListCorrelate[data1,data2,1]

will give you a list with peaks at the locations corresponding to the shifts required to "match" the two signals.

rbj,
From my experimentation, your method gives $d_1 - d_2$, but other than that, it's perfect, thanks.

If we generalize it a little, so that

$$f_1(t) = a_1\cos(\omega t + d_1)$$
$$f_2(t) = a_2\cos(\omega t + d_2)$$

Is there a way to get the relative amplitude $a_1/a_2$ from this method?

rbj,
From my experimentation, your method gives $d_1 - d_2$, but other than that, it's perfect, thanks.

yeah, other than the fact that it was meant to work in opposite world where "up" is down, "left" is right, "positive" is negative, "right" equals wrong (but i thought it equalled left), "good" is evil. whatever. other than that, it's perfect.

If we generalize it a little, so that

$$f_1(t) = a_1\cos(\omega t + d_1)$$
$$f_2(t) = a_2\cos(\omega t + d_2)$$

Is there a way to get the relative amplitude $a_1/a_2$ from this method?

the same simple way you would measure magnitude gain of a filter. rectify, smooth, and divide. or you can square, smooth, divide, and square root. or you can Hilbert transform both. for each quadrature pair, square each and add, no need to smooth, divide, and square root. i dunno.