22 template<class DataType_t>
34 template<
class DataType_t>
47 template<
class DataType_t>
60 template<
class DataType_t>
75 template<
class DataType_t>
78 if(this->pWavelet) this->pWavelet->release();
79 if(this->pWavelet)
delete this->pWavelet;
84 template<
class DataType_t>
89 if(pWavelet->allocate())
90 maxlevel = pWavelet->getMaxLevel();
98 template<
class DataType_t>
105 int I = maxLayer()+1;
106 double df = this->
rate()/I/4.;
107 if(I==1)
return this->
rate()/2.;
108 if(pWavelet->m_WaveType==
WDMT)
return i*this->
rate()/(I-1)/2.;
109 if(pWavelet->BinaryTree())
return 2*df*(i+0.5);
110 df = this->
rate()/(1<<I);
112 while(i--){ f += 2*
df; df*=2; }
116 template<
class DataType_t>
120 int I = this->maxLayer()+1;
133 if(j==J) {error =
true;
break;}
138 cout<<
"wseries::mul(): "<<x.
size()<<
" "<<y.
size()<<endl;
142 this->putLayer(y,i++);
144 if(error) {cout<<
"wseries::mul() error!\n";
exit(1);}
148 template<
class DataType_t>
156 int I =
int(maxLayer()+1);
157 if(I<2 || f>=this->
rate()/2.)
return I;
158 double df = this->
rate()/(I-1)/2.;
159 if(pWavelet->m_WaveType==
WDMT)
return int((f+df/2.)/df);
160 df = this->
rate()/I/2.;
161 if(pWavelet->BinaryTree())
return int(f/df);
162 df = this->
rate()/(1<<I);
165 while(ff<f){ ff += 2*
df; df*=2; i++; }
174 template<
class DataType_t>
178 if(n > maxLayer()) index = maxLayer();
179 slice s = pWavelet->getSlice(index);
181 if(this->limit(s) <= this->
size()){
183 value.
rate(this->wrate());
186 value.
edge(this->edge());
192 cout<<
"WSeries::getLayer(): data length mismatch: "<<this->limit(s)<<
" "<<this->
size()<<
"\n";
200 template<
class DataType_t>
203 slice s = this->pWavelet->getSlice(index);
205 if( (s.
size() < value.
size()) || (this->limit(s) > this->
size()) ){
206 cout<<
"WSeries::putLayer(): invalid array size.\n";
215 template<
class DataType_t>
224 pWavelet->allocate(this->
size(), this->
data);
227 template<
class DataType_t>
230 if(pWavelet->allocate()){
231 pWavelet->nSTS = pWavelet->nWWS;
233 if(pWavelet->pWWS != this->data || pWavelet->nWWS != this->Size){
234 this->
data = pWavelet->pWWS;
235 this->Size = pWavelet->nWWS;
242 throw std::invalid_argument
243 (
"WSeries::Forward(): data is not allocated");
247 template<
class DataType_t>
251 if(pWavelet->allocate()) pWavelet->release();
253 this->wrate(x.
rate());
254 f_high = x.
rate()/2.;
255 pWavelet->allocate(this->
size(), this->
data);
260 template<
class DataType_t>
264 if(pWavelet->allocate()) pWavelet->release();
266 this->wrate(x.
rate());
267 f_high = x.
rate()/2.;
272 template<
class DataType_t>
275 if(pWavelet->allocate()){
277 if(pWavelet->pWWS != this->data || pWavelet->nWWS != this->Size) {
278 this->
data = pWavelet->pWWS;
279 this->Size = pWavelet->nWWS;
281 this->wrate(this->
rate());
289 throw std::invalid_argument
290 (
"WSeries::Inverse(): data is not allocated");
294 template<
class DataType_t>
302 if(!this->pWavelet) {
303 cout<<
"WSeries::bandpass ERROR: no transformation is specified"<<endl;
311 size_t I = this->maxLayer()+1;
312 for(
size_t i=0;
i<I;
i++) {
314 if(freq<flow || freq>fhigh) {
316 this->putLayer(wa,
i);
321 template<
class DataType_t>
334 double fl =
fabs(f1)>0. ?
fabs(f1) : this->getlow();
335 double fh =
fabs(f2)>0. ?
fabs(f2) : this->gethigh();
336 size_t n = this->pWavelet->m_WaveType==
WDMT ? size_t((fl+dF/2.)/dF+0.1) : size_t(fl/dF+0.1);
337 size_t m = this->pWavelet->m_WaveType==
WDMT ? size_t((fh+dF/2.)/dF+0.1)-1 : size_t(fh/dF+0.1)-1;
338 size_t M = this->maxLayer()+1;
343 for(i=0; i<
int(M); i++) {
345 if((f1>=0 && i>n) && (f2>=0 && i<=
m))
continue;
346 if((f1<0 && i<n) || (f2<0 && i>
m))
continue;
347 if((f1<0 && f2>=0 && i<n))
continue;
348 if((f1>=0 && f2<0 && i>=m))
continue;
350 this->
getLayer(w,i+0.01); w=
a; this->putLayer(w,i+0.01);
351 this->
getLayer(w,-i-0.01); w=
a; this->putLayer(w,-i-0.01);
357 template<
class DataType_t>
381 double aa,ee,EE,uu,UU,dd,DD,em,
ss,cc,nn;
382 double shape,mean,alp;
384 int J = this->
size()/2;
385 if(pattern==0)
return 1.;
386 if(!this->isWDM() || (2*J)/this->xsize()==1)
exit(0);
393 int p[] = {0,0,0,0,0,0,0,0,0};
398 if(c!=
'a'||c!=
'A') this->
resize(J);
405 else if(pattern==2) {
409 else if(pattern==3) {
414 else if(pattern==4) {
419 else if(pattern==5) {
420 p[1]=M+1; p[2]=-M-1; p[3]=2*M+2; p[4]=-2*M-2;
424 else if(pattern==6) {
425 p[1]=-M+1; p[2]=M-1; p[3]=-2*M+2; p[4]=2*M-2;
429 else if(pattern==7) {
430 p[1]=1; p[2]=-1; p[3]=
M; p[4]=-
M;
434 else if(pattern==8) {
435 p[1]=M+1; p[2]=-M+1; p[3]= M-1; p[4]=-M-1;
439 else if(pattern==9) {
440 p[1]=1; p[2]=-1; p[3]=
M; p[4]=-
M;
441 p[5]=M+1; p[6]=M-1; p[7]=-M+1; p[8]=-M-1;
444 }
else { shape = mean = 1.; }
449 for(
int j=jb;
j<je;
j++) {
451 if(m<mL || m>mH) {this->
data[
j]=0.;
continue;}
452 q = in.
data+
j; Q = q+J;
453 ss = q[p[1]]*Q[p[1]]+q[p[2]]*Q[p[2]]+q[p[3]]*Q[p[3]]+q[p[4]]*Q[p[4]]
454 + q[p[5]]*Q[p[5]]+q[p[6]]*Q[p[6]]+q[p[7]]*Q[p[7]]+q[p[8]]*Q[p[8]];
455 ee = q[p[1]]*q[p[1]]+q[p[2]]*q[p[2]]+q[p[3]]*q[p[3]]+q[p[4]]*q[p[4]]
456 + q[p[5]]*q[p[5]]+q[p[6]]*q[p[6]]+q[p[7]]*q[p[7]]+q[p[8]]*q[p[8]];
457 EE = Q[p[1]]*Q[p[1]]+Q[p[2]]*Q[p[2]]+Q[p[3]]*Q[p[3]]+Q[p[4]]*Q[p[4]]
458 + Q[p[5]]*Q[p[5]]+Q[p[6]]*Q[p[6]]+Q[p[7]]*Q[p[7]]+Q[p[8]]*Q[p[8]];
459 ss+= q[p[0]]*Q[p[0]]*(mean-8);
460 ee+= q[p[0]]*q[p[0]]*(mean-8);
461 EE+= Q[p[0]]*Q[p[0]]*(mean-8);
464 nn = sqrt(cc*cc+ss*ss);
465 if(ee+EE<nn) nn=ee+EE;
466 aa = sqrt((ee+EE+nn)/2) + sqrt((ee+EE-nn)/2);
467 em = (c==
'e'||c==
'l'||mean==1.) ? (ee+EE)/2. : aa*aa/4;
468 alp = shape-
log(shape)/3;
471 if(em<alp) {em=0;
continue;}
472 em-= alp*(1+
log(em/alp));
476 nn = sqrt((1+cc)*(1+cc)+ss*ss)/2;
477 this->
data[
j]=aa*cc/nn; this->
data[
j+J]=aa*ss/nn/2;
478 }
else {this->
data[
j] = em;}
479 if(hist && em>0.01) hist->Fill(sqrt(em));
485 template<
class DataType_t>
501 cout<<
"wseries::maxEnergy(): illegal wavelet\n";
510 tmp.
setlow(this->getlow());
514 for(
int k=N;
k<=
K;
k+=
N) {
517 tmp.
setlow(this->getlow());
523 tmp.
setlow(this->getlow());
529 for(
int k=N;
k<=
K;
k+=
N) {
541 this->putLayer(xx,0.1);
543 this->putLayer(xx,M-1);
545 int m = abs(pattern);
546 if(m==5 || m==6 || m==9) {
548 this->putLayer(xx,1);
550 this->putLayer(xx,M-2);
553 if(!pattern)
return 1.;
554 return Gamma2Gauss(hist);
557 template<
class DataType_t>
568 double fff = (nR-nL)*tmp.
wavecount(0.001)/double(nn);
569 double med = tmp.
waveSplit(nL,nR,nR-
int(0.5*fff));
570 double amp, aaa, bbb, rms, alp;
572 aaa = bbb = 0.; nn = 0;
573 for(
int i=nL;
i<nR;
i++) {
574 amp = (double)this->
data[
i];
575 if(amp>0.01 && amp<20*med) {
576 aaa += amp; bbb +=
log(amp); nn++;
579 alp =
log(aaa/nn)-bbb/nn;
580 alp = (3-alp+sqrt((alp-3)*(alp-3)+24*alp))/12./alp;
581 double avr = med*(3*alp+0.2)/(3*alp-0.8);
584 double ALP = med*alp/
avr;
585 for(
int i=0;
i<this->
size();
i++) {
586 amp = (double)this->
data[
i]*alp/avr;
587 if(amp<ALP) {this->
data[
i]=0.;
continue;}
588 this->
data[
i] = amp-ALP*(1+
log(amp/ALP));
593 rms = 1./tmp.
waveSplit(nL,nR,nR-
int(0.3173*fff));
595 for(
int i=0;
i<this->
size();
i++) {
596 this->
data[
i] *= rms;
597 if(hist &&
i>nL &&
i<nR) hist->Fill(sqrt(this->
data[
i]));
603 template<
class DataType_t>
614 std::vector<double>
vR;
615 std::vector<double> vF;
616 int mBand, mRate, level;
618 double time, freq,
a,
A,ee,EE,
ss,cc,gg,b,
B;
619 double mean = 2.*N-1;
620 double shape = N-
log(N)/
N;
622 mBand = 0; mRate = 0;
623 for(
int n=0;
n<
N;
n++) {
624 vF.push_back(pws[
n]->resolution());
625 vR.push_back(pws[
n]->wrate());
626 vM.push_back(pws[
n]->maxLayer()+1);
627 vJ.push_back(pws[
n]->
size()/2);
628 if(vM[
n]>mBand) {mBand = vM[
n]; nn=
n;}
629 if(pws[
n]->wrate()>mRate) mRate = pws[
n]->
wrate();
634 time = pws[nn]->
size()/pws[nn]->
wrate()/2;
635 J = mRate*
int(time+0.1);
636 cout<<
"debug1: "<<mBand<<
" "<<mRate<<
" "<<J<<
" "<<pws[nn]->
size()<<endl;
643 this->pWavelet->nWWS = J;
644 this->pWavelet->nSTS = (J/mBand)*(mBand-1);
646 int nL =
int(this->Edge*mRate*mBand);
647 int nR = this->
size()-nL-1;
649 for(
int i=0;
i<this->
size();
i++) {
650 time = double(
i/mBand)/mRate;
651 freq = (
i%mBand)*vF[nn];
653 for(
int n=1;
n<
N;
n++) {
654 j =
int(time*vR[
n-1]+0.1)*vM[
n-1];
655 j+=
int(freq/vF[
n-1]+0.1);
657 A = pws[
n]->
data[j+vJ[
n-1]];
658 j =
int(time*vR[
n]+0.1)*vM[
n];
659 j+=
int(freq/vF[
n]+0.1);
661 B = pws[
n]->
data[j+vJ[
n]];
667 gg = sqrt(cc*cc+4*ss*ss);
668 a = sqrt((ee+EE+gg)/2);
669 A = sqrt((ee+EE-gg)/2);
673 this->
data[
i] = (a+
A)*(a+A)/mean/2.;
674 if(hist &&
i>nL &&
i<nR) hist->Fill(this->
data[
i]);
683 template<
class DataType_t>
687 if( pWavelet->allocate() ) pWavelet->release();
688 if(p->
size() != a.
size()) pWavelet->reset();
692 f_high = a.
rate()/2.;
693 pWavelet->allocate(this->
size(), this->
data);
697 template<
class DataType_t>
712 template<
class DataType_t>
716 if(this->limit() > this->
size()){
717 cout <<
"WSeries::operator[]: Illegal argument: "<<this->limit()<<
" "<<this->
size()<<
"\n";
723 template<
class DataType_t>
727 template<
class DataType_t>
731 template<
class DataType_t>
735 template<
class DataType_t>
744 template<
class DataType_t>
748 template<
class DataType_t>
752 template<
class DataType_t>
761 if(pWavelet->m_TreeType != a.
pWavelet->m_TreeType){
762 cout<<
"WSeries::operator* : wavelet tree type mismatch."<<endl;
771 for(i=0; i<= max_layer; i++)
772 (*p)[pWavelet->getSlice(i)] *= (*pa)[a.
pWavelet->getSlice(i)];
777 template<
class DataType_t>
785 if(pWavelet->m_TreeType != a.
pWavelet->m_TreeType){
786 cout<<
"WSeries::operator+ : wavelet tree type mismatch."<<endl;
795 for(i=0; i<= max_layer; i++)
796 (*p)[pWavelet->getSlice(i)] += (*pa)[a.
pWavelet->getSlice(i)];
801 template<
class DataType_t>
809 if(pWavelet->m_TreeType != a.
pWavelet->m_TreeType){
810 cout<<
"WSeries::operator- : wavelet tree type mismatch."<<endl;
819 for(i=0; i<= max_layer; i++)
820 (*p)[pWavelet->getSlice(i)] -= (*pa)[a.
pWavelet->getSlice(i)];
825 template<
class DataType_t>
830 size_t max_layer = maxLayer()+1;
832 if(max_layer == a.
size()) {
833 for(i=0; i< max_layer; i++) {
834 (*p)[pWavelet->getSlice(i)] *= a.
data[
i];
841 else cout<<
"WSeries::operator* - no operation is performed"<<endl;
848 template<
class DataType_t>
853 int n = this->
size();
854 int m = this->maxLayer()+1;
856 if (app == 1)
strcpy (mode,
"a");
860 if ( (fp = fopen(fname, mode)) == NULL ) {
861 cout <<
" Dump() error: cannot open file " << fname <<
". \n";
866 fprintf( fp,
"# start time: -start %lf \n", this->Start );
867 fprintf( fp,
"# sampling rate: -rate %lf \n", this->
Rate );
868 fprintf( fp,
"# number of samples: -size %d \n", (
int)this->Size );
869 fprintf( fp,
"# number of layers: -n %d \n", m );
872 for (i = 0; i <
m; i++) {
875 for(j = 0; j <
n; j++)
fprintf( fp,
"%e ", (
float)a.
data[j]);
882 template<
class DataType_t>
886 if( pWavelet->allocate() ) pWavelet->release();
887 p->wavearray<DataType_t>::resize(n);
888 pWavelet->allocate(this->
size(), this->
data);
892 wRate = this->
rate();
893 f_high = this->
rate()/2.;
896 template<
class DataType_t>
900 if( pWavelet->allocate() ) pWavelet->release();
901 p->wavearray<DataType_t>::resample(f,nF);
902 pWavelet->allocate(this->
size(), this->
data);
906 wRate = p->wavearray<DataType_t>::rate();
907 f_high = p->wavearray<DataType_t>::rate()/2.;
910 template<
class DataType_t>
913 #if !defined (__SUNPRO_CC)
918 register float* q = NULL;
919 register float* p = NULL;
928 if(!pWavelet->BinaryTree())
return 1.;
930 int ni = maxLayer()+1;
931 int nj = this->
size()/ni;
935 bool CROSS = t<0 || f<0;
950 p[
j] = (float)x.
data[j];
951 q[j] = (
float)y.
data[
j];
963 if(p[j]==0. && q[j]==0.)
continue;
965 is = i-f<0 ? 0 : i-
f;
966 js = j-t<0 ? 0 : j-
t;
967 ie = i+f>n ? n : i+
f;
968 je = j+t>m ? m : j+
t;
972 for(u=is; u<=ie; u++)
973 for(v=js; v<=je; v++){
974 if(CROSS && !(i==u || j==v))
continue;
975 if(B[u][v]!=0.) energy +=
log(
fabs(B[u][v]));
977 if(energy < threshold) x.
data[
j]=0;
982 for(u=is; u<=ie; u++)
983 for(v=js; v<=je; v++){
984 if(CROSS && !(i==u || j==v))
continue;
985 if(A[u][v]!=0.) energy +=
log(
fabs(A[u][v]));
987 if(energy < threshold) y.
data[
j]=0;
1001 template<
class DataType_t>
1007 size_t N = a.
size();
1014 DataType_t* p = NULL;
1015 DataType_t* q = NULL;
1016 DataType_t*
P = NULL;
1017 DataType_t*
Q = NULL;
1019 if(pWavelet->m_TreeType != a.
pWavelet->m_TreeType){
1020 cout<<
"WSeries::operator- : wavelet tree type mismatch."<<endl;
1024 size_t max_layer = (maxLayer() > a.
maxLayer()) ? a.
maxLayer() : maxLayer();
1026 for(k=0; k<= max_layer; k++){
1038 if(!n && w>=0.) n++;
1047 for(i=x.
start(); i<
N; i+=step)
1049 if(this->
data[i] == 0.)
continue;
1064 if(*p>0) { E += *p; m++; }
1069 if(
gammaCL(E,m) > S-
log(
double(m))) pass =
true;
1073 else this->
data[
i] = 0;
1079 if(
size_t(maxLayer())>max_layer){
1081 for(k=max_layer+1; k<= size_t(maxLayer()); k++) {
1082 (*pw)[getSlice(k)] = 0;
1086 return double(event)/this->
size();
1090 template<
class DataType_t>
1094 int M = maxLayer()+1;
1097 this->setSlice(getSlice(i));
1107 template<
class DataType_t>
1110 if(offset<T) offset =
T;
1112 size_t M = maxLayer()+1;
1118 if(mode<2) a.
lprFilter(T,mode,stride,offset);
1119 else a.
spesla(T,stride,offset);
1127 template<
class DataType_t>
1147 double segT = this->
stop()-this->
start();
1148 if(t <= 0.) t = segT-2.*
offset;
1149 int M = maxLayer()+1;
1150 int m = mode<0 ? -1 : 1;
1152 this->w_mode = abs(mode);
1153 if(stride > t || stride<=0.) stride =
t;
1154 size_t K = size_t((segT-2.*offset)/stride);
1166 getLayer(aa,-(i+0.01)); aa*=aa; a*=
a; a+=aa;
1168 b = a.
white(t,w_mode,offset,stride);
1169 if(b.
size() != K+1) cout<<
"wseries::white(): array size mismatch\n";
1171 if(w_mode) putLayer(a,i*m);
1185 template<
class DataType_t>
1194 if(!nRMS.
size())
return false;
1199 int m = mode<0 ? -1 : 1;
1201 size_t I = this->maxLayer()+1;
1204 double R = this->wrate();
1209 this->w_mode = abs(mode);
1212 cout<<
"wseries::white error: mismatch between WSeries and nRMS\n";
1216 for(i=(1-m)/2; i<I; i++){
1223 for(j=0; j<J; j++) {
1227 if(t >= To+K*dT) r = na.
data[
K];
1228 else if(t <= To) r = na.
data[0];
1230 if(t > T) {k++; T+=
dT;}
1231 r = (na.
data[k-1]*(dT-T+
t) + na.
data[k]*(T-t))/dT;
1233 xx.
data[
j] *= abs(mode)<=1 ? DataType_t(1./r) : DataType_t(1./r/r);
1235 this->putLayer(xx,(i+0.01)*m);
1241 template<
class DataType_t>
1245 size_t M = maxLayer()+1;
1253 if(!pWavelet->BinaryTree()) { v = 1.;
return v; }
1262 v.
data[i/
k] += x>0. ? 1./x/x : 0.;
1268 for(i=0; i<v.
size(); i++)
1277 template<
class DataType_t>
1280 if(fl<0.) fl = this->getlow();
1281 if(fh<0.) fh = this->gethigh();
1282 double tsRate = this->
rate();
1285 size_t M = maxLayer()+1;
1286 size_t N = this->
size()/
M;
1287 size_t ml = size_t(2.*M*fl/tsRate+0.5);
1288 size_t mh = size_t(2.*M*fh/tsRate+0.5);
1292 size_t nL = ml+
int((mh-
int(ml))/4.+0.5);
1293 size_t nR = mh-
int((mh-
int(ml))/4.+0.5);
1294 size_t n = size_t(
fabs(t)*tsRate/M);
1305 if(!pWavelet->BinaryTree() || mh<ml+8 || nL<1)
1306 { v = 1.;
return v; }
1314 laYer[inDex[
j]] =
j;
1320 p = this->
data + i*
M;
1321 for(j=0; j<
M; j++){ pp[
j] = &(p[inDex[
j]]); }
1322 this->waveSplit(pp,ml,mh-1,nL-1);
1323 this->waveSplit(pp,nL,mh-1,nR);
1324 v.
data[
i] = float(*pp[nR] - *pp[nL-1])/2./0.6745;
1328 v.
rate(this->wrate());
1346 if(laYer[j]>=ml && laYer[j]<mh) *p /= (DataType_t)v.
data[i];
1355 template<
class DataType_t>
1360 register DataType_t*
P=NULL;
1361 register DataType_t** pp;
1362 DataType_t
A, aL, aR;
1364 size_t nS, kS, nL, nR, lS;
1367 size_t nsub = t>0. ? size_t(this->
size()/this->wrate()/t+0.1) : 1;
1373 if((f>1. ||
bpp!=1.) && mode) {
1374 cout<<
"WSeries fraction(): invalid bpp: "<<
bpp<<
" fraction="<<f<<endl;
1379 size_t M = maxLayer()+1;
1382 pp = (DataType_t **)malloc(
sizeof(DataType_t*));
1393 lS = nS*nsub<S.
size() ? S.
size()-nS*nsub : 0;
1397 for(k=0; k<nsub; k++) {
1401 if(k+1 == nsub) nS += lS;
1403 nL = nS&1 ? nS/2 : nS/2-1;
1406 if(nL<1 || nR>nS-1) {
1407 cout<<
"WSeries::fraction() error: too short wavelet layer"<<endl;
1412 pp = (DataType_t **)realloc(pp,nS*
sizeof(DataType_t*));
1417 for(j=0; j<nS; j++) pp[j] = p + j*kS;
1419 this->waveSplit(pp,0,nS-1,nL);
1420 this->waveSplit(pp,nL,nS-1,nR);
1421 aL = *pp[nL]; aR = *pp[nR];
1423 for(j=0; j<nS; j++){
1426 if(j<nL) *P = (DataType_t)
fabs(A - aL);
1427 else if(j>nR) *P = (DataType_t)
fabs(A - aR);
1428 else { *P = 0; nZero++; }
1436 if(mode == 1)
continue;
1440 for(j=0; j<nS; j++){
1441 if(a.
data[j] == 0.)
continue;
1442 do{ r =
int(nS*drand48()-0.1);}
1443 while(p[r*kS] != 0);
1444 p[r*kS] = a.
data[
j];
1453 if(drand48() > f) { this->
data[
i] = 0; nZero++; }
1458 for(i=0; i<
M; i++) {
1459 if(this->
data[i]==0.) nZero++;
1464 return double(this->
size()-nZero)/double(this->
size());
1468 template<
class DataType_t>
1474 double tsRate = this->
rate();
1477 size_t M = maxLayer()+1;
1478 size_t il = size_t(2.*M*getlow()/tsRate);
1479 size_t ih = size_t(2.*M*gethigh()/tsRate+0.5);
1481 double ratio = double(this->
size());
1485 cout<<
"WSeries::significance(): invalid low and high: ";
1486 cout<<
"low = "<<il<<
" high = "<<ih<<endl;
1495 if(i>=il && i<=ih)
continue;
1503 for(j=0; j<
k; j++) p[j*m] = 0;
1505 ratio /= this->
size();
1511 size_t n = size_t(
fabs(T)*tsRate/S.
stride()/ratio+0.1);
1512 if(n<1) n = S.
size();
1521 nB = size_t(
bpp*nS*ratio);
1523 if(!nS || !nB || tsRate<=0. || m*S.
size()!=this->
size()) {
1524 cout<<
"WSeries::significance() error: invalid parameters"<<endl;
1528 l = (S.
size()-k*
n)*m;
1533 pp = (DataType_t **)malloc(nS*
sizeof(DataType_t*));
1543 for(j=0; j<nS; j++) {
1544 if(*p == 0.) {p++;
continue;}
1545 *p = (DataType_t)
fabs(
double(*p));
1550 if(nP>2) this->waveSort(pp,0,nP-1);
1552 for(j=0; j<nP; j++) {
1553 if(!i && l && pp[j]>=this->
data+l)
continue;
1554 *pp[
j] = nP<nB ? (DataType_t)
log(
double(nP)/(nP-
j)) :
1555 (DataType_t)
log(
double(nB)/(nP-
j));
1562 p = this->
data+i*nS+
l;
1567 return double(nZero)/ratio/this->
size();
1571 template<
class DataType_t>
1575 DataType_t*
px=NULL;
1590 size_t N = S.
size();
1598 nB = size_t(
bpp*nS);
1603 if(!nS || !nB || this->
rate()<=0. || m*S.
size()!=this->
size()) {
1604 cout<<
"WSeries::significance() error: invalid WSeries"<<endl;
1608 pp = (DataType_t **)malloc(nS*
sizeof(DataType_t*));
1609 xx = (DataType_t *)malloc(nS*
sizeof(DataType_t));
1610 qq = (DataType_t **)malloc(nS*
sizeof(DataType_t*));
1611 yy = (DataType_t *)malloc(nS*
sizeof(DataType_t));
1614 for(j=0; j<nS; j++){
1625 this->waveSplit(pp,0,nS-1,nL-1);
1626 this->waveSplit(pp,nL,nS-1,nR);
1627 aL = *pp[nL]; aR = *pp[nR];
1629 for(j=0; j<nL; j++) yy[j] = (DataType_t)
fabs(*pp[j] - aL);
1630 for(j=nR; j<nS; j++) yy[j+nL-nR] = (DataType_t)
fabs(*pp[j] - aR);
1632 this->waveSort(qq,0,nB-1);
1634 for(j=0; j<nB; j++) {
1636 if(index>nL) index+=nR-nL;
1638 if(next != index/m)
continue;
1639 this->
data[index-next*m+i*
m] = (DataType_t)
log(
double(nB)/(nB-
j));
1645 for(j=0; j<
m; j++) { *(px++) = *p; *p++ = 0;}
1650 if(next>2*n) next = 0;
1651 if(last>2*n) last = 0;
1660 return double(nBlack)/double(this->
size());
1664 template<
class DataType_t>
1677 double tsRate = this->
rate();
1685 size_t M = maxLayer()+1;
1686 size_t il = size_t(2.*M*getlow()/tsRate);
1687 size_t ih = size_t(2.*M*gethigh()/tsRate+0.5);
1691 cout<<
"WSeries::significance(): invalid low and high: ";
1692 cout<<
"low = "<<il<<
" high = "<<ih<<endl;
1699 for(j=0; j<il; j++){ this->
getLayer(wa,j); wa=1234567891.; this->putLayer(wa,j); }
1700 for(j=ih; j<
M; j++){ this->
getLayer(wa,j); wa=1234567891.; this->putLayer(wa,j); }
1706 size_t N = S.
size();
1707 size_t k = size_t(t*this->wrate());
1708 size_t n = size_t(T*this->wrate()/2.);
1710 if(t<=0. || k<1) k = 1;
1711 if(T<=0. || n<1) n = 1;
1714 size_t nW = (2*n+
k)*M;
1720 nB = size_t(
bpp*nS);
1725 if(!nS || !nB || this->
rate()<=0. || M*S.
size()!=this->
size()) {
1726 cout<<
"WSeries::significance() error: invalid WSeries"<<endl;
1730 pp = (DataType_t **)malloc(nS*
sizeof(DataType_t*));
1731 px = (DataType_t **)malloc(nS*
sizeof(DataType_t*));
1732 py = (DataType_t **)malloc(nS*
sizeof(DataType_t*));
1733 xx = (DataType_t *)malloc(nS*
sizeof(DataType_t));
1734 yy = (DataType_t *)malloc(nS*
sizeof(DataType_t));
1738 for(i=0; i<nW; i++){
1739 if(*p != 1234567891.){
1752 cout<<
"wseries::rSignificance() error 1 - illegal sample count"<<endl;
1756 for(i=0; i<
N; i+=
k){
1758 this->waveSplit(px,0,nS-1,nL-1);
1759 this->waveSplit(px,nL,nS-1,nR);
1760 aL = *px[nL]; aR = *px[nR];
1762 for(j=0; j<nL; j++) yy[j] = (DataType_t)
fabs(*px[j] - aL);
1763 for(j=nR; j<nS; j++) yy[j+nL-nR] = (DataType_t)
fabs(*px[j] - aR);
1765 if(nB != nS-nR+nL) {
1766 cout<<
"wseries::rSignificance: nB="<<nB<<
", N="<<nS-nR+nL<<endl;
1770 this->waveSort(py,0,nB-1);
1772 for(j=0; j<nB; j++) {
1774 if(index>nL) index+=nR-nL;
1778 if(index>=i*M && index<(i+k)*M) {
1779 if(this->data[index]!=0) {
1780 cout<<
"WSeries::rSignificance error: "<<this->data[
index]<<endl;
1782 this->data[
index] = DataType_t(
log(
double(nB)/(nB-j)));
1787 for(l=i; l<i+
k; l++){
1790 for(j=0; j<
M; j++) {
1791 if(*p != 1234567891.) { xx[J] = *p; pp[J++]=p; }
1796 cout<<
"wseries::rSignificance() error 2 - illegal sample count"<<endl;
1802 if(next==2*n+k) next = 0;
1803 if(last==2*n+k) last = 0;
1813 return double(nBlack)/double(this->
size());
1818 template<
class DataType_t>
1838 size_t M = maxLayer()+1;
1839 size_t il = size_t(2.*getlow()/this->wrate());
1840 size_t ih = size_t(2.*gethigh()/this->wrate()+0.5);
1844 cout<<
"WSeries::significance(): invalid low and high: ";
1845 cout<<
"low = "<<il<<
" high = "<<ih<<endl;
1852 for(j=0; j<il; j++){ this->
getLayer(wa,j); wa=1234567891.; this->putLayer(wa,j); }
1853 for(j=ih; j<
M; j++){ this->
getLayer(wa,j); wa=1234567891.; this->putLayer(wa,j); }
1859 size_t N = S.
size();
1860 size_t k = size_t(t*this->wrate());
1861 size_t n = size_t(T*this->wrate()/2.);
1863 if(t<=0. || k<1) k = 1;
1864 if(T<=0. || n<1) n = 1;
1867 size_t nW = (2*n+
k)*M;
1873 nB = size_t(
bpp*nS);
1878 if(!nS || !nB || this->
rate()<=0. || M*S.
size()!=this->
size()) {
1879 cout<<
"WSeries::gSignificance() error: invalid WSeries"<<endl;
1883 pp = (DataType_t **)malloc(nS*
sizeof(DataType_t*));
1884 px = (DataType_t **)malloc(nS*
sizeof(DataType_t*));
1885 py = (DataType_t **)malloc(nS*
sizeof(DataType_t*));
1886 xx = (DataType_t *)malloc(nS*
sizeof(DataType_t));
1887 yy = (DataType_t *)malloc(nS*
sizeof(DataType_t));
1891 for(i=0; i<nW; i++){
1892 if(*p != 1234567891.){
1905 cout<<
"wseries::gSignificance() error 1 - illegal sample count"<<endl;
1909 for(i=0; i<
N; i+=
k){
1911 this->waveSplit(px,0,nS-1,nL-1);
1912 this->waveSplit(px,nL,nS-1,nR);
1913 aL = *px[nL]; aR = *px[nR];
1915 for(j=0; j<nL; j++) yy[j] = (DataType_t)
fabs(*px[j]);
1916 for(j=nR; j<nS; j++) yy[j+nL-nR] = (DataType_t)
fabs(*px[j]);
1918 if(nB != nS-nR+nL) {
1919 cout<<
"wseries::gSignificance: nB="<<nB<<
", N="<<nS-nR+nL<<endl;
1923 this->waveSort(py,0,nB-1);
1925 for(j=0; j<nB; j++) {
1927 if(index>nL) index+=nR-nL;
1929 *(pp[
index]) = pow(*py[j]+1.11/2,2)/2./1.07 +
log(
bpp);
1933 for(l=i; l<i+
k; l++){
1936 for(j=0; j<
M; j++) {
1937 if(*p != 1234567891.) { xx[J] = *p; pp[J++]=p; }
1942 cout<<
"wseries::gSignificance() error 2 - illegal sample count"<<endl;
1948 if(next==2*n+k) next = 0;
1949 if(last==2*n+k) last = 0;
1959 return double(nBlack)/double(this->
size());
1964 template<
class DataType_t>
1983 size_t max_layer = maxLayer()+1;
1985 pc = ∾ pp = ≈ pm = p = NULL;
1990 for(k=1; k<=max_layer; k++){
1995 if(pp!=NULL) mp = pp->
size()/pc->
size();
1996 if(pm!=NULL) mm = pc->
size()/pm->
size();
2000 for(i=0; i<=
n; i++) {
2003 if(pc->
data[i] == 0.)
continue;
2004 if(pc->
data[i] > 9.7) cout<<
"pixclean: "<<pc->
data[
i]<<endl;
2006 if(i>0 && pc->
data[i-1]!=0.)
continue;
2007 if(i<n && pc->
data[i+1]!=0.)
continue;
2017 for(j=nm; j<
np; j++) {
2018 if(pp->
data[j] != 0) {
2034 for(j=nm; j<
np; j++) {
2035 if(pm->
data[j] != 0) {
2043 if(pc->
data[i]<S) {a.
data[
i]=0;
event--;}
2053 p = pm==NULL ? &am :
pm;
2059 return double(event)/this->
size();
2063 template<
class DataType_t>
2068 register DataType_t*
P=NULL;
2069 register DataType_t** pp;
2073 size_t nS, kS, mS, nL, nR;
2078 if((f>=1. ||
bpp!=1.) && mode) {
2079 cout<<
"WSeries percentile(): invalid bpp: "<<
bpp<<
" fraction="<<f<<endl;
2084 if(pin) *
this = *pin;
2086 size_t M = maxLayer()+1;
2090 size_t n0 = S.
size();
2091 if(n0) pp = (DataType_t **)malloc(n0*
sizeof(DataType_t*));
2105 nL = size_t(f*nS/2.+0.5);
2108 if(nL<2 || nR>nS-2) {
2109 cout<<
"WSeries::percentile() error: too short wavelet layer"<<endl;
2114 pp = (DataType_t **)realloc(pp,nS*
sizeof(DataType_t*));
2118 for(j=0; j<nS; j++) pp[j] = p + j*kS;
2120 this->waveSplit(pp,0,nS-1,nL-1);
2121 this->waveSplit(pp,nL,nS-1,nR);
2122 aL = double(*pp[nL-1]); aR = double(*pp[nR]);
2124 for(j=0; j<nS; j++){
2125 P = pp[
j]; A = double(*P);
2130 if(j<nL) *P = DataType_t(
fabs(A - aL));
2131 else if(j>nR) *P = DataType_t(
fabs(A - aR));
2132 else { *P = 0; nZero++; }
2134 if(mode == -1)
continue;
2135 if(pin) pin->
data[mS+P-p] = *
P;
2136 if(j>nL && j<nR)
continue;
2137 a.
data[(P-p)/kS] = *P;
2139 if(j>=nR) pp[nL+j-nR] =
P;
2142 if(mode == -1)
continue;
2145 this->waveSort(pp,0,nL-1);
2147 if(abs(mode)!=1) b =
a;
2149 for(j=0; j<nL; j++){
2152 x =
log(
double(nL)/(nL-j));
2153 *pp[
j] = mode==1 ? DataType_t(x) : 0;
2154 if(mode>1) a.
data[
r]=DataType_t(x);
2157 if(abs(mode)==1)
continue;
2161 for(j=0; j<nL; j++){
2163 do{ r =
int(nS*drand48()-0.1);}
2164 while(p[r*kS] != 0);
2165 p[r*kS] = a.
data[(P-p)/kS];
2166 if(pin) pin->
data[mS+r*kS] = b.
data[(P-p)/kS];
2174 if(drand48() > f) { this->
data[
i] = 0; nZero++; }
2180 if(this->
data[i]==0) nZero++;
2184 return double(this->
size()-nZero)/double(this->
size());
2188 template<
class DataType_t>
2197 size_t M = maxLayer()+1;
2199 double tsRate = this->
rate();
2201 double tleft, trigh;
2204 double cstep = 1./a.
rate();
2205 double sTart = this->
start();
2206 double sTop = this->
start()+this->
size()/tsRate;
2226 for(i=0; i<a.
size(); i++) {
2227 if(a.
start()+i/a.
rate() < sTart)
continue;
2236 for(i=0; i<g.
size(); i++) {
2237 if(g.
start()+i/g.
rate() < sTart)
continue;
2251 cout<<
"WSeries<DataType_t>::calibrate() no calibration data\n";
2264 righ += tsRate/2./S.
stride();
2265 if(righ > n*df)
break;
2270 while(left+count*df < righ){
2275 count++; pR++; pC++;
2283 for(i=0; i<alp.
size(); i++){
2285 if(alp.
data[i]<=0. || gam.
data[i]<=0.) {
2286 cout<<
"WSeries<DataType_t>::calibrate() zero alpha error\n";
2293 reH = 1.+(reH-1.)*gam.
data[
i];
2295 x.data[
i] = sqrt(reH*reH+imH*imH);
2309 sTart = alp.
start();
2310 sTop = alp.
start()+(alp.
size()-1)*cstep;
2313 trigh = sTart+cstep;
2318 if(t <= sTart) *p *= (DataType_t)
x.data[0];
2319 else if(t >= sTop) *p *= (DataType_t)
x.data[alp.
size()-1];
2321 if(t>trigh) { tleft=trigh; trigh+=cstep; m++; }
2322 c = (t-tleft)/cstep;
2323 *p *= DataType_t(
x.data[m]*(1-c) +
x.data[m+1]*c);
2333 template<
class DataType_t>
2341 template <
class DataType_t>
2347 if (R__b.IsReading()) {
2348 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
if (R__v) { }
2360 if(!bWWS) m_WaveType=-1;
2362 switch(m_WaveType) {
2385 this->setWavelet(*pW);
2398 R__b << pWavelet->m_WaveType;
2399 bool bWWS = ((pWavelet->pWWS==NULL)&&(pWavelet->m_WaveType==
HAAR)) ?
false :
true;
2401 pWavelet->Streamer(R__b);
2402 R__b.SetByteCount(R__c, kTRUE);
2408 #define CLASS_INSTANTIATION(class_) template class WSeries< class_ >;
2413 #undef CLASS_INSTANTIATION
wavearray< double > t(hp.size())
void mul(WSeries< DataType_t > &)
virtual void resize(unsigned int)
virtual size_t size() const
double gammaCL(double x, double n)
static double g(double e)
fprintf(stdout,"start=%f duration=%f rate=%f\n", x.start(), x.size()/x.rate(), x.rate())
virtual void rate(double r)
wavearray< double > a(hp.size())
virtual double percentile(double=0., int=0, WSeries< DataType_t > *=NULL)
param: f - black pixel fraction param: m - mode options: f = 0 - returns black pixel occupancy m = 1 ...
virtual void edge(double s)
virtual WSeries< float > variability(double=0., double=-1., double=-1.)
param: first - time window to calculate normalization constants second - low frequency boundary for c...
virtual WSeries< DataType_t > & operator+=(WSeries< DataType_t > &)
void bandpass(wavearray< DataType_t > &ts, double flow, double fhigh, int n=-1)
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\t layers : "<< nLAYERS<< "\t dF(hz) : "<< dF<< "\t dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1)*itime+ifreq;double time=itime *dT;double freq=(ifreq >0)?ifreq *dF:dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
wavearray< DataType_t > & operator=(const wavearray< DataType_t > &)
virtual double rSignificance(double, double=1., double=0.)
param: T - sliding window duration in seconds param: f - black pixel fraction param: t - sliding step...
std::slice getSlice(double n)
virtual wavearray< DataType_t > & operator*=(wavearray< DataType_t > &)
virtual double median(size_t=0, size_t=0) const
bool allocate(size_t, DataType_t *)
virtual void start(double s)
virtual WSeries< DataType_t > & operator*=(WSeries< DataType_t > &)
cout<< "Selected Pixels : "<< nPix<< endl;wc.cluster(1, 1);SSeries< double > ss
virtual double significance(double, double=1.)
param: n - sub-interval duration in seconds param: f - black pixel fraction options: f = 0 - returns ...
void setWavelet(const Wavelet &w)
size_t wavecount(double x, int n=0)
virtual double gSignificance(double, double=1., double=0.)
param: T - sliding window duration in seconds param: f - black pixel fraction param: t - sliding step...
virtual wavearray< double > filter(size_t)
param: n - number of decomposition steps algorithm: 1) do forward wavelet transform with n decomposit...
double maxEnergy(wavearray< DataType_t > &ts, Wavelet &w, double=0, int=1, int=0, TH1F *=NULL)
virtual WSeries< DataType_t > & operator-=(WSeries< DataType_t > &)
virtual double fraction(double=0., double=0., int=0)
param: t - sub interval duration. If can not divide on integer
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
TIter next(twave->GetListOfBranches())
#define CLASS_INSTANTIATION(class_)
virtual double coincidence(WSeries< DataType_t > &, int=0, int=0, double=0.)
param: WSeries object used for coincidence param: coincidence window in seconds return pixel occupanc...
virtual wavearray< DataType_t > & operator+=(wavearray< DataType_t > &)
virtual WSeries< DataType_t > & operator[](const std::slice &)
void wavescan(WSeries< DataType_t > **, int, TH1F *=NULL)
virtual Wavelet * Clone() const
return: Wavelet* - duplicate of *this, allocated on heap
virtual void lprFilter(wavearray< double > &)
virtual std::slice getSlice(const double)
wavearray< double > ts(N)
virtual void resample(double, int=6)
virtual wavearray< double > white(double, int=0, double=0., double=0.) const
virtual void stop(double s)
virtual wavearray< DataType_t > & operator-=(wavearray< DataType_t > &)
double fabs(const Complex &x)
void print()
param: int n if n<0, zero pixels defined in mask (regression) if n>=0, zero all pixels except ones de...
virtual void spesla(double, double, double=0.)
double Gamma2Gauss(TH1F *=NULL)
virtual double Coincidence(WSeries< DataType_t > &, double=0., double=0.)
param: WSeries object used for coincidence param: coincidence window in seconds param: threshold on s...
virtual void lprFilter(double, int=0, double=0., double=0.)
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
strcpy(RunLabel, RUN_LABEL)
Meyer< double > S(1024, 2)
virtual void median(double t, bool norm=false)
WSeries< DataType_t > & operator=(const wavearray< DataType_t > &)
virtual double pixclean(double=0.)
param: S - threshold on pixel significance return pixel occupancy.
WaveDWT< DataType_t > * pWavelet
double wdmPacket(int pattern, char opt='L', TH1F *=NULL)
virtual void resize(unsigned int)
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
virtual void waveSplit(DataType_t **pp, size_t l, size_t r, size_t m) const
void putLayer(wavearray< DataType_t > &, double n)
param: n - layer number
virtual void Dump(const char *, int=0)
virtual double rsignificance(size_t=0, double=1.)
param: n - sub-interval duration in domain units param: f - black pixel fraction options: f = 0 - ret...
virtual WSeries< double > calibrate(size_t, double, d_complex *, d_complex *, wavearray< double > &, wavearray< double > &, size_t ch=0)
param: number of samples in calibration arrays R & C param: frequency resolution param: pointer to re...
void cpf(const wavearray< DataType_t > &, int=0, int=0, int=0)
virtual WSeries< double > white(double, int, double=0., double=0.)
what it does: each wavelet layer is devided into k intervals.