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
Generated on Thu Apr 19 14:18:29 2007 for IT++ by Doxygen 1.5.1