extracting phase from sinusoidal data

Question

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]

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"}]

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

Many thanks!

Answer 1

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}]

Also, the phase in degrees gives

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

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"}}]

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.

Answer 2

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]

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:

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"]

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"]

Answer 3

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]

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

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