Stan Math Library  2.20.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 
6 #include <unsupported/Eigen/FFT>
7 #include <complex>
8 #include <vector>
9 
10 namespace stan {
11 namespace math {
12 namespace internal {
17 inline size_t fft_next_good_size(size_t N) {
18  if (N <= 2)
19  return 2;
20  while (true) {
21  size_t m = N;
22  while ((m % 2) == 0)
23  m /= 2;
24  while ((m % 3) == 0)
25  m /= 3;
26  while ((m % 5) == 0)
27  m /= 5;
28  if (m <= 1)
29  return N;
30  N++;
31  }
32 }
33 } // namespace internal
34 
55 template <typename T>
56 void autocorrelation(const std::vector<T>& y, std::vector<T>& ac,
57  Eigen::FFT<T>& fft) {
58  using std::complex;
59  using std::vector;
60 
61  size_t N = y.size();
62  size_t M = internal::fft_next_good_size(N);
63  size_t Mt2 = 2 * M;
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  for (size_t i = 0; i < N; ++i) {
82  ac[i] /= (N - i);
83  }
84  T var = ac[0];
85  for (size_t i = 0; i < N; ++i)
86  ac[i] /= var;
87 }
88 
109 template <typename T, typename DerivedA, typename DerivedB>
110 void autocorrelation(const Eigen::MatrixBase<DerivedA>& y,
111  Eigen::MatrixBase<DerivedB>& ac, Eigen::FFT<T>& fft) {
112  size_t N = y.size();
113  size_t M = internal::fft_next_good_size(N);
114  size_t Mt2 = 2 * M;
115 
116  // centered_signal = y-mean(y) followed by N zeros
117  Eigen::Matrix<T, Eigen::Dynamic, 1> centered_signal(Mt2);
118  centered_signal.setZero();
119  centered_signal.head(N) = y.array() - y.mean();
120 
121  Eigen::Matrix<std::complex<T>, Eigen::Dynamic, 1> freqvec(Mt2);
122  fft.SetFlag(fft.HalfSpectrum);
123  fft.fwd(freqvec, centered_signal);
124  // cwiseAbs2 == norm
125  freqvec = freqvec.cwiseAbs2();
126 
127  Eigen::Matrix<T, Eigen::Dynamic, 1> ac_tmp(Mt2);
128  fft.inv(ac_tmp, freqvec);
129  fft.ClearFlag(fft.HalfSpectrum);
130 
131  for (size_t i = 0; i < N; ++i)
132  ac_tmp(i) /= (N - i);
133 
134  ac = ac_tmp.head(N).array() / ac_tmp(0);
135 }
136 
153 template <typename T>
154 void autocorrelation(const std::vector<T>& y, std::vector<T>& ac) {
155  Eigen::FFT<T> fft;
156  size_t N = y.size();
157  ac.resize(N);
158 
159  const Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, 1> > y_map(&y[0], N);
160  Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1> > ac_map(&ac[0], N);
161  autocorrelation<T>(y_map, ac_map, fft);
162 }
163 
180 template <typename T, typename DerivedA, typename DerivedB>
181 void autocorrelation(const Eigen::MatrixBase<DerivedA>& y,
182  Eigen::MatrixBase<DerivedB>& ac) {
183  Eigen::FFT<T> fft;
184  autocorrelation(y, ac, fft);
185 }
186 
187 } // namespace math
188 } // namespace stan
189 #endif
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:33
size_t fft_next_good_size(size_t N)
Find the optimal next size for the FFT so that a minimum number of zeros are padded.
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:21
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–2018, Stan Development Team.