7
$\begingroup$

I am dealing with data which are essentially are sinusoidal with various degrees of noise. I have an interest in extracting the phase of the data and rather than simply fitting a sine wave, I would like to extract it from the fourier space as there is evidence that there may be some level of phase drift in our data and I want to avoid fitting a chirp.

I am simulating the system with a pure tone version. Although I can extract the frequency just fine, I struggle to get the value for the phase when dealing with the ArcTan[Re,Im]. Am I missing something? Many thanks!

A = 0.5;             (* amplitude *)
fc = 10;              (*frequency *)
\[Phi] = 30  \[Pi]/180;    (* phase in radians *)
fs = 32 fc;       (* sample rate *)
t = Table[ A Sin[2 \[Pi] fc  x + \[Phi]], {x, 0, 2, 1/fs}];  (* the data *)

(* extract relevant parts from the fourier: Re, Im, Magnitude and Arctan[Im/Re]*)
{re, im, abs, arg} = #@Fourier[t, FourierParameters -> {-1, 1}] & /@ {Re, Im, Abs, Arg};  

(* generate the frequency space based on the sample rate used *)
freqSpace = Table[(n - 1.) fs/Length[abs], {n, Length[abs]}]; 
   
ListPlot[{freqSpace, abs}\[Transpose][[;; 50]], FrameLabel -> {"frequency Hz", "|X(f)|"}, Filling -> Axis]

enter image description here

unitiseArg = arg* Unitize@Threshold[abs,0.2];  (* extract the arctan based on located freq peak *)
ListPlot[{freqSpace, 180/\[Pi]  unitiseArg}\[Transpose][[;; 50]], GridLines -> {{10}, {\[Phi] 180 /\[Pi], -\[Phi] 180 /\[Pi]}}, FrameLabel -> {"Frequency", "Phase , deg"}]

enter image description here

Wouldn't one expect to read the phase here to be 30 deg rather than 55 deg?

Many thanks!

$\endgroup$
9
$\begingroup$

This is a good question and I may add this concern to my notes on Fourier here

Two points that may be helpful.

  1. Fourier takes cosine waves as the reference for phase.
  2. You need to adjust your frequency spacing to have one less point. The frequency goes from zero to one point less than the sample rate.

If we repeat your calculations using a cosine wave we get

A = 0.5;(*amplitude*)
fc = 10;(*frequency*)
ϕ = 30 π/180;(*phase in radians*)
fs = 32 fc;(*sample rate*)
t = Table[
  A Cos[2 π fc x + ϕ], {x, 0, 2 - 1/fs, 1/fs}];(*the data*)
(*extract relevant parts from the fourier:Re,Im,Magnitude and \
Arctan[Im/Re]*){re, im, abs, 
  arg} = #@Fourier[t, FourierParameters -> {-1, 1}] & /@ {Re, Im, Abs,
    Arg};

freq = Table[(n - 1.) fs/Length[abs], {n, Length[abs]}];

Plotting the magnitude gives

ListLinePlot[Transpose[{freq, abs}][[1 ;; 50]], 
 PlotRange -> {All, All}]

enter image description here

Also, the phase in degrees gives

ListLinePlot[Transpose[{freq, 360/(2 \[Pi]) arg}][[1 ;; 50]], 
 PlotRange -> {All, All}]

enter image description here

There are spikes in the phase due to numerical noise but at 10 Hz the phase is -30 deg as you expected. This may be seen from

TableForm[Transpose[{freq, 360/(2 \[Pi]) arg}][[19 ;; 25]], 
 TableHeadings -> {None, {"Freq./Hz", "Phase/deg"}}]

enter image description here

The general problem of finding the amplitude and phase from a sine wave is not trivial. See here for various techniques.

Further, fitting a sine wave to data is also not straightforward see here for details.

Hope that helps.

$\endgroup$
1
  • 1
    $\begingroup$ That is great! Thank you so much! I did stumble on your helper but didn't find what I wanted. What a tiny change and yet what a massive difference. I wonder now how would one go about taking the 'cleaning' the fourier of a data that has a high level of noise (but whose frequency is distinguishable). How do you ensure even number of cycles there? Many thanks as well for you links! $\endgroup$
    – alex
    Aug 26 at 21:40
5
$\begingroup$

Your data should cover an even multiple of periods. However, consider

t =Table[A Sin[2 \[Pi] fc x + \[Phi]], {x, 0, 2 , 1/fs}];

The last point for x==2 starts a new period. This leads to some leaking at the bottom of the peak as following plot shows:

ListPlot[{freqSpace, abs}\[Transpose][[;; 50]], 
 FrameLabel -> {"frequency Hz", "|X(f)|"}, Filling -> Axis, 
 PlotRange -> All]

enter image description here

However, if your data cover a multiple of periods:

t = Table[A Sin[2 \[Pi] fc x + \[Phi]], {x, 0, 2 - 1/fs, 1/fs}];

we get the following clean amplitude versus frequency plot that shows no leaking:

and the phase versus frequency plot:

enter image description here

This shows a phase of 60 degrees. But why is it not 30 degrees? The answer is, that you are looking at the argument of the Fourier coefficients and not the phase shift of the original signal.

If we call the first Fourier coefficient different from zero cp and the second coef. different from zero cm and if we remember (see the help of InverseFourier) that the base function for the pos. freq. is Exp[-I x] and for the neg. freq. Exp[I x] (note, the sign and note that the analyze function is Exp[I x] and the synthesize function is Exp[-I x]), we can plot the arguments of the pos. and neg. parts together with the original function:

{cp, cm} = 
 Select[Fourier[t, FourierParameters -> {-1, 1}] // Chop, # != 0 &]

Plot[{cp Exp[- I x] + cm Exp[I x], Arg[cp Exp[- I x]], 
  Arg[cm Exp[ I x]]}, {x, -Pi, 2 Pi}, PlotLabels -> "Expressions"]

enter image description here

As can be seen (look at the zero crossings), the original function has a phase shift of Pi/6, however, the Fourier coefficients have a shift of -Pi/3

The question remains, why is the phase shift of the Fourier coefficients not equal to the phase shift of the original function. The answer is that only the real part of cp Exp[-I x] and cm Exp[I x] contribute to the original function and only these contributions have the same phase shift as the original function. This can be seen in the following plot:

Plot[{cp Exp[-I x] + cm Exp[I x], Re[cp Exp[-I x]], 
  Re[cm Exp[I x]]}, {x, -Pi, 2 Pi}, PlotLabels -> "Expressions"]

enter image description here

$\endgroup$
3
  • 1
    $\begingroup$ Hi Daniel, thank you for the thorough explanation! May I ask how you get the values for cp and cm in mathematica? Also from @Hugh 's answer below, he mentioned that mathematica uses Cosine for the fourier, would that have resolved the discrepancy of \pi/2? -also there may be a typo in your text about the sign of Exp[+I x] for the negative coefficient. :) $\endgroup$
    – alex
    Aug 26 at 21:45
  • 2
    $\begingroup$ HI, I added the code for cp and cm. Then note, as you can read in the help of Fourier and InverseFourier, that the function used to get the Fourier coefficients is Exp[I x] and the function used for the back transformation is Exp[-I x], the conjugate complex. I also fixed the typos. $\endgroup$ Aug 27 at 9:27
  • 1
    $\begingroup$ Try also to set the phase shift to zero and look at the argument of the Fourier coefficients. You will see that the real part has no phase shit, but the imaginary part has a phase shift of Pi/2. And as inn the Fourier synthesis the imaginary get canceled, it is the phase shift of the real part that matters in the result. $\endgroup$ Aug 27 at 9:59
4
$\begingroup$

You and those that answered have shown that extracting the phase using Fourier functions is at best complicated. I understand you would prefer not to use NonlinearModelFit but why not just use it as it is so much more direct and gives you the ability to check the fit (both visually and measures of precision are available).

If there is phase drift, you'll see that in the fit and then you can attempt to model for that. So how could a method that ignores the error do better? It can't, of course.

As stated by @Hugh, there can be issues with fitting sine curves. But all you really need is a reasonable estimate of the frequency which you can get from a plot of the data. (There are also more automatic ways but that is for a different time.)

a = 0.5;(*amplitude*)
fc = 10;(*frequency*)
ϕ = 30 π/180;(*phase in radians*)
fs = 32 fc;(*sample rate*)
data = Table[{x, a Sin[2 π fc x + ϕ]}, {x, 0, 2, 1/fs}];(*the data*)

For a reasonable test we really need to add some noise to the party. I generally deal with natural populations and not seeing at least some noise is foreign to me.

SeedRandom[12345];
data[[All, 2]] = data[[All, 2]] + RandomVariate[NormalDistribution[0, 0.1], Length[data]];
ListPlot[data, AspectRatio -> 1/6, ImageSize -> Large]

Plot of data

nlm = NonlinearModelFit[data, {a0 Sin[2 π fc0 x + ϕ0], a0 > 0}, 
  {{a0, Mean[data[[All, 1]]]}, {fc0, 10.4}, {ϕ0, 0}}, x];
nlm["ParameterTable"]

Parameter table

phase = (ϕ0 /. nlm["BestFitParameters"])*180/π
(* 31.69364076930844` *)

Show[ListPlot[data, AspectRatio -> 1/6, ImageSize -> Large],
 Plot[nlm[x], {x, Min[data[[All, 1]]], Max[data[[All, 1]]]}, PlotStyle -> Red]]

Data and fit

$\endgroup$
6
  • $\begingroup$ Hi @JimB! NonlinearFit is my go-to when the dataset is generally speaking clean of noise. Unfortuntaely the data I'm dealing with are quite noise heavy (by nature of experiment not negligence! :P) and in those cases Mathematica cannot always make a good fit, especially when many of the variables for the fit can vary by more than 2 orders of magnitude and there are multiple tones. you can have a 2-step fitting process, sure, but it adds to computational time when you have 200GB+ of data to deal with. $\endgroup$
    – alex
    Aug 28 at 11:18
  • $\begingroup$ By all means, this is a fantastic answer and definitely completes this page! $\endgroup$
    – alex
    Aug 28 at 11:19
  • $\begingroup$ Good. But now I play Devil's Advocate: You want to model the signal but think that a technique that ignores the noise (variability in phase, frequency, response, etc.) is better than a technique (NonlinearModelFit) that at least deals with one aspect of the noise? Why not build in something that would account for the phase drift? Now NonlinearModelFit is unlikely to do that (it really only hands one type of additive error). In other words, without modeling the noise in addition to the signal, how do you know the estimates of the signal aren't influenced by the ignoring of the noise? $\endgroup$
    – JimB
    Aug 28 at 19:07
  • $\begingroup$ However I would argue that the main benefit of the fourier is that you can extract the raw values of the phase and you can then deal with them separately. By contrast, in the nonlinearfit you will fit and then you need to check whether the fit is a good one or not, which adds to an already computationally heavy task (for many and long files). As far as I understand there is no information that you can extract from fouier that you couldn't do from a nonlinearfit (assuming sinusoidal). it's just a speed thing for me. What do you think? $\endgroup$
    – alex
    Aug 29 at 22:25
  • 1
    $\begingroup$ "you need to check whether the fit is a good one or not, which adds to an already computationally heavy task." And a Fourier analysis does not need a check? That's why you asked the question in the first place. My religion (Statistics) says that all analyses need a check. That such a check might be onerous is not a reason to avoid checking (although many do ignore checks). What do I think? I think a Fourier analysis can be very useful. But having it figure out the error structure? Not so much. $\endgroup$
    – JimB
    Aug 30 at 5:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.