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

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.

• This is an interesting problem. How many samples ahead do you wish to predict? The function seems to be doing something quite complicated. The example expression has 8 unknows if you tried to fit the available data. I assume you don't know the functional form of the data -is that correct?
– Hugh
Commented Feb 15, 2022 at 16:09
• I am interested in predicting the next six samples. We also can't assume all the sines have amplitude of 1 as in my example, and we don't know how many sine terms are in the sum. We only know all the frequencies are 0.4 radians/sample or less. Commented Feb 15, 2022 at 17:10
• I don't think you fully answered @Hugh 's question: What exactly do you know about the process that generates the data? And what about that process that remains the same for future predictions? You give a single dataset and ask for good estimation method but methods are judged on multiple datasets rather than on single datasets. If some sort of a model for generating data isn't specified, then maybe using the last known observation as the prediction is the best one can do. Clearly you have in mind what the data generation process is. I just don't see that articulated in the question.
– JimB
Commented Feb 16, 2022 at 3:16
• We know the data is generated by a sum or sines and possibly a constant. The frequency of each sine does not exceed 0.4 radians/sample. We don't know the amplitude, frequency, or phase of the sines. Commented Feb 17, 2022 at 0:19

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.

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.

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]


• I learned a lot from your code, but sometimes it gives a poor prediction. I demonstrate that in an "answer" I posted. Commented Feb 16, 2022 at 0:26

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.