Discrete Fourier Transform to find phase shift - Mathematica

[NeoDevin]

If I use the following code in Mathematica

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]

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

 

[NeoDevin]

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]

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)

 

[NeoDevin]

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]

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 I am components at the dominant frequency, and compare those to the phase shift in your original data. Do the Re and I am 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 I am (or whichever direction it rotates)?

 

[rbj]

i wouldn't even bother with the Fourier Transform. but if neither f₁[.] or f₂[.] 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 d₁ to 0?) then you can also generate (from the same parameters)

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

now, if you cannot generate f₁(.) (and thus not g₁(.)), and are given f₁ as data, then you need to Hilbert transform f₁ to get g₁. from f₁ and g₁, you cross-correlate these against f₂ and you will get relative amplitudes of the components of f₂ that is in-phase with f₁ and g₁. 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]

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.

 

[NeoDevin]

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]

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.