Stan Math Library  2.10.0
reverse mode automatic differentiation
autocorrelation.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUN_AUTOCORRELATION_HPP
2 #define STAN_MATH_PRIM_MAT_FUN_AUTOCORRELATION_HPP
3 
5 #include <unsupported/Eigen/FFT>
6 #include <complex>
7 #include <vector>
8 
9 
10 namespace stan {
11 
12  namespace math {
13 
14  namespace {
19  size_t fft_next_good_size(size_t N) {
20  if (N <= 2) return 2;
21  while (true) {
22  size_t m = N;
23  while ((m % 2) == 0) m /= 2;
24  while ((m % 3) == 0) m /= 3;
25  while ((m % 5) == 0) m /= 5;
26  if (m <= 1)
27  return N;
28  N++;
29  }
30  }
31  }
32 
53  template <typename T>
54  void autocorrelation(const std::vector<T>& y,
55  std::vector<T>& ac,
56  Eigen::FFT<T>& fft) {
57  using std::vector;
58  using std::complex;
59 
60  size_t N = y.size();
61  size_t M = fft_next_good_size(N);
62  size_t Mt2 = 2 * M;
63 
64 
65  vector<complex<T> > freqvec;
66 
67  // centered_signal = y-mean(y) followed by N zeroes
68  vector<T> centered_signal(y);
69  centered_signal.insert(centered_signal.end(), Mt2-N, 0.0);
70  T mean = stan::math::mean(y);
71  for (size_t i = 0; i < N; i++)
72  centered_signal[i] -= mean;
73 
74  fft.fwd(freqvec, centered_signal);
75  for (size_t i = 0; i < Mt2; ++i)
76  freqvec[i] = complex<T>(norm(freqvec[i]), 0.0);
77 
78  fft.inv(ac, freqvec);
79  ac.resize(N);
80 
81  /*
82  vector<T> mask_correction_factors;
83  vector<T> mask;
84  mask.insert(mask.end(), N, 1.0);
85  mask.insert(mask.end(), N, 0.0);
86 
87  freqvec.resize(0);
88  fft.fwd(freqvec, mask);
89  for (size_t i = 0; i < Nt2; ++i)
90  freqvec[i] = complex<T>(norm(freqvec[i]), 0.0);
91 
92  fft.inv(mask_correction_factors, freqvec);
93 
94  for (size_t i = 0; i < N; ++i) {
95  ac[i] /= mask_correction_factors[i];
96  }
97  */
98  for (size_t i = 0; i < N; ++i) {
99  ac[i] /= (N - i);
100  }
101  T var = ac[0];
102  for (size_t i = 0; i < N; ++i)
103  ac[i] /= var;
104  }
105 
122  template <typename T>
123  void autocorrelation(const std::vector<T>& y,
124  std::vector<T>& ac) {
125  Eigen::FFT<T> fft;
126  return autocorrelation(y, ac, fft);
127  }
128 
129 
130  }
131 }
132 
133 #endif
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:31
boost::math::tools::promote_args< T >::type mean(const std::vector< T > &v)
Returns the sample mean (i.e., average) of the coefficients in the specified standard vector...
Definition: mean.hpp:23
void autocorrelation(const std::vector< T > &y, std::vector< T > &ac, Eigen::FFT< T > &fft)
Write autocorrelation estimates for every lag for the specified input sequence into the specified res...

     [ Stan Home Page ] © 2011–2016, Stan Development Team.