IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
sigfun.cpp
Go to the documentation of this file.
1 
29 #include <itpp/signal/sigfun.h>
30 #include <itpp/signal/transforms.h>
31 #include <itpp/signal/window.h>
32 #include <itpp/base/converters.h>
34 #include <itpp/base/matfunc.h>
35 #include <itpp/base/specmat.h>
36 #include <itpp/base/itcompat.h>
37 #include <itpp/stat/misc_stat.h>
38 
39 
40 namespace itpp
41 {
42 
43 vec xcorr_old(const vec &x, const int max_lag, const std::string scaleopt)
44 {
45  vec out;
46  xcorr_old(x, x, out, max_lag, scaleopt);
47  return out;
48 }
49 
50 vec xcorr(const vec &x, const int max_lag, const std::string scaleopt)
51 {
52  cvec out(2*x.length() - 1); //Initial size does ont matter, it will get adjusted
53  xcorr(to_cvec(x), to_cvec(x), out, max_lag, scaleopt, true);
54 
55  return real(out);
56 }
57 
58 cvec xcorr(const cvec &x, const int max_lag, const std::string scaleopt)
59 {
60  cvec out(2*x.length() - 1); //Initial size does ont matter, it will get adjusted
61  xcorr(x, x, out, max_lag, scaleopt, true);
62 
63  return out;
64 }
65 
66 vec xcorr(const vec &x, const vec &y, const int max_lag, const std::string scaleopt)
67 {
68  cvec out(2*x.length() - 1); //Initial size does ont matter, it will get adjusted
69  xcorr(to_cvec(x), to_cvec(y), out, max_lag, scaleopt, false);
70 
71  return real(out);
72 }
73 
74 cvec xcorr(const cvec &x, const cvec &y, const int max_lag, const std::string scaleopt)
75 {
76  cvec out(2*x.length() - 1); //Initial size does ont matter, it will get adjusted
77  xcorr(x, y, out, max_lag, scaleopt, false);
78 
79  return out;
80 }
81 
82 void xcorr(const vec &x, const vec &y, vec &out, const int max_lag, const std::string scaleopt)
83 {
84  cvec xx = to_cvec(x);
85  cvec yy = to_cvec(y);
86  cvec oo = to_cvec(out);
87  xcorr(xx, yy, oo, max_lag, scaleopt, false);
88 
89  out = real(oo);
90 }
91 
92 void xcorr_old(const vec &x, const vec &y, vec &out, const int max_lag, const std::string scaleopt)
93 {
94  int m, n;
95  double s_plus, s_minus, M_double, xcorr_0, coeff_scale = 0.0;
96  int M, N;
97 
98  M = x.size();
99  M = std::max(x.size(), y.size());
100  M_double = double(M);
101 
102  if (max_lag == -1) {
103  N = std::max(x.size(), y.size());
104  }
105  else {
106  N = max_lag + 1;
107  }
108 
109  out.set_size(2*N - 1, false);
110 
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.");
112 
113  if (scaleopt == "coeff") {
114  coeff_scale = std::sqrt(energy(x)) * std::sqrt(energy(y));
115  }
116 
117  for (m = 0; m < N; m++) {
118  s_plus = 0;
119  s_minus = 0;
120 
121  for (n = 0;n < M - m;n++) {
122  s_minus += index_zero_pad(x, n) * index_zero_pad(y, n + m);
123  s_plus += index_zero_pad(x, n + m) * index_zero_pad(y, n);
124  }
125 
126  if (m == 0) { xcorr_0 = s_plus; }
127 
128  if (scaleopt == "none") {
129  out(N + m - 1) = s_plus;
130  out(N - m - 1) = s_minus;
131  }
132  else if (scaleopt == "biased") {
133  out(N + m - 1) = s_plus / M_double;
134  out(N - m - 1) = s_minus / M_double;
135  }
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);
139  }
140  else if (scaleopt == "coeff") {
141  out(N + m - 1) = s_plus / coeff_scale;
142  out(N - m - 1) = s_minus / coeff_scale;
143  }
144  else
145  it_error("Incorrect scaleopt specified.");
146  }
147 }
148 
149 
150 vec xcorr_old(const vec &x, const vec &y, const int max_lag, const std::string scaleopt)
151 {
152  vec out;
153  xcorr_old(x, y, out, max_lag, scaleopt);
154  return out;
155 }
156 
157 //Correlation
158 void xcorr(const cvec &x, const cvec &y, cvec &out, const int max_lag, const std::string scaleopt, bool autoflag)
159 {
160  int N = std::max(x.length(), y.length());
161 
162  //Compute the FFT size as the "next power of 2" of the input vector's length (max)
163  int b = ceil_i(::log2(2.0 * N - 1));
164  int fftsize = pow2i(b);
165 
166  int end = fftsize - 1;
167 
168  cvec temp2;
169  if (autoflag == true) {
170  //Take FFT of input vector
171  cvec X = fft(zero_pad(x, fftsize));
172 
173  //Compute the abs(X).^2 and take the inverse FFT.
174  temp2 = ifft(elem_mult(X, conj(X)));
175  }
176  else {
177  //Take FFT of input vectors
178  cvec X = fft(zero_pad(x, fftsize));
179  cvec Y = fft(zero_pad(y, fftsize));
180 
181  //Compute the crosscorrelation
182  temp2 = ifft(elem_mult(X, conj(Y)));
183  }
184 
185  // Compute the total number of lags to keep. We truncate the maximum number of lags to N-1.
186  int maxlag;
187  if ((max_lag == -1) || (max_lag >= N))
188  maxlag = N - 1;
189  else
190  maxlag = max_lag;
191 
192 
193  //Move negative lags to the beginning of the vector. Drop extra values from the FFT/IFFt
194  if (maxlag == 0) {
195  out.set_size(1, false);
196  out = temp2(0);
197  }
198  else
199  out = concat(temp2(end - maxlag + 1, end), temp2(0, maxlag));
200 
201 
202  //Scale data
203  if (scaleopt == "biased")
204  //out = out / static_cast<double_complex>(N);
205  out = out / static_cast<std::complex<double> >(N);
206  else if (scaleopt == "unbiased") {
207  //Total lag vector
208  vec lags = linspace(-maxlag, maxlag, 2 * maxlag + 1);
209  cvec scale = to_cvec(static_cast<double>(N) - abs(lags));
210  out /= scale;
211  }
212  else if (scaleopt == "coeff") {
213  if (autoflag == true) // Normalize by Rxx(0)
214  out /= out(maxlag);
215  else { //Normalize by sqrt(Rxx(0)*Ryy(0))
216  double rxx0 = sum(abs(elem_mult(x, x)));
217  double ryy0 = sum(abs(elem_mult(y, y)));
218  out /= std::sqrt(rxx0 * ryy0);
219  }
220  }
221  else if (scaleopt == "none") {}
222  else
223  it_warning("Unknow scaling option in XCORR, defaulting to <none> ");
224 
225 }
226 
227 
228 mat cov(const mat &X, bool is_zero_mean)
229 {
230  int d = X.cols(), n = X.rows();
231  mat R(d, d), m2(n, d);
232  vec tmp;
233 
234  R = 0.0;
235 
236  if (!is_zero_mean) {
237  // Compute and remove mean
238  for (int i = 0; i < d; i++) {
239  tmp = X.get_col(i);
240  m2.set_col(i, tmp - mean(tmp));
241  }
242 
243  // Calc corr matrix
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);
248  }
249  R(j, i) = R(i, j); // When i=j this is unnecassary work
250  }
251  }
252  }
253  else {
254  // Calc corr matrix
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);
259  }
260  R(j, i) = R(i, j); // When i=j this is unnecassary work
261  }
262  }
263  }
264  R /= n;
265 
266  return R;
267 }
268 
269 vec spectrum(const vec &v, int nfft, int noverlap)
270 {
271  it_assert_debug(pow2i(levels2bits(nfft)) == nfft,
272  "nfft must be a power of two in spectrum()!");
273 
274  vec P(nfft / 2 + 1), w(nfft), wd(nfft);
275 
276  P = 0.0;
277  w = hanning(nfft);
278  double w_energy = nfft == 1 ? 1 : (nfft + 1) * .375; // Hanning energy
279 
280  if (nfft > v.size()) {
281  P = sqr(abs(fft(to_cvec(elem_mult(zero_pad(v, nfft), w)))(0, nfft / 2)));
282  P /= w_energy;
283  }
284  else {
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);
288  P += sqr(abs(fft(to_cvec(wd))(0, nfft / 2)));
289  idx += nfft - noverlap;
290  }
291  P /= k * w_energy;
292  }
293 
294  P.set_size(nfft / 2 + 1, true);
295  return P;
296 }
297 
298 vec spectrum(const vec &v, const vec &w, int noverlap)
299 {
300  int nfft = w.size();
301  it_assert_debug(pow2i(levels2bits(nfft)) == nfft,
302  "The window size must be a power of two in spectrum()!");
303 
304  vec P(nfft / 2 + 1), wd(nfft);
305 
306  P = 0.0;
307  double w_energy = energy(w);
308 
309  if (nfft > v.size()) {
310  P = sqr(abs(fft(to_cvec(elem_mult(zero_pad(v, nfft), w)))(0, nfft / 2)));
311  P /= w_energy;
312  }
313  else {
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);
317  P += sqr(abs(fft(to_cvec(wd))(0, nfft / 2)));
318  idx += nfft - noverlap;
319  }
320  P /= k * w_energy;
321  }
322 
323  P.set_size(nfft / 2 + 1, true);
324  return P;
325 }
326 
327 vec filter_spectrum(const vec &a, int nfft)
328 {
329  vec s = sqr(abs(fft(to_cvec(a), nfft)));
330  s.set_size(nfft / 2 + 1, true);
331  return s;
332 }
333 
334 vec filter_spectrum(const vec &a, const vec &b, int nfft)
335 {
336  vec s = sqr(abs(elem_div(fft(to_cvec(a), nfft), fft(to_cvec(b), nfft))));
337  s.set_size(nfft / 2 + 1, true);
338  return s;
339 }
340 
341 } // namespace itpp
342 
343 
344 
345 
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:14 for IT++ by Doxygen 1.8.1.2