43 vec
xcorr_old(
const vec &x,
const int max_lag,
const std::string scaleopt)
50 vec
xcorr(
const vec &x,
const int max_lag,
const std::string scaleopt)
52 cvec out(2*x.length() - 1);
58 cvec
xcorr(
const cvec &x,
const int max_lag,
const std::string scaleopt)
60 cvec out(2*x.length() - 1);
61 xcorr(x, x, out, max_lag, scaleopt,
true);
66 vec
xcorr(
const vec &x,
const vec &y,
const int max_lag,
const std::string scaleopt)
68 cvec out(2*x.length() - 1);
74 cvec
xcorr(
const cvec &x,
const cvec &y,
const int max_lag,
const std::string scaleopt)
76 cvec out(2*x.length() - 1);
77 xcorr(x, y, out, max_lag, scaleopt,
false);
82 void xcorr(
const vec &x,
const vec &y, vec &out,
const int max_lag,
const std::string scaleopt)
87 xcorr(xx, yy, oo, max_lag, scaleopt,
false);
92 void xcorr_old(
const vec &x,
const vec &y, vec &out,
const int max_lag,
const std::string scaleopt)
95 double s_plus, s_minus, M_double, xcorr_0, coeff_scale = 0.0;
100 M_double = double(M);
109 out.set_size(2*N - 1,
false);
111 it_assert(N <=
std::max(x.size(), y.size()),
"max_lag cannot be as large as, or larger than, the maximum length of x and y.");
113 if (scaleopt ==
"coeff") {
117 for (m = 0; m < N; m++) {
121 for (n = 0;n < M - m;n++) {
126 if (m == 0) { xcorr_0 = s_plus; }
128 if (scaleopt ==
"none") {
129 out(N + m - 1) = s_plus;
130 out(N - m - 1) = s_minus;
132 else if (scaleopt ==
"biased") {
133 out(N + m - 1) = s_plus / M_double;
134 out(N - m - 1) = s_minus / M_double;
136 else if (scaleopt ==
"unbiased") {
137 out(N + m - 1) = s_plus / double(M - m);
138 out(N - m - 1) = s_minus / double(M - m);
140 else if (scaleopt ==
"coeff") {
141 out(N + m - 1) = s_plus / coeff_scale;
142 out(N - m - 1) = s_minus / coeff_scale;
145 it_error(
"Incorrect scaleopt specified.");
150 vec
xcorr_old(
const vec &x,
const vec &y,
const int max_lag,
const std::string scaleopt)
158 void xcorr(
const cvec &x,
const cvec &y, cvec &out,
const int max_lag,
const std::string scaleopt,
bool autoflag)
160 int N =
std::max(x.length(), y.length());
164 int fftsize =
pow2i(b);
166 int end = fftsize - 1;
169 if (autoflag ==
true) {
187 if ((max_lag == -1) || (max_lag >= N))
195 out.set_size(1,
false);
199 out =
concat(temp2(end - maxlag + 1, end), temp2(0, maxlag));
203 if (scaleopt ==
"biased")
205 out = out /
static_cast<std::complex<double>
>(N);
206 else if (scaleopt ==
"unbiased") {
208 vec lags =
linspace(-maxlag, maxlag, 2 * maxlag + 1);
209 cvec scale =
to_cvec(static_cast<double>(N) -
abs(lags));
212 else if (scaleopt ==
"coeff") {
213 if (autoflag ==
true)
221 else if (scaleopt ==
"none") {}
223 it_warning(
"Unknow scaling option in XCORR, defaulting to <none> ");
228 mat
cov(
const mat &X,
bool is_zero_mean)
230 int d = X.cols(), n = X.rows();
231 mat R(d, d), m2(n, d);
238 for (
int i = 0; i < d; i++) {
240 m2.set_col(i, tmp -
mean(tmp));
244 for (
int i = 0; i < d; i++) {
245 for (
int j = 0; j <= i; j++) {
246 for (
int k = 0; k < n; k++) {
247 R(i, j) += m2(k, i) * m2(k, j);
255 for (
int i = 0; i < d; i++) {
256 for (
int j = 0; j <= i; j++) {
257 for (
int k = 0; k < n; k++) {
258 R(i, j) += X(k, i) * X(k, j);
272 "nfft must be a power of two in spectrum()!");
274 vec P(nfft / 2 + 1), w(nfft), wd(nfft);
278 double w_energy = nfft == 1 ? 1 : (nfft + 1) * .375;
280 if (nfft > v.size()) {
285 int k = (v.size() - noverlap) / (nfft - noverlap), idx = 0;
286 for (
int i = 0; i < k; i++) {
287 wd =
elem_mult(v(idx, idx + nfft - 1), w);
289 idx += nfft - noverlap;
294 P.set_size(nfft / 2 + 1,
true);
298 vec
spectrum(
const vec &v,
const vec &w,
int noverlap)
302 "The window size must be a power of two in spectrum()!");
304 vec P(nfft / 2 + 1), wd(nfft);
307 double w_energy =
energy(w);
309 if (nfft > v.size()) {
314 int k = (v.size() - noverlap) / (nfft - noverlap), idx = 0;
315 for (
int i = 0; i < k; i++) {
316 wd =
elem_mult(v(idx, idx + nfft - 1), w);
318 idx += nfft - noverlap;
323 P.set_size(nfft / 2 + 1,
true);
330 s.set_size(nfft / 2 + 1,
true);
337 s.set_size(nfft / 2 + 1,
true);