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

Generated on Sat Apr 19 10:59:23 2008 for IT++ by Doxygen 1.5.5