# Power spectrum density calulation with ARProcess

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:

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?

• Can you please advise me how this can be solved with mathematica. Probably a suitable function is ARProcess, but from the help description I am not able to find out how to continue. – mrz Jul 8 '16 at 10:18

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.

• How can the power spectrum be plotted with that data ? – mrz Jul 13 '16 at 10:01
• @mrz Using your S(f) or perhaps Mathematica's PowerSpectralDensity? – Taiki Jul 13 '16 at 10:32

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"]

ListLinePlot[Normal[tsm][[3]]/Abs[Table[