2
$\begingroup$

I am trying to fit the following tide data (taken from here) with a sinusoidal curve that looks like the following: displacement+amplitude * sin(t w - ph), where t is the independent variable, w is the angular frequency, ph is the phase, displacement is the average offset from zero, and amplitude is the sine curve's amplitude.

(* Times measured from midnight, 28 September (hours) *)
xd = {19.4000, 25.6333, 31.7833, 38., 44.1167, 50.3667, 56.4667, 
   62.6833, 68.8167, 75.0500, 81.1333, 87.3333, 93.5000, 99.7500, 
   105.8000, 112.0167};
(* Tidal heights offset from standard water level (m) *)
yd = {0.88, 12.71, 0.88, 12.63, 0.78, 12.59, 1.06, 12.57, 0.88, 12.32,
    1.40, 12.35, 1.15, 11.94, 1.84, 12.02};
samples = Transpose[{xd, yd}];

Thus far, I have tried four methods (LinearModelFit, Solve, NonlinearModelFit, and Fourier). This question concerns the last.

Here is the code I have used to produce the sine function using the discrete Fourier transform:

  SinusoidFitDFT[data_] := 
    Module[{xd, yd, deltaT, numberOfSamples, samplingPeriod, 
        samplingRate, displacement, amplitude, \[Omega], \[CurlyPhi]},
       xd = data[[1 ;; All, 1]];
 yd = data[[1 ;; All, 2]]; 
       deltaT = Mean[Differences[xd]]; numberOfSamples = Length[xd]; 
       samplingPeriod = numberOfSamples deltaT; samplingRate = 1/deltaT; 
       dft = Fourier[yd]; 
       maximumFrequencyPosition = (1. + Position[#1, Max[#1]][[1, 1]] &)[
         Abs[Rest[dft]]]; 
       frequencies = 
        Table[((i - 1) samplingRate)/
         numberOfSamples, {i, numberOfSamples}]; 
       maximumFrequencyContributor = 
        frequencies[[maximumFrequencyPosition]]; displacement = Mean[yd]; 
       amplitude = 
        Norm[Standardize[yd, Mean, 1 &], \[Infinity]]; \[Omega] = 
        2 \[Pi] maximumFrequencyContributor; \[CurlyPhi] = 
        Arg[dft[[maximumFrequencyPosition]]]; {displacement, 
        amplitude, \[Omega], \[CurlyPhi]}
       ];

This solution produces reasonable parameters, with the exception of the phase ([CurlyPhi]). Here is a plot of the original data (red dots) with the overlaid Sine curve (blue line). Included also are the values of the four fitted parameters and the root-mean square error.

{displacement, amplitude, \[Omega], \[CurlyPhi]} = 
 SinusoidFitDFT[samples]
fittedFunction = (displacement + 
     amplitude*Sin[\[Omega]*# - \[CurlyPhi]] &);
RootMeanSquare[yd - fittedFunction /@ xd]
Show[
 ListPlot[samples, PlotStyle -> Red],
 Plot[fittedFunction[x], {x, samples[[1, 1]], samples[[Length@xd, 1]]}]
 ]

ListPlot of Fourier curve fit

It is my understanding that the phase angle of the curve can be obtained by taking the ArcTangent of the result of the Fourier transform at the position of the maximum frequency contributor. Here are the intermediate results. First, the list of frequencies used to construct the transform:

Frequencies

Here is the Fourier transform of the tidal data:

Fourier Transform

Now the frequency spectrum plot:

Frequency Plot Spectrum

Everything looks okay to this stage. Note that the biggest frequency contributor (besides the zero-frequency component) occurs at the ninth position in the transform. The value of the frequency at this position is at about .08, which naturally corresponds to one high tide about every twelve hours.

Maximum Frequency Contributor

It is my understanding that the phase information for the curve can be obtained by looking at the (two-argument) arctangent of the value of the Fourier transform at this position (available via the Arg function). In the case of the maximum frequency, the phase is [Pi].

Phase

When I look at the original graph, however, it appears that the phase should be -[Pi]/2. That is, the first data point is at near the curve's minimum value. When that value is added into the fitted function, the error is reduced by a factor of twenty.

My questions are: What have I misunderstood about the way that Fourier works? How should I use the Fourier function to fit the data to a sine curve?

BTW, I am aware of this, this and this previous SE questions.

Edit 1:

I ran across an interesting mystery when I tried to refactor the solution provide below by Ulrich Neumann. First, I wrapped his code into the following function call:

MyFunc[] :=
  Module[{J, a, b, om, c, opt},
    J = Sum[(a + b Sin[om xy[[1]] + c] - xy[[2]])^2, {xy, samples}];
    opt = NMinimize[J, {a, b, c, om}];
    {a, b, c, om} /. opt[[2]]
  ];

Show[
  {displacement, amplitude, \[CurlyPhi], \[Omega]} = MyFunc[];
  Plot[
    displacement + amplitude Sin[\[Omega] t + \[CurlyPhi]],
    {t, Min[xd], Max[xd]},
    PlotRange -> All
  ],
  ListPlot[samples]
]

After doing so, the produced plot is identical to that shown in his answer. The anomaly comes when I try to do a simple renaming of a local variable in MyFunc:

MyFunc[] :=
  Module[{J, a, b, cx, c, opt},
    J = Sum[(a + b Sin[cx xy[[1]] + c] - xy[[2]])^2, {xy, samples}];
    opt = NMinimize[J, {a, b, c, cx}];
    {a, b, c, cx} /. opt[[2]]
  ];

Note that the variable om has become cx. When the code is now run, the results look much different:

enter image description here

This strongly violates my intuition that changing a local variable name should not affect the computed results. BTW, I did a Remove["*"];` before each execution to ensure a clean environment. Can anyone explain to me what is happening here?

Edit 2:

I tried a somewhat different approach to get around the problems mentioned by the respondants (unequal increments, incomplete cycles, and too few data points). I took the original data and constructed an interpolation function. I found the roots of this function and extracted values at equally spaced intervals, starting at the first root and continuing through a full number of cycles. I then supplied the extracted values to the Fourier function. Unfortunately the results were no better: the phase angle is still off. The bottom line is that I still don't know how to get a meaningful phase angle out of the Fourier function.

$\endgroup$
  • $\begingroup$ I think the problem is that your original data is not on a uniform grid, especially because it contains an offset 19.4000 (the first entry of xd) which should dictate the phase of the signal, but for your phase analysis you just do a discrete Fourier transform of your yd values, which means from the perspective of Fourier the entries switch from minimum value to maximum value on the grid at Nyquist frequency, which results in the phase angle of Pi. You need to incorporate your xd values to get a meaningful result. $\endgroup$ – Thies Heidecke Oct 30 '18 at 16:02
  • $\begingroup$ I agree with your observation. Is there a standard technique for doing the "incorporation"? $\endgroup$ – Spencer Rugaber Oct 31 '18 at 1:58
3
$\begingroup$

I'm unsure, wether the number of samples is sufficient to apply FourierSeries.

But you can identify a sine-function directly:

J = Sum[(a + b Sin[om xy[[1]] + c] - xy[[2]])^2, {xy, samples}];
opt = NMinimize[J, {a, b, c, om}]
(*{0.293904, {a -> 6.79905, b -> 5.9243, c -> 7.19992, om -> 0.515526}}*)    

The results fits your samples very good

Show[{Plot[a + b Sin[om t + c] /. opt[[2]], {t, Min[xd], Max[xd]},PlotRange -> All], ListPlot[samples]}]

enter image description here

$\endgroup$
  • $\begingroup$ Great solution! I think because of the uneven spacing in xd a global optimization like you did is a better approach than a Fourier transform in this case. If the number of points is much larger this might change, but for the problem as it is stated here it's perfect. $\endgroup$ – Thies Heidecke Oct 30 '18 at 21:34
  • $\begingroup$ @Ulrich Neumann: Is there a rule of thumb for the number of samples required to produce a quality solution? I am aware of other approaches. I have tried LinearModelFit, Solve, and NonlinearModelFit, but I was unaware of NMinimize. I used RMSError to compare the outcomes, and NonIinearModelFit was the best--somewhat better than NMinimize. I am still interested in getting Fourier to work by addressing the problem Thies Heidecke identified. Any suggestions would be appreciated. $\endgroup$ – Spencer Rugaber Oct 31 '18 at 2:22
  • $\begingroup$ @Spencer Rugaber The rule is the "Nyquist–Shannon sampling theorem" which states that you have to sample with at least the double of the frequence you want to detect. Just one point: The uneven spacing is unusual. But if you interprete the FourierSeries as a trigonometric approximation that makes no problem. $\endgroup$ – Ulrich Neumann Oct 31 '18 at 7:17
  • $\begingroup$ @Ulrich Neumann I noticed that your solution has changed the sign on the phase angle from what I had originally used. When I tried changing it back, I got a much different "solution", specifically with a frequency approximately three times as high. I am guessing that NMinimize is sensitive to the specific situation that it is given, but I am interested to learn how you knew to change the sign in order to get such a nice fit. Is there a rule of thumb that I should use in situations like this? $\endgroup$ – Spencer Rugaber Nov 5 '18 at 16:17
  • $\begingroup$ @Spencer Rugaber What means "tried to change the sign"? You should be aware, that the phase angle isn't unique. Perhaps you can change the optimization to opt = NMinimize[{J,0<=c<= 2 Pi}, {a, b, c, om}] $\endgroup$ – Ulrich Neumann Nov 6 '18 at 8:05
2
$\begingroup$

One way to approach this uses extensions of the Fourier transform, or of periodograms, to irregularly spaced data. I use some code from a talk that took place a few years ago:

https://materials.dagstuhl.de/files/15/15251/15251.DanielLichtblau.Slides.pdf

I will repeat the input data for purposes of making this self contained.

xd = {19.4000, 25.6333, 31.7833, 38., 44.1167, 50.3667, 56.4667, 
   62.6833, 68.8167, 75.0500, 81.1333, 87.3333, 93.5000, 99.7500, 
   105.8000, 112.0167};
yd = {0.88, 12.71, 0.88, 12.63, 0.78, 12.59, 1.06, 12.57, 0.88, 12.32,
    1.40, 12.35, 1.15, 11.94, 1.84, 12.02};
samples = Transpose[{xd, yd}];
mean = Mean[yd];

The irregularly spaced discrete FT can be defined as below.

idft[w_] := Total[(yd - mean) Exp[-I w xd]]/Length[samples]

We'll take a look at the picture.

Plot[Abs[idft[w]], {w, .1, 1}, PlotRange -> All]

enter image description here

We want to extract the peak value, the point at which it is hit, and also we need phase information.

{amp, f} = NMaximize[{Abs[idft[w]], 0 <= w <= 1}, w]
freq = (w /. f)

(* Out[138]= {5.64033025475, {w -> 0.509018083073}}

Out[139]= 0.509018083073 *)

There are probably better ways to get the phase, but I'll use the one from those slides. We first define and compute the Lomb-Scargle periodogram.

stot[w_] := Total[Sin[2*w*xd]]
ctot[w_] := Total[Cos[2*w*xd]]
tau[w_] := ArcTan[stot[w]/ctot[w]]
mean = Mean[yd];
var = Variance[yd];
periodogramLS[w_?NumberQ] := 
 Abs[Total[(yd - mean)*Exp[-I*(w*xd - tau[w])]]]^2/((2*var)*
  Total[Abs[Exp[-I*w*(xd - tau[w])]^2]])

Again we'll take a look at the plot.

Plot[periodogramLS[w], {w, .01, 1}, PlotRange -> All]

enter image description here

We now use the frequency to get the phase.

phase = tau[w /. f]

(* Out[153]= 0.952885595031 *)

Now plot this along with the data points. Here I borrow and slightly modify code by @UlrighNeumann, posted in a different response.

Show[{Plot[mean + amp*Sin[freq*t + phase], {t, Min[xd], Max[xd]}, 
   PlotRange -> All], ListPlot[samples, ColorFunction -> (Red &)]}]

enter image description here

Based on the plot, the values found all seem to be quite reasonable.

$\endgroup$
  • $\begingroup$ I think that "values" in the definition of periodogramLS should by yd. $\endgroup$ – Spencer Rugaber Nov 14 '18 at 3:39
  • $\begingroup$ Thanks @spencerrugaber, I knew I was going to forget some of the variable renamings in going from the talk slides to the current usage. I fixed it. $\endgroup$ – Daniel Lichtblau Nov 14 '18 at 14:50

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.