# extracting phase from sinusoidal data

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!

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


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.

• 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!
– alex
Aug 26, 2021 at 21:40

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


• 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. :)
– alex
Aug 26, 2021 at 21:45
• 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. Aug 27, 2021 at 9:27
• 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. Aug 27, 2021 at 9:59

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


• 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.
– alex
Aug 28, 2021 at 11:18
• 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?