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
Generated on Sun Sep 14 18:52:28 2008 for IT++ by Doxygen 1.5.6