IT++ Logo

sigfun.cpp

Go to the documentation of this file.
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 
SourceForge Logo

Generated on Fri May 1 11:09:19 2009 for IT++ by Doxygen 1.5.8