How to predict future samples from discrete samples of a low frequency signal?

Question

Suppose I have discrete samples of a low frequency signal as in: For example:

data=Table[Sin[0.12 n-0.3]+Sin[0.23 n+0.1]+Sin[0.34 n-0.5]+Sin[0.4 n+0.21],{n,0,52,1}];
ListLinePlot[data,Mesh->All]

What is a good way of predicting the next few samples when the frequencies and phases of the sine terms are unkown but we know the highest frequency does not exceed 0.4 radians/sample? It seems Predict, TimeSeriesForcast or AdjustTimeSeriesForecast would be good for that, but I am not familiar with Time Series Processes.

**** Update **** I made the problem more general by indicating that the sines may have different aplitudes.

Answer 1

If you know the form of the model, why not just use regular nonlinear regression?

data = Table[{n, Sin[0.12 n - 0.3] + Sin[0.23 n + 0.1] + Sin[0.34 n - 0.5] + 
     Sin[0.4 n + 0.21]}, {n, 1, 100, 1}];
nlm = NonlinearModelFit[data[[1 ;; 50]], Sum[a[i] Sin[w[i] n + θ[i]], {i, 4}],
  Flatten[Table[{a[i], w[i], θ[i]}, {i, 4}], 1], n, MaxIterations -> 10000, Method -> "NMinimize"];
ListPlot[{data, Table[{n, nlm[n]}, {n, 50, 100}]}, Mesh -> All, Joined -> {True, False}, 
 PlotStyle -> {{Blue, PointSize[0.015]}, {Red, PointSize[0.03]}}]

If you don't know the form of the model or more specifics as to how the data is generated, then no algorithm will be adequate or robust. (Also, you haven't added any noise to the data. Why?)

If you do know the form of the model as described below

$$y(n)=\sum_{i=1}^{k} a_i \sin(\omega_i n+\theta_i)+\epsilon_n$$

but you just don't know $k$ (the number of sine terms), then fit several plausible values of $k$ and use AIC (or AICc) to choose the best model. Then make the predictions.

Answer 2

Here is one approach. I take the Fourier transform of your data and then generate time histories with the same spectrum modulus but random phases. This means that the generated time histories have similar frequency content to your original data. Here is an example of that process.

nn = Length@data;
sr = 0.34;
ft = Fourier[data];
ff = Table[(n - 1) sr/nn, {n, nn}];
ListPlot[Transpose[{ff, Abs[ft]}], PlotRange -> {{0, sr/2}, All}]
SeedRandom[123];
a = RandomReal[{-π, π}, (nn - 1)/2];
ph = Join[{0}, a, -Reverse[a]];
th = InverseFourier[Table[Abs[ft[[n]]] E^(-I ph[[n]]), {n, 1, nn}]];
ListPlot[th, Mesh -> All]

The first plot is the spectrum and the second the regenerated data. We now have to find 6 points in the generated data that look like the end of the original data. I do this by using Nearest. There may be better methods. I put the nearest fit onto the original data as the extrapolation. I now run through this process a number of times to see what happens.

npts = 6;(* Number of points to extrapolate *)
nex = 4; (* Number of extrapolations to try *)
SeedRandom[123];
interps = Table[
   a = RandomReal[{-π, π}, (nn - 1)/2];
   ph = Join[{0}, a, -Reverse[a]];
   th = InverseFourier[Table[Abs[ft[[n]]] E^(-I ph[[n]]), {n, 1, nn}]];
   de = data[[-npts ;; -1]];
   {n} = Nearest[Partition[th, npts] -> "Index", de];
   ex1 = th[[npts (n - 1) + 1 ;; npts n]];
   diff = ex1[[1]] - de[[-1]];
   ex2 = ex1 - diff;
   ex = Join[data, ex2];
   ListPlot[ex, PlotStyle -> Red, ImageSize -> 6 72],
   {nex}];

Plot the extrapolated data and the original data



Column[Table[
      Show[interps[[n]], ListLinePlot[data]], {n, Length@interps}]]

If you follow the code I have moved the selected extrapolation so that it begins level with the last point of the original data.

I am sure this can be improved. Suggestions please.

Answer 3

dataAll = Table[Sin[0.12 n - 0.3] + Sin[0.23 n + 0.1] + Sin[0.34 n - 0.5] + 
    Sin[0.4 n + 0.21], {n, 0, 80, 1}];
data = dataAll[[;; 53]];

ts = TimeSeries[data];
tsm = TimeSeriesModelFit[ts, "ARIMA"];

ListLinePlot[{TimeSeries[dataAll], TimeSeriesForecast[tsm, {15}]}, 
 PlotRange -> {-3, 4}, Mesh -> All]

Answer 4

Doman gave an answer that gave very good results in one case. However, it does a poor job below.

dataAll=Table[Sin[0.12 n-0.3]+Sin[0.23n+0.1]+Sin[0.34 n-0.5]+Sin[0.4 n+0.21],{n,0,85}];
data=dataAll[[;;78]];
ts=TimeSeries[data];
tsm=TimeSeriesModelFit[ts,"ARIMA"];
ListLinePlot[{TimeSeries[dataAll],TimeSeriesForecast[tsm,{7}]},PlotRange->All,Mesh->All]

Perhaps there is no robust algorithm.