Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Using Fourier Transform to get frequency and phase

Suppose I have samples of a periodic function, what is a good way to get frequency and phase information out of it?

In particular, I'd like to get a form like

a+b Cos[c x + d]

Here's a part of the sample

{255,255,255,249,64,0,0,0,0,0,0,0,0,0,0,0,0,233,255,255,255,255,255,255,255,255,255,209,0,0,0,0,0,0,0,0,0,0,0,0,118,255,255,255,255,255,255,255,255,255,255,132,0,0,0,0,0,0,0,0,0,0,0,0,200,255,255,255,255,255,255,255,255,255,239,19,0,0,0,0,0,0,0,0,0,0,0,46,245,255,255,255,255,255,255,255,255,255,186,0}
like image 723
Yaroslav Bulatov Avatar asked Feb 24 '23 15:02

Yaroslav Bulatov


1 Answers

Using Autocorrelation and FindFit[ ]

(*Your list*)
ListPlot@l

enter image description here

(*trim the list*)
l1 = Drop[l, (First@Position[l, 0])[[1]] - 1];
l2 = Drop[l1, Length@l1 - (Last@Position[l1, 0])[[1]] - 1];
(*autocorrelate*)
ListLinePlot@(ac = ListConvolve[l2, l2, {1, 1}])

enter image description here

(*Find Period by taking means of maxs and mins spacings*)
period = Mean@
   Join[
    Differences@(maxs = Table[If[ac[[i - 1]] < ac[[i]] > ac[[i + 1]], i, 
                               Sequence @@ {}], {i, 2, Length@ac - 1}]), 
    Differences@(mins = Table[If[ac[[i - 1]] > ac[[i]] < ac[[i + 1]], i, 
                               Sequence @@ {}], {i, 2, Length@ac - 1}])];

(*Show it*)
Show[ListLinePlot[(ac = ListConvolve[l2, l2, {1, 1}]), 
  Epilog -> 
   Inset[Framed[Style["Mean Period = " <> ToString@N@period, 20], 
     Background -> LightYellow]]], 
 Graphics[Join[{Arrowheads[{-.05, .05}]}, {Red}, 
   Sequence @@@ Arrow[{{{#[[1]], Min@ac}, {#[[2]], Min@ac}}}] & /@ 
    Partition[mins, 2, 1], {Blue}, 
   Sequence @@@ Arrow[{{{#[[1]], Max@ac}, {#[[2]], Max@ac}}}] & /@ 
    Partition[maxs, 2, 1]]]]

enter image description here

(*Now let's fit the Cos[ ] to find the phase*)
model = a + b Cos[x (2 Pi)/period + phase];
ff = FindFit[l, model, {a, b, phase}, x, 
             Method -> NMinimize, MaxIterations -> 100];

(*Show results*)
Show[ListPlot[l, PlotRange -> All, 
  Epilog -> 
   Inset[Framed[Style["Phase = " <> ToString@N@(phase /. ff), 20], 
     Background -> LightYellow]]], Plot[model /. ff, {x, 1, 100}]]

enter image description here

like image 121
Dr. belisarius Avatar answered Mar 19 '23 04:03

Dr. belisarius