Power spectrum density calulation with ARProcess

Question

I have a measurement of $x$ coordinates of an oscillating particle at presence of noise taken at constant time steps of 1/60 s.

The corresponding data set can be obtained here:

https://drive.google.com/open?id=0B9wKP6yNcpyfTjE3UzNUZWlyaFk

I would like to measure the power spectrum density (PSD) using the autoregressive (AR) estimation method.

In the paper http://www.measurement.sk/2011/Angrisani.pdf (equation 1) it is supposed that the $x$ coordinates are the output of a linear system:

$x(n)=-\sum\limits_{m=1}^{p} a_{p,m} x(n-m) + \epsilon(n)$

where $x(n)$ is the analyzed signal sample at the time interval $n$, $a_{p,1}$, $a_{p,2}$, ... , $a_{p,p}$ are the model coefficients, $\{\epsilon(n)\}$ is a white noise process with variance $\sigma^2_p$, and $p$ is the model order. The PSD of a modeled variation of the $x$ coordinates in this way is totally described by the model parameters and the variance of the white noise process (equation 2):

$S(f)=\frac{\sigma^2_p T_S}{\left|1+ \sum\limits_{m=1}^{p} a_{p,m}e^{-j2\pi m f T_S}\right|}$ with $|f|\leq f_N$

where $T_S=1/f_S$ is the sampling interval and $f_N=1/(2T_S$) is the Nyquist rate.

How can I find the model coefficients ($a_{p,1}$, $a_{p,2}$, ... , $a_{p,p}$) and the noise variance $\sigma$? with mathematica?

I think that ARProcess is the right function, but I do not know how to use it?

Answer 1

data = Import["data.dat", "List"]
tsm = TimeSeriesModelFit[data, "AR"]
Normal[tsm]

Result:

ARProcess[0.00435059, {1.41701, -0.0804567, -0.245416, -0.0753475, -0.0415635}, 5.10131*10^-7]

Then you can read off the coefficients and noise variance.

Answer 2

Using this code,you can plot Power spectrum density

data = Flatten@Import["C:\\Users\\xiaoz\\Downloads\\20160708_coordinates_x_mm.dat"];
tsm = TimeSeriesModelFit[data, "AR"]
TraditionalForm[Normal[tsm][[3]]/Abs[1 - Total@
  MapIndexed[#1^(E^(-I*w*First@#2)) &, Normal[tsm][[2]]]]^2]

plot it

ListLinePlot[Normal[tsm][[3]]/Abs[Table[
1 - Total@MapIndexed[#1^(E^(-I*w*First@#2)) &, Normal[tsm][[2]]], 
{w, 0, Pi, 0.01}]]^2, DataRange -> {0, 60/2}]