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

enter image description here

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?

  • $\begingroup$ 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. $\endgroup$
    – mrz
    Commented Jul 8, 2016 at 10:18

2 Answers 2

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


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.

  • $\begingroup$ How can the power spectrum be plotted with that data ? $\endgroup$
    – mrz
    Commented Jul 13, 2016 at 10:01
  • 1
    $\begingroup$ @mrz Using your S(f) or perhaps Mathematica's PowerSpectralDensity? $\endgroup$
    – Taiki
    Commented Jul 13, 2016 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"]
TraditionalForm[Normal[tsm][[3]]/Abs[1 - Total@
  MapIndexed[#1^(E^(-I*w*First@#2)) &, Normal[tsm][[2]]]]^2]

enter image description here

plot it

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

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.