Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Discrete Fourier Transform to find phase shift - Mathematica

  1. May 27, 2008 #1
    If I use the following code in Mathematica

    Code (Text):
    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?
     
  2. jcsd
  3. May 28, 2008 #2

    berkeman

    User Avatar

    Staff: 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...
     
  4. May 28, 2008 #3
    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.
     
  5. May 28, 2008 #4

    berkeman

    User Avatar

    Staff: 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)
     
  6. May 28, 2008 #5
    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?
     
  7. May 28, 2008 #6

    berkeman

    User Avatar

    Staff: Mentor

    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)?
     
  8. May 28, 2008 #7

    rbj

    User Avatar

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

    [tex] f_1(t) = \cos(\omega t + d_1) [/tex]

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

    [tex] g_1(t) = \sin(\omega t + d_1) [/tex]

    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

    [tex] A = \sum_n f_1[n] f_2[n] w[n] [/tex]

    [tex] B = \sum_n g_1[n] f_2[n] w[n] [/tex]

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

    [tex] d_2-d_1 = \arg\{A+iB\} = \arctan\left( \frac{B}{A} \right) [/tex]

    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.
     
  9. May 28, 2008 #8

    Dale

    Staff: Mentor

    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.
     
  10. May 28, 2008 #9
    rbj,
    From my experimentation, your method gives [itex]d_1 - d_2[/itex], but other than that, it's perfect, thanks.

    If we generalize it a little, so that

    [tex]f_1(t) = a_1\cos(\omega t + d_1)[/tex]
    [tex]f_2(t) = a_2\cos(\omega t + d_2)[/tex]

    Is there a way to get the relative amplitude [itex]a_1/a_2[/itex] from this method?
     
  11. May 29, 2008 #10

    rbj

    User Avatar

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

    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: Discrete Fourier Transform to find phase shift - Mathematica
Loading...