16 #include <gnusstream.h>
24 #include "TObjArray.h"
25 #include "TObjString.h"
37 template<class DataType_t>
39 data(NULL),Size(0),
Rate(1.),Start(0.),Stop(0.),Edge(0.),fftw(NULL), ifftw(NULL)
45 template<
class DataType_t>
47 Rate(1.),Start(0.),Edge(0.),fftw(NULL),ifftw(NULL)
50 data = (DataType_t *)malloc(n*
sizeof(DataType_t));
58 template<
class DataType_t>
60 data(NULL),Size(0),
Rate(1.),Start(0.),Stop(0.),Edge(0.),fftw(NULL), ifftw(NULL)
64 template<
class DataType_t>
template<
class T>
67 data(NULL),Size(0),
Rate(1.),Start(0.),Edge(0.),fftw(NULL), ifftw(NULL)
70 if(n != 0 && p != NULL){
71 data = (DataType_t *)malloc(n*
sizeof(DataType_t));
72 for (i=0; i <
n; i++)
data[i] = p[i];
81 template<
class DataType_t>
90 template<
class DataType_t>
94 if (
this==&a)
return *
this;
105 for (i=0; i <
N; i++) { data[
i] = *
pm; pm+=
m; }
136 template<
class DataType_t>
141 unsigned int N = limit(a);
142 unsigned int n = Slice.stride();
147 for (i=Slice.start(); i<
N; i+=
n){ data[
i] = *p; p +=
m; }
154 template<
class DataType_t>
159 unsigned int n = Slice.stride();
160 unsigned int N = limit();
163 for (i=Slice.start(); i <
N; i+=
n) data[i] = c;
171 template<
class DataType_t>
176 unsigned int N = limit(a);
177 unsigned int n = Slice.stride();
182 for (i=Slice.start(); i<
N; i+=
n){ data[
i] += *p; p +=
m; }
189 template<
class DataType_t>
194 unsigned int n = Slice.stride();
195 unsigned int N = limit();
198 for (i=Slice.start(); i <
N; i+=
n) data[i] += c;
206 template<
class DataType_t>
211 unsigned int N = limit(a);
212 unsigned int n = Slice.stride();
217 for (i=Slice.start(); i<
N; i+=
n){ data[
i] -= *p; p +=
m; }
224 template<
class DataType_t>
229 unsigned int n = Slice.stride();
230 unsigned int N = limit();
233 for (i=Slice.start(); i <
N; i+=
n) data[i] -= c;
240 template<
class DataType_t>
245 unsigned int n = Slice.stride();
246 unsigned int N = limit();
249 for (i=Slice.start(); i <
N; i+=
n) data[i] *= c;
256 template<
class DataType_t>
261 unsigned int N = limit(a);
262 unsigned int n = Slice.stride();
267 for (i=Slice.start(); i<
N; i+=
n){ data[
i] *= *p; p +=
m; }
275 template<
class DataType_t>
280 if(limit() >
size()){
281 cout <<
"wavearray::operator[slice]: Illegal argument "<<limit()<<
" "<<
size()<<
"\n";
288 template<
class DataType_t>
293 cout <<
"wavearray::operator[int]: Illegal argument\n";
300 template<
class DataType_t>
305 if (app == 1)
strcpy (mode,
"a");
309 if ( (fp = fopen(fname, mode)) == NULL ) {
310 cout <<
" Dump() error: cannot open file " << fname <<
". \n";
315 fprintf( fp,
"# start time: -start %lf \n", this->Start );
316 fprintf( fp,
"# sampling rate: -rate %lf \n", this->
Rate );
317 fprintf( fp,
"# number of samples: -size %d \n", (
int)this->Size );
321 double dt = 1./this->
rate();
322 for (
int i = 0;
i <
n;
i++) {
324 fprintf( fp,
"%14f %e\n", time, data[
i]);
327 for (
int i = 0;
i <
n;
i++)
fprintf( fp,
"%e ", (
float)data[
i]);
334 template<
class DataType_t>
337 int n =
size() *
sizeof(DataType_t);
341 if (app == 1)
strcpy (mode,
"a");
345 if ( (fp=fopen(fname, mode)) == NULL ) {
346 cout <<
" DumpBinary() error : cannot open file " << fname <<
" \n";
350 fwrite(data, n, 1, fp);
355 template<
class DataType_t>
360 if (app == 1)
strcpy (mode,
"a");
363 if ( (fp = fopen(fname, mode)) == NULL ) {
364 cout <<
" DumpShort() error : cannot open file " << fname <<
". \n";
371 for (
int i=0;
i<
n;
i++ ) dtemp[
i]=
short(data[
i]) ;
373 n = n *
sizeof(short);
375 fwrite(dtemp, n, 1, fp);
381 template<
class DataType_t>
384 TFile* rfile =
new TFile(fname,
"RECREATE");
385 this->
Write(
"wavearray");
391 template<
class DataType_t>
394 int step =
sizeof(DataType_t);
398 if ( (fp=fopen(fname,
"r")) == NULL ) {
399 cout <<
" ReadBinary() error : cannot open file " << fname <<
". \n";
414 int n =
size() *
sizeof(DataType_t);
416 fread(data, n, 1, fp);
421 template<
class DataType_t>
425 dtmp =
new short[
size()];
426 int n =
size() *
sizeof(short);
429 if ( (fp=fopen(fname,
"r")) == NULL ) {
430 cout <<
" ReadShort() error : cannot open file " << fname <<
". \n";
434 cout <<
" Reading binary record, size="<< n <<
"\n";
436 fread(dtmp, n, 1, fp);
437 for (
unsigned int i = 0;
i <
size();
i++)
438 data[
i] = DataType_t(dtmp[
i]);
444 template<
class DataType_t>
447 DataType_t *p = NULL;
457 p = (DataType_t *)realloc(data,n*
sizeof(DataType_t));
462 Stop = Start + n/
Rate;
466 cout<<
"wavearray::resize(): memory allocation failed.\n";
483 template<
class DataType_t>
495 const DataType_t *p = a.
data;
500 double *temp=
new double[nF];
510 int nL =
int(nP2/ratio);
512 for (i = 0; i < nL; i++)
513 data[i] = (DataType_t)
Nevill(i*ratio, nP, p, temp);
517 int nM =
int((a.
size()-nP2)/ratio);
522 iL =
int(x) - nP2 + 1 ;
523 data[
i] = (DataType_t)
Nevill(x-iL, nP, p+iL, temp);
526 for (i = nL; i < nM; i+=2) {
528 iL =
int(x) - nP2 + 1 ;
529 data[
i] = (DataType_t)
Nevill(x-iL, nP, p+iL, temp);
531 iL =
int(x) - nP2 + 1 ;
532 data[i+1] = (DataType_t)
Nevill(x-iL, nP, p+iL, temp);
537 int nR = a.
size() - nP;
539 for (i = nM; i <
N; i++)
540 data[i] = (DataType_t)
Nevill(i*ratio-nR, nP, p, temp);
545 template<
class DataType_t>
553 template<
class DataType_t>
557 cout<<
"wavearray::Resample(): error - resample frequency must be >0\n";
561 double rsize = f/this->
rate()*this->
size();
563 if(fmod(rsize,2)!=0) {
564 cout<<
"wavearray::Resample(): error - resample frequency (" << f <<
") not allowed\n";
572 for(
int i=size;
i<rsize;
i++) this->data[
i]=0;
577 template<
class DataType_t>
583 int ibeg=0;
int iend=0;
585 if(this->data[
i]!=0 && ibeg==0) ibeg=
i;
586 if(this->data[
i]!=0) iend=
i;
588 int ilen=iend-ibeg+1;
591 int isize = 2*ilen+2*ishift;
592 isize = isize + (isize%4 ? 4 - isize%4 : 0);
594 w.rate(this->
rate()); w=0;
596 for(
int i=0;
i<ilen;
i++) {w[
i+ishift+ilen/2]=this->data[ibeg+
i];this->data[ibeg+
i]=0;}
602 double df = w.rate()/w.size();
604 for (
int ii=0;ii<(
int)w.size()/2;ii++) {
605 TComplex X(w[2*ii],w[2*ii+1]);
606 X=X*C.Exp(TComplex(0.,-2*pi*ii*df*T));
613 for(
int i=0;
i<(
int)w.size();
i++) {
614 int j=ibeg-(ishift+ilen/2)+
i;
615 if((j>=0)&&(j<(
int)this->
size())) this->data[
j]=w[
i];
620 template<
class DataType_t>
624 int i1,
N, n1, n2,
n, np2 = np/2;
632 N =
int(n/ratio + 0.5);
637 for (
int i = 0;
i <
np;
i++) {
639 for (
int j = 0;
j <
np;
j++)
if (
j !=
i) m *= (
i -
j);
643 for (
int i = 0;
i <
N;
i++)
648 x = x - i1 + np2 - 1;
652 n2 = i1 + np2 + 1 -
n;
662 for (
int j = 0;
j <
np;
j++) v[
j] = c[
j]*a.
data[i1 +
j - np2 + 1];
676 for (
int k = 0;
k <
np;
k++) {
677 for (
int j = 0;
j <
np;
j++)
if (
j !=
k) v[
j] *=
x;
681 for (
int j = 0;
j <
np;
j++) s += v[
j];
683 data[
i] = (DataType_t)s;
707 length = ((
size() - pos) < (a.
size() - a_pos))?
708 (
size() - pos) : (a.
size() - a_pos);
710 if( length > (
int)(
size() - pos) ) length =
size() - pos;
711 if( length > (
int)(a.
size() - a_pos) ) length = a.
size() - a_pos;
714 data[
i + pos] = a.
data[
i + a_pos];
734 length = ((
size() - pos) < (a.
size() - a_pos))?
735 (
size() - pos) : (a.
size() - a_pos);
737 if( length > (
int)(
size()- pos) ) length =
size() - pos;
738 if( length > (
int)(a.
size() - a_pos) ) length = a.
size() - a_pos;
741 data[
i + pos] += a.
data[
i + a_pos];
760 length = ((
size() - pos) < (a.
size() - a_pos))?
761 (
size() - pos) : (a.
size() - a_pos);
763 if( length > (
int)(
size() - pos) ) length =
size() - pos;
764 if( length > (
int)(a.
size() - a_pos) ) length = a.
size() - a_pos;
767 data[
i + pos] -= a.
data[
i + a_pos];
777 size_t n = this->
size();
781 cout <<
"wavearray::append() warning: sample rate mismatch.\n";
783 if(m == 0 )
return this->
size();
797 size_t n = this->
size();
800 this->Stop += 1./
rate();
813 template<
class DataType_t>
828 for (
int i=0;
i<
N;
i++) { a[
i] = data[
i]; b[
i]=0.;}
833 for (
int i=0;
i<n2;
i++)
834 { data[2*
i] = (DataType_t)a[
i]/N;
835 data[2*
i+1] = (DataType_t)b[
i]/N;
840 data[1] = (DataType_t)a[n2]/N;
844 if ((N&1) == 1) data[N-1] = (DataType_t)b[n2]/
N;
852 for (
int i=1;
i<n2;
i++)
863 { a[n2]=data[1]; b[n2]=data[N-1]; }
865 { a[n2]=data[1]; b[n2]=0.; }
869 for (
int i=0;
i<
N;
i++)
870 data[
i] = (DataType_t)a[
i];
877 template<
class DataType_t>
893 fftw =
new TFFTRealComplex(1,&N,
false);
894 fftw->Init(
"LS",sign,&kind);
899 for (
int i=0;
i<
N;
i++) { a[
i] = data[
i]; b[
i]=0.;}
903 fftw->GetPointsComplex(a, b);
906 for (
int i=0;
i<n2;
i++)
907 { data[2*
i] = (DataType_t)a[
i]/N;
908 data[2*
i+1] = (DataType_t)b[
i]/N;
913 data[1] = (DataType_t)a[n2]/N;
917 if ((N&1) == 1) data[N-1] = (DataType_t)b[n2]/
N;
925 ifftw =
new TFFTComplexReal(1,&N,
false);
926 ifftw->Init(
"LS",sign,&kind);
930 for (
int i=1;
i<n2;
i++)
941 { a[n2]=data[1]; b[n2]=data[N-1]; }
943 { a[n2]=data[1]; b[n2]=0.; }
945 ifftw->SetPointsComplex(a,b);
949 for (
int i=0;
i<
N;
i++)
950 data[
i] = (DataType_t)a[
i];
958 template<
class DataType_t>
961 if(fftw) {
delete fftw;fftw=NULL;}
962 if(ifftw) {
delete ifftw;ifftw=NULL;}
971 template<
class DataType_t>
981 cout <<
" Stack() error: data length too short to contain \n"
982 << length <<
" samples\n";
986 if (
size() != (
unsigned int)length)
resize(length);
993 data[
i] = (DataType_t)ss/
k;
1001 data[
i] -= (DataType_t)s0;
1002 s2 += data[
i]*data[
i];
1015 template<
class DataType_t>
1026 cout <<
" Stack() error: data length too short to contain \n"
1027 << length <<
" samples\n";
1033 *
this *= DataType_t(1./k);
1034 getStatistics(avr,rms);
1035 *
this -= DataType_t(avr);
1044 template<
class DataType_t>
1048 return this->Stack(td,
int(td.
rate() * window));
1052 template<
class DataType_t>
1056 register double x = 0;
1058 if(!
size())
return 0.;
1060 for(i=0; i<
size(); i++) x += data[i];
1064 template<
class DataType_t>
1070 if(!
size())
return 0;
1072 size_t nn = size_t(this->Edge*
rate());
1074 if(!nn || 2*nn>=
size()-2) {
return mean();}
1078 size_t N =
size()-2*nn;
1079 size_t mL = f<0. ? (N-
int(N*ff))/2 : 0;
1080 size_t mR = f<0. ? (N+
int(N*ff))/2-1 :
int(N*ff)-1;
1085 if(f>0) wtmp *= wtmp;
1086 DataType_t *p = wtmp.
data;
1087 DataType_t **pp = (DataType_t **)malloc(N*
sizeof(DataType_t*));
1088 for(i=nn; i<N+nn; i++) pp[i-nn] = p + i;
1090 if(mL>0) waveSplit(pp,0,N-1,mL);
1091 if(mR<N-2) waveSplit(pp,0,N-1,mR);
1092 for(i=mL; i<=mR; i++) x += this->data[pp[i]-p];
1095 return x/(mR-mL+1.);
1098 template<
class DataType_t>
1102 register double x = 0.;
1103 register DataType_t *p = data + s.
start();
1104 size_t N = s.
size();
1106 if(
size()<limit(s)) N = (limit(s) - s.
start() - 1)/m;
1107 for(i=0; i<
N; i++) { x += *p; p +=
m; }
1108 return (N==0) ? 0. : x/
N;
1112 template<
class DataType_t>
1129 size_t step = Slice.stride();
1130 size_t N = Slice.size();
1131 size_t n = size_t(t*
rate()/step);
1134 cout<<
"wavearray<DataType_t>::mean() short time window"<<endl;
1149 xx = (DataType_t *)malloc((n+1)*
sizeof(DataType_t));
1151 p = data+Slice.start();
1152 q = data+Slice.start();
1153 for(i=0; i<=
n; i++) {
1163 pm->
data[i/skip] = DataType_t(sum/(n+1.));
1164 if(clean) q[i*step] -= DataType_t(sum/(n+1.));
1167 if(clean) q[i*step] -= DataType_t(sum/(n+1.));
1168 else q[i*step] = DataType_t(sum/(n+1.));
1178 if(last>n) last = 0;
1187 template<
class DataType_t>
1191 register double x = 0.;
1192 register double y = 0.;
1193 size_t N = (
size()>>2)<<2;
1194 register DataType_t *p = data +
size() -
N;
1196 if(!
size())
return 0.;
1198 for(i=0; i<
size()-
N; i++) { x += data[
i]; y += data[
i]*data[
i]; }
1199 for(i=0; i<
N; i+=4){
1200 x += p[
i] + p[i+1] + p[i+2] + p[i+3];
1201 y += p[
i]*p[
i] + p[i+1]*p[i+1] + p[i+2]*p[i+2] + p[i+3]*p[i+3];
1204 return sqrt(y/
size()-x*x);
1208 template<
class DataType_t>
1222 size_t step = Slice.stride();
1223 size_t N = Slice.size();
1224 size_t n = size_t(t*
rate()/step);
1227 cout<<
"wavearray<DataType_t>::median() short time window"<<endl;
1242 pp = (DataType_t **)malloc((n+1)*
sizeof(DataType_t*));
1243 xx = (DataType_t *)malloc((n+1)*
sizeof(DataType_t));
1245 p = data+Slice.start();
1246 q = data+Slice.start();
1247 for(i=0; i<=
n; i++) {
1248 xx[
i] = *p>0 ? *p : -(*p);
1256 if(i==(i/skip)*skip) {
1257 waveSplit(pp,0,n,nM);
1262 pm->
data[i/skip] = DataType_t(rm/0.6745);
1263 if(clean) q[i*step] *= DataType_t(0.6745/rm);
1266 if(clean) q[i*step] *= DataType_t(0.6745/rm);
1267 else q[i*step] = DataType_t(rm/0.6745);
1271 xx[last++] = *p>0 ? *p : -(*p);
1275 if(last>n) last = 0;
1286 template<
class DataType_t>
1290 register double a = 0.;
1291 register double x = 0.;
1292 register double y = 0.;
1293 register DataType_t *p = data + s.
start();
1294 size_t n = s.
size();
1297 if(
size()<limit(s)) n = (limit(s) - s.
start() - 1)/m;
1299 size_t N = (n>>2)<<2;
1301 for(i=0; i<n-
N; i++) {
1302 a = *p; x +=
a; y += a*
a; p +=
m;
1305 for(i=0; i<
N; i+=4) {
1306 a = *p; x +=
a; y += a*
a; p +=
m;
1307 a = *p; x +=
a; y += a*
a; p +=
m;
1308 a = *p; x +=
a; y += a*
a; p +=
m;
1309 a = *p; x +=
a; y += a*
a; p +=
m;
1312 return sqrt(y/N-x*x);
1315 template<
class DataType_t>
1318 if(!
size())
return 0;
1319 register unsigned int i;
1320 register DataType_t
x = data[0];
1321 for(i=1; i<
size(); i++) {
if(x<data[i]) x=data[
i]; }
1325 template<
class DataType_t>
1332 cout<<
"wavearray::max(): illegal imput array\n";
1337 size_t K = Slice.stride();
1338 size_t I = Slice.start();
1339 size_t N = Slice.size();
1341 DataType_t* pin = x.
data+I;
1342 DataType_t* pou = this->data+I;
1343 for(n=0; n<
N; n+=
K) {
1344 if(pou[n]<pin[n]) pou[
n]=pin[
n];
1349 template<
class DataType_t>
1352 if(!
size())
return 0;
1354 register DataType_t
x = data[0];
1355 for(i=1; i<
size(); i++) {
if(x>data[i]) x=data[
i]; }
1360 template<
class DataType_t>
1364 register int i = l-1;
1368 data[
n]=data[
r]; data[
r]=
v;
1372 while(data[++i]<v && i<j);
1373 while(data[--j]>v && i<j);
1375 data[
r]=data[
n]; data[
n]=
v;
1380 template<
class DataType_t>
1384 register int i = l-1;
1388 data[
n]=data[
r]; data[
r]=
v;
1393 while((data[++i]>0 ? data[i] : -data[i])<vv && i<j);
1394 while((data[--j]>0 ? data[j] : -data[j])>vv && i<j);
1396 data[
r]=data[
n]; data[
n]=
v;
1402 template<
class DataType_t>
1411 register DataType_t
v;
1412 register DataType_t* p;
1413 register DataType_t** pp = pin;
1415 if(pp==NULL || !this->
size())
return;
1416 if(!r) r = this->
size()-1;
1419 register size_t i = (r+
l)/2;
1420 register size_t j = r-1;
1423 if(*pp[l] > *pp[i]) {p=pp[
l]; pp[
l]=pp[
i]; pp[
i]=p;}
1424 if(*pp[l] > *pp[r]) {p=pp[
l]; pp[
l]=pp[
r]; pp[
r]=p;}
1425 if(*pp[i] > *pp[r]) {p=pp[
i]; pp[
i]=pp[
r]; pp[
r]=p;}
1429 p=pp[
i]; pp[
i]=pp[
j]; pp[
j]=p;
1434 while(*pp[++i] < v);
1435 while(*pp[--j] > v);
1437 p=pp[
i]; pp[
i]=pp[
j]; pp[
j]=p;
1440 p=pp[
i]; pp[i++]=pp[r-1]; pp[r-1]=p;
1442 if(j-l > 2) waveSort(pp,l,j);
1445 if(*pp[l] > *pp[k]) {p=pp[
l]; pp[
l]=pp[
k]; pp[
k]=p;}
1446 if(*pp[l] > *pp[j]) {p=pp[
l]; pp[
l]=pp[
j]; pp[
j]=p;}
1447 if(*pp[k] > *pp[j]) {p=pp[
k]; pp[
k]=pp[
j]; pp[
j]=p;}
1450 if(r-i > 2) waveSort(pp,i,r);
1453 if(*pp[i] > *pp[k]) {p=pp[
i]; pp[
i]=pp[
k]; pp[
k]=p;}
1454 if(*pp[i] > *pp[r]) {p=pp[
i]; pp[
i]=pp[
r]; pp[
r]=p;}
1455 if(*pp[k] > *pp[r]) {p=pp[
k]; pp[
k]=pp[
r]; pp[
r]=p;}
1461 template<
class DataType_t>
1466 size_t N = this->
size();
1467 DataType_t *pd = (DataType_t *)malloc(N*
sizeof(DataType_t));
1468 DataType_t **pp = (DataType_t **)malloc(N*
sizeof(DataType_t*));
1469 for(
size_t i=0;
i<
N;
i++) {
1470 pd[
i] = this->data[
i];
1474 for(
size_t i=0;
i<
N;
i++) this->data[
i] = *pp[
i];
1480 template<
class DataType_t>
1492 if(*pp[l] > *pp[i]) {p=pp[
l]; pp[
l]=pp[
i]; pp[
i]=p;}
1493 if(*pp[l] > *pp[r]) {p=pp[
l]; pp[
l]=pp[
r]; pp[
r]=p;}
1494 if(*pp[i] > *pp[r]) {p=pp[
i]; pp[
i]=pp[
r]; pp[
r]=p;}
1498 p=pp[
i]; pp[
i]=pp[
j]; pp[
j]=p;
1503 while(*pp[++i] < v);
1504 while(*pp[--j] > v);
1506 p=pp[
i]; pp[
i]=pp[
j]; pp[
j]=p;
1508 p=pp[
i]; pp[
i]=pp[r-1]; pp[r-1]=p;
1510 if(i > m) waveSplit(pp,l,i,m);
1511 else if(i < m) waveSplit(pp,i,r,m);
1516 template<
class DataType_t>
1522 size_t N = this->
size();
1523 DataType_t *pd = (DataType_t *)malloc(N*
sizeof(DataType_t));
1524 DataType_t **pp = (DataType_t **)malloc(N*
sizeof(DataType_t*));
1525 for(
size_t i=0;
i<
N;
i++) {
1526 pd[
i] = this->data[
i];
1529 waveSplit(pp,l,r,m);
1530 for(
size_t i=0;
i<
N;
i++) this->data[
i] = *pp[
i];
1534 return this->data[
m];
1537 template<
class DataType_t>
1543 size_t N = this->
size();
1544 DataType_t *pd = (DataType_t *)malloc(N*
sizeof(DataType_t));
1545 DataType_t **pp = (DataType_t **)malloc(N*
sizeof(DataType_t*));
1546 for(
size_t i=0;
i<
N;
i++) {
1547 pd[
i] = this->data[
i];
1550 waveSplit(pp,0,N-1,m);
1551 for(
size_t i=0;
i<
N;
i++) this->data[
i] = *pp[
i];
1557 template<
class DataType_t>
1560 if(!r) r =
size()-1;
1565 size_t m = N/2+(N&1);
1569 DataType_t **pp = (DataType_t **)malloc(N*
sizeof(DataType_t*));
1570 for(i=l; i<=
r; i++) pp[i-l] = data + i;
1572 waveSplit(pp,0,N-1,m);
1580 template<
class DataType_t>
1598 size_t step = Slice.stride();
1599 size_t N = Slice.size();
1600 size_t n = size_t(t*
rate()/step);
1603 cout<<
"wavearray<DataType_t>::median() short time window"<<endl;
1618 pp = (DataType_t **)malloc((n+1)*
sizeof(DataType_t*));
1619 xx = (DataType_t *)malloc((n+1)*
sizeof(DataType_t));
1621 p = data+Slice.start();
1622 q = data+Slice.start();
1623 for(i=0; i<=
n; i++) {
1632 if(i==(i/skip)*skip) {
1633 waveSplit(pp,0,n,nM);
1638 pm->
data[i/skip] = am;
1639 if(clean) q[i*step] -= am;
1642 if(clean) q[i*step] -= am;
1643 else q[i*step] = am;
1651 if(last>n) last = 0;
1661 template<
class DataType_t>
1671 size_t N = Slice.size();
1672 size_t step = Slice.stride();
1673 size_t n = size_t(t*
rate()/step);
1676 cout<<
"wavearray<DataType_t>::median() short time window"<<endl;
1685 DataType_t** pp = (DataType_t **)malloc((n+1)*
sizeof(DataType_t*));
1688 p = data+Slice.start();
1689 q = data+Slice.start();
1690 for(i=0; i<=
n; i++) {
1701 q[i*step] = DataType_t(r>0. ? -
log(1.-r) :
log(1.+r));
1704 xx.
data[last++] = *p; p+=step;
1708 if(next>n) next = 0;
1709 if(last>n) last = 0;
1718 template<
class DataType_t>
1724 DataType_t **pp = NULL;
1729 pp = (DataType_t **)malloc(n*
sizeof(DataType_t*));
1733 for(i=0; i<
n; i++) pp[i] = data + i;
1734 qsort(pp, n,
sizeof(DataType_t *),
1738 if(i==0) out = *(pp[0]);
1739 else if(i>=n-1) out = *(pp[n-1]);
1740 else out = (*(pp[
i])+*(pp[i+1]))/2;
1742 for(i=0; i<
n; i++) *(pp[i]) = n-
i;
1751 template<
class DataType_t>
1773 x1 *= DataType_t(0.5);
1777 for(k=1; k<
K; k++) {
1780 for(j=offset; j<offset+
M; j++) {
1785 for(j=1; j<
N; j++) {
1789 p += (xp-xm)*(xp+xm);
1792 q += (xp-xm)*(xp+xm);
1794 xp = k==1 ? 2*p/q : p/q;
1807 template<
class DataType_t>
1816 for(i=0; i<
N; i++) {
1817 for(j=1; j<m && (i-
j)>=0; j++) {
1818 data[
i] += DataType_t(w.
data[j]*x.
data[i-j]);
1827 template<
class DataType_t>
1837 int k = stride>0 ?
int(duration/stride) : 1;
1840 int m =
int((N-2.*oFFset*
rate())/k);
1844 if(offset&1) offset++;
1851 for(l=0; l<
k; l++) {
1858 nl = l==k-1 ? N : n+
m;
1861 if(mode == 0 || mode == 1) {
1862 for(i=nf; i<nl; i++) {
1864 b = (!mode && i<N-
n) ? 0.5 : 1.;
1865 for(j=1; j<
M; j++) {
1867 data[
i] += DataType_t(a*b);
1872 if(mode == 0 || mode == -1) {
1873 for(i=nf; i<nl; i++) {
1874 if(i >= N-n)
continue;
1875 b = (!mode && i>=
n) ? 0.5 : 1.;
1876 for(j=1; j<
M; j++) {
1878 data[
i] += DataType_t(a*b);
1892 template<
class DataType_t>
1898 size_t nL = size_t(f*N+0.5);
1903 if(
size()<=offset) {
1904 cout<<
"wavearray<DataType_t>::getLPRFilter() invalid input parameters\n";
1915 DataType_t ** pp = (DataType_t **)malloc(N*
sizeof(DataType_t*));
1917 for(i=offset; i<
n; i++) pp[i-offset] = x.
data + i;
1919 waveSplit(pp,0,N-1,nn);
1920 waveSplit(pp,0,nn,nL);
1921 waveSplit(pp,nn,N-1,nR);
1924 for(i=0; i<nL; i++) *pp[i] = 0;
1925 for(i=nR; i<
N; i++) *pp[i] = 0;
1932 for (m=0; m<
M; m++) {
1934 for (i=offset; i<
n; i++) {
1943 double s = r.
data[0];
1945 a = 0.; a.
data[0] = 1.;
1947 for (m=1; m<
M; m++) {
1950 for (i=0; i<
m; i++) q -= a.
data[i] * r.
data[m-i];
2000 template<
class DataType_t>
2008 if(t<=0.) t = segT-2.*oFFset;
2010 if(offset&1) offset--;
2012 if(stride > t || stride<=0.) stride =
t;
2013 int K =
int((segT-2.*oFFset)/stride);
2022 int mL =
int(0.15865*m+0.5);
2025 int jR = N-offset-
m;
2041 if(m<3 || mL<2 || mR>m-2) {
2042 cout<<
"wavearray::white(): too short input array."<<endl;
2043 return mode!=0 ? norm50 : meDIan;
2046 register DataType_t *p = data;
2051 DataType_t ** pp = (DataType_t **)malloc(m*
sizeof(DataType_t*));
2053 for(j=0; j<=
K; j++) {
2055 if(jj < offset) p = data +
offset;
2056 else if(jj >= jR) p = data + jR;
2060 if(p+m>data+N) {cout<<
"wavearray::white(): error1\n";
exit(1);}
2062 for(i=0; i<
m; i++) pp[i] = p + i;
2063 waveSplit(pp,0,m-1,mm);
2064 waveSplit(pp,0,mm,mL);
2065 waveSplit(pp,mm,m-1,mR);
2066 meDIan[
j] = mode ? *pp[mm] : sqrt((*pp[mm])*0.7191);
2067 norm50[
j] = (*pp[mR] - *pp[mL])/2.;
2076 for(i=0; i<mm; i++){
2077 x = double(*p)-meDIan.
data[0];
2079 *(p++) = mode==1 ? DataType_t(x/r) : DataType_t(x/r/r);
2082 for(j=0; j<
K; j++) {
2083 for(i=0; i<
k; i++) {
2084 x = double(*p)-(meDIan.
data[j+1]*i + meDIan.
data[
j]*(k-
i))/
k;
2085 r = (norm50.
data[j+1]*i + norm50.
data[
j]*(k-
i))/
k;
2086 *(p++) = mode==1 ? DataType_t(x/r) : DataType_t(x/r/r);
2091 for(i=0; i<mm; i++){
2092 x = double(*p)-meDIan.
data[
K];
2094 *(p++) = mode==1 ? DataType_t(x/r) : DataType_t(x/r/r);
2100 return mode!=0 ? norm50 : meDIan;
2104 template<
class DataType_t>
2110 DataType_t *p =
const_cast<DataType_t *
>(data);
2112 size_t N =
size() - 1 + size_t(
size()&1);
2114 if(!
size())
return 0.;
2124 for(i=1; i<
N; i+=2) {
2125 a = p[
i]; b = p[i+1];
2133 r = r/double(N) - m*
m;
2135 y = 4.*(y/N - m*m + m*(p[0]+p[
i]-
m)/N);
2136 y /= 4.*r - 2.*((p[0]-
m)*(p[0]-m)+(p[
i]-
m)*(p[i]-m))/N;
2139 a = (
fabs(y) < 1) ? sqrt(0.5*(1.-
fabs(y))) : 0.;
2141 return y>0 ? -a :
a;
2144 template <
class DataType_t>
2149 if (R__b.IsReading()) {
2150 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
if (R__v) { }
2151 TNamed::Streamer(R__b);
2156 if(R__v>1) R__b >> Edge;
2157 R__b.StreamObject(&(Slice),
typeid(
slice));
2159 R__b.ReadFastArray((
char*)data,Size*
sizeof(DataType_t));
2163 TNamed::Streamer(R__b);
2169 R__b.StreamObject(&(Slice),
typeid(
slice));
2170 R__b.WriteFastArray((
char*)data,Size*
sizeof(DataType_t));
2171 R__b.SetByteCount(R__c, kTRUE);
2176 template <
class DataType_t>
2180 name.ReplaceAll(
" ",
"");
2182 TObjString* ext_tok = (TObjString*)token->At(token->GetEntries()-1);
2196 cout <<
"wavearray operator (>>) : file type " << ext.Data() <<
" not supported" << endl;
2202 template <
class DataType_t>
2205 cout << endl << endl;
2207 cout <<
"Size\t= " << this->Size <<
" samples" << endl;
2208 cout <<
"Rate\t= " << this->Rate <<
" Hz" << endl;
2209 cout <<
"Start\t= " << this->Start <<
" sec" << endl;
2210 cout <<
"Stop\t= " << this->Stop <<
" sec" << endl;
2211 cout <<
"Length\t= "<< this->Stop-this->Start <<
" sec" << endl;
2212 cout <<
"Edge\t= " << this->Edge <<
" sec" << endl;
2213 cout <<
"Mean\t= " << this->mean() << endl;
2214 cout <<
"RMS\t= " << this->rms() << endl;
2222 template<
class DataType_t>
2226 Interval ti = ts.getTStep();
2227 double Tsample = ti.GetSecs();
2229 unsigned int n=ts.getNSample();
2233 rate(
double(
int(1./Tsample + 0.5)));
2235 cout <<
" Invalid sampling interval = 0 sec.\n";
2238 start(
double(ts.getStartTime().totalS()));
2240 TSeries
r(ts.getStartTime(), ti, ts.getNSample());
2243 vec_ref= (
float*)(r.refData());
2245 for (
unsigned int i=0;
i<
n;
i++ )
2246 data[
i] = (DataType_t) (vec_ref[
i]);
2254 #define CLASS_INSTANTIATION(class_) template class wavearray< class_ >;
2266 #undef CLASS_INSTANTIATION
2269 #if !defined (__SUNPRO_CC)
2271 wavearray(
const double *,
unsigned int,
double);
2274 wavearray(
const float *,
unsigned int,
double );
2277 wavearray(
const short *,
unsigned int,
double );
2280 wavearray(
const double *,
unsigned int,
double );
2283 wavearray(
const float *,
unsigned int,
double );
2286 wavearray(
const short *,
unsigned int,
double );
2291 static void fake_instatination_SUN_CC ()
wavearray< double > t(hp.size())
size_t append(const wavearray< DataType_t > &)
virtual size_t size() const
fprintf(stdout,"start=%f duration=%f rate=%f\n", x.start(), x.size()/x.rate(), x.rate())
virtual void rate(double r)
virtual void ReadShort(const char *)
wavearray< double > a(hp.size())
double getStatistics(double &mean, double &rms) const
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 int getSampleRankE(size_t n, size_t l, size_t r) const
double median(std::vector< double > &vec)
void add(const wavearray< DataType_t > &, int=0, int=0, int=0)
virtual void ReadBinary(const char *, int=0)
virtual wavearray< DataType_t > & operator*=(wavearray< DataType_t > &)
double Nevill(const double x0, int n, DataType_t *p, double *q)
virtual double median(size_t=0, size_t=0) const
virtual void start(double s)
virtual DataType_t min() const
virtual double mean() const
cout<< "Selected Pixels : "<< nPix<< endl;wc.cluster(1, 1);SSeries< double > ss
DataType_t rank(double=0.5) const
virtual DataType_t max() const
virtual wavearray< double > getLPRFilter(size_t, size_t=0)
virtual void DumpObject(const char *)
void wavefft(double a[], double b[], int ntot, int n, int nspan, int isn)
void sub(const wavearray< DataType_t > &, int=0, int=0, int=0)
virtual char * operator>>(char *)
void Resample(const wavearray< DataType_t > &, double, int=6)
TIter next(twave->GetListOfBranches())
void Dump(WSeries< double > &w, double t1, double t2, const char *fname)
virtual wavearray< DataType_t > & operator+=(wavearray< DataType_t > &)
virtual int getSampleRank(size_t n, size_t l, size_t r) const
virtual wavearray< DataType_t > & operator<<(wavearray< DataType_t > &)
virtual void DumpBinary(const char *, int=0)
virtual void lprFilter(wavearray< double > &)
virtual void waveSort(DataType_t **pp, size_t l=0, size_t r=0) const
wavearray< double > ts(N)
virtual void delay(double T)
virtual void DumpShort(const char *, int=0)
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)
virtual void spesla(double, double, double=0.)
void resample(const wavearray< DataType_t > &, double, int=6)
virtual wavearray< DataType_t > & operator[](const std::slice &)
strcpy(RunLabel, RUN_LABEL)
virtual void Dump(const char *, int=0)
virtual void exponential(double)
int(* qsort_func)(const void *, const void *)
#define CLASS_INSTANTIATION(class_)
double Stack(const wavearray< DataType_t > &, int)
for(int i=0;i< 101;++i) Cos2[2][i]=0
virtual void resize(unsigned int)
virtual void waveSplit(DataType_t **pp, size_t l, size_t r, size_t m) const
void cpf(const wavearray< DataType_t > &, int=0, int=0, int=0)