00001 00030 #include <itpp/signal/sigfun.h> 00031 #include <itpp/signal/transforms.h> 00032 #include <itpp/signal/window.h> 00033 #include <itpp/base/converters.h> 00034 #include <itpp/base/math/elem_math.h> 00035 #include <itpp/base/matfunc.h> 00036 #include <itpp/base/specmat.h> 00037 #include <itpp/stat/misc_stat.h> 00038 00039 00040 namespace itpp 00041 { 00042 00043 vec xcorr_old(const vec &x, const int max_lag, const std::string scaleopt) 00044 { 00045 vec out; 00046 xcorr_old(x, x, out, max_lag, scaleopt); 00047 return out; 00048 } 00049 00050 vec xcorr(const vec &x, const int max_lag, const std::string scaleopt) 00051 { 00052 cvec out(2*x.length() - 1); //Initial size does ont matter, it will get adjusted 00053 xcorr(to_cvec(x), to_cvec(x), out, max_lag, scaleopt, true); 00054 00055 return real(out); 00056 } 00057 00058 cvec xcorr(const cvec &x, const int max_lag, const std::string scaleopt) 00059 { 00060 cvec out(2*x.length() - 1); //Initial size does ont matter, it will get adjusted 00061 xcorr(x, x, out, max_lag, scaleopt, true); 00062 00063 return out; 00064 } 00065 00066 vec xcorr(const vec &x, const vec &y, const int max_lag, const std::string scaleopt) 00067 { 00068 cvec out(2*x.length() - 1); //Initial size does ont matter, it will get adjusted 00069 xcorr(to_cvec(x), to_cvec(y), out, max_lag, scaleopt, false); 00070 00071 return real(out); 00072 } 00073 00074 cvec xcorr(const cvec &x, const cvec &y, const int max_lag, const std::string scaleopt) 00075 { 00076 cvec out(2*x.length() - 1); //Initial size does ont matter, it will get adjusted 00077 xcorr(x, y, out, max_lag, scaleopt, false); 00078 00079 return out; 00080 } 00081 00082 void xcorr(const vec &x, const vec &y, vec &out, const int max_lag, const std::string scaleopt) 00083 { 00084 cvec xx = to_cvec(x); 00085 cvec yy = to_cvec(y); 00086 cvec oo = to_cvec(out); 00087 xcorr(xx, yy, oo, max_lag, scaleopt, false); 00088 00089 out = real(oo); 00090 } 00091 00092 void xcorr_old(const vec &x, const vec &y, vec &out, const int max_lag, const std::string scaleopt) 00093 { 00094 int m, n; 00095 double s_plus, s_minus, M_double, xcorr_0, coeff_scale = 0.0; 00096 int M, N; 00097 00098 M = x.size(); 00099 M = std::max(x.size(), y.size()); 00100 M_double = double(M); 00101 00102 if (max_lag == -1) { 00103 N = std::max(x.size(), y.size()); 00104 } 00105 else { 00106 N = max_lag + 1; 00107 } 00108 00109 out.set_size(2*N - 1, false); 00110 00111 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."); 00112 00113 if (scaleopt == "coeff") { 00114 coeff_scale = std::sqrt(energy(x)) * std::sqrt(energy(y)); 00115 } 00116 00117 for (m = 0; m < N; m++) { 00118 s_plus = 0; 00119 s_minus = 0; 00120 00121 for (n = 0;n < M - m;n++) { 00122 s_minus += index_zero_pad(x, n) * index_zero_pad(y, n + m); 00123 s_plus += index_zero_pad(x, n + m) * index_zero_pad(y, n); 00124 } 00125 00126 if (m == 0) { xcorr_0 = s_plus; } 00127 00128 if (scaleopt == "none") { 00129 out(N + m - 1) = s_plus; 00130 out(N - m - 1) = s_minus; 00131 } 00132 else if (scaleopt == "biased") { 00133 out(N + m - 1) = s_plus / M_double; 00134 out(N - m - 1) = s_minus / M_double; 00135 } 00136 else if (scaleopt == "unbiased") { 00137 out(N + m - 1) = s_plus / double(M - m); 00138 out(N - m - 1) = s_minus / double(M - m); 00139 } 00140 else if (scaleopt == "coeff") { 00141 out(N + m - 1) = s_plus / coeff_scale; 00142 out(N - m - 1) = s_minus / coeff_scale; 00143 } 00144 else 00145 it_error("Incorrect scaleopt specified."); 00146 } 00147 } 00148 00149 00150 vec xcorr_old(const vec &x, const vec &y, const int max_lag, const std::string scaleopt) 00151 { 00152 vec out; 00153 xcorr_old(x, y, out, max_lag, scaleopt); 00154 return out; 00155 } 00156 00157 //Correlation 00158 void xcorr(const cvec &x, const cvec &y, cvec &out, const int max_lag, const std::string scaleopt, bool autoflag) 00159 { 00160 int N = std::max(x.length(), y.length()); 00161 00162 //Compute the FFT size as the "next power of 2" of the input vector's length (max) 00163 int b = ceil_i(::log2(2.0 * N - 1)); 00164 int fftsize = pow2i(b); 00165 00166 int end = fftsize - 1; 00167 00168 cvec temp2; 00169 if (autoflag == true) { 00170 //Take FFT of input vector 00171 cvec X = fft(zero_pad(x, fftsize)); 00172 00173 //Compute the abs(X).^2 and take the inverse FFT. 00174 temp2 = ifft(elem_mult(X, conj(X))); 00175 } 00176 else { 00177 //Take FFT of input vectors 00178 cvec X = fft(zero_pad(x, fftsize)); 00179 cvec Y = fft(zero_pad(y, fftsize)); 00180 00181 //Compute the crosscorrelation 00182 temp2 = ifft(elem_mult(X, conj(Y))); 00183 } 00184 00185 // Compute the total number of lags to keep. We truncate the maximum number of lags to N-1. 00186 int maxlag; 00187 if ((max_lag == -1) || (max_lag >= N)) 00188 maxlag = N - 1; 00189 else 00190 maxlag = max_lag; 00191 00192 00193 //Move negative lags to the beginning of the vector. Drop extra values from the FFT/IFFt 00194 if (maxlag == 0) { 00195 out.set_size(1, false); 00196 out = temp2(0); 00197 } 00198 else 00199 out = concat(temp2(end - maxlag + 1, end), temp2(0, maxlag)); 00200 00201 00202 //Scale data 00203 if (scaleopt == "biased") 00204 //out = out / static_cast<double_complex>(N); 00205 out = out / static_cast<std::complex<double> >(N); 00206 else if (scaleopt == "unbiased") { 00207 //Total lag vector 00208 vec lags = linspace(-maxlag, maxlag, 2 * maxlag + 1); 00209 cvec scale = to_cvec(static_cast<double>(N) - abs(lags)); 00210 out /= scale; 00211 } 00212 else if (scaleopt == "coeff") { 00213 if (autoflag == true) // Normalize by Rxx(0) 00214 out /= out(maxlag); 00215 else { //Normalize by sqrt(Rxx(0)*Ryy(0)) 00216 double rxx0 = sum(abs(elem_mult(x, x))); 00217 double ryy0 = sum(abs(elem_mult(y, y))); 00218 out /= std::sqrt(rxx0 * ryy0); 00219 } 00220 } 00221 else if (scaleopt == "none") {} 00222 else 00223 it_warning("Unknow scaling option in XCORR, defaulting to <none> "); 00224 00225 } 00226 00227 00228 mat cov(const mat &X, bool is_zero_mean) 00229 { 00230 int d = X.cols(), n = X.rows(); 00231 mat R(d, d), m2(n, d); 00232 vec tmp; 00233 00234 R = 0.0; 00235 00236 if (!is_zero_mean) { 00237 // Compute and remove mean 00238 for (int i = 0; i < d; i++) { 00239 tmp = X.get_col(i); 00240 m2.set_col(i, tmp - mean(tmp)); 00241 } 00242 00243 // Calc corr matrix 00244 for (int i = 0; i < d; i++) { 00245 for (int j = 0; j <= i; j++) { 00246 for (int k = 0; k < n; k++) { 00247 R(i, j) += m2(k, i) * m2(k, j); 00248 } 00249 R(j, i) = R(i, j); // When i=j this is unnecassary work 00250 } 00251 } 00252 } 00253 else { 00254 // Calc corr matrix 00255 for (int i = 0; i < d; i++) { 00256 for (int j = 0; j <= i; j++) { 00257 for (int k = 0; k < n; k++) { 00258 R(i, j) += X(k, i) * X(k, j); 00259 } 00260 R(j, i) = R(i, j); // When i=j this is unnecassary work 00261 } 00262 } 00263 } 00264 R /= n; 00265 00266 return R; 00267 } 00268 00269 vec spectrum(const vec &v, int nfft, int noverlap) 00270 { 00271 it_assert_debug(pow2i(levels2bits(nfft)) == nfft, 00272 "nfft must be a power of two in spectrum()!"); 00273 00274 vec P(nfft / 2 + 1), w(nfft), wd(nfft); 00275 00276 P = 0.0; 00277 w = hanning(nfft); 00278 double w_energy = nfft == 1 ? 1 : (nfft + 1) * .375; // Hanning energy 00279 00280 if (nfft > v.size()) { 00281 P = sqr(abs(fft(to_cvec(elem_mult(zero_pad(v, nfft), w)))(0, nfft / 2))); 00282 P /= w_energy; 00283 } 00284 else { 00285 int k = (v.size() - noverlap) / (nfft - noverlap), idx = 0; 00286 for (int i = 0; i < k; i++) { 00287 wd = elem_mult(v(idx, idx + nfft - 1), w); 00288 P += sqr(abs(fft(to_cvec(wd))(0, nfft / 2))); 00289 idx += nfft - noverlap; 00290 } 00291 P /= k * w_energy; 00292 } 00293 00294 P.set_size(nfft / 2 + 1, true); 00295 return P; 00296 } 00297 00298 vec spectrum(const vec &v, const vec &w, int noverlap) 00299 { 00300 int nfft = w.size(); 00301 it_assert_debug(pow2i(levels2bits(nfft)) == nfft, 00302 "The window size must be a power of two in spectrum()!"); 00303 00304 vec P(nfft / 2 + 1), wd(nfft); 00305 00306 P = 0.0; 00307 double w_energy = energy(w); 00308 00309 if (nfft > v.size()) { 00310 P = sqr(abs(fft(to_cvec(elem_mult(zero_pad(v, nfft), w)))(0, nfft / 2))); 00311 P /= w_energy; 00312 } 00313 else { 00314 int k = (v.size() - noverlap) / (nfft - noverlap), idx = 0; 00315 for (int i = 0; i < k; i++) { 00316 wd = elem_mult(v(idx, idx + nfft - 1), w); 00317 P += sqr(abs(fft(to_cvec(wd))(0, nfft / 2))); 00318 idx += nfft - noverlap; 00319 } 00320 P /= k * w_energy; 00321 } 00322 00323 P.set_size(nfft / 2 + 1, true); 00324 return P; 00325 } 00326 00327 vec filter_spectrum(const vec &a, int nfft) 00328 { 00329 vec s = sqr(abs(fft(to_cvec(a), nfft))); 00330 s.set_size(nfft / 2 + 1, true); 00331 return s; 00332 } 00333 00334 vec filter_spectrum(const vec &a, const vec &b, int nfft) 00335 { 00336 vec s = sqr(abs(elem_div(fft(to_cvec(a), nfft), fft(to_cvec(b), nfft)))); 00337 s.set_size(nfft / 2 + 1, true); 00338 return s; 00339 } 00340 00341 } // namespace itpp 00342 00343 00344 00345
Generated on Sun Jul 26 08:36:51 2009 for IT++ by Doxygen 1.5.9