IT++ Logo Newcom Logo

sigfun.cpp

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

Generated on Thu Apr 19 14:14:57 2007 for IT++ by Doxygen 1.5.1