Logo coherent WaveBurst  
Library Reference Guide
Logo
CompareWDMvsTime.C
Go to the documentation of this file.
1 // Read waveform from file
2 // resample data
3 // WDM transfom
4 // compare energy in time domain vs energy in time-freq domain
5 // plot WDM
6 // plot waveform in time/fft domain
7 // save plots to disk
8 // compare Time/WDM energies for each time bin
9 
10 #define WF_HP "Wolfy_hplus.txt" // Wolfy waveform
11 #define NZEROES 50000
12 
13 #define USE_WDM_00_90 // use 00 + 90 components
14  // if commented then only 00 component is used
15 
16 //#define SAVE_PLOTS // save plots to disk
17 
23 
25 
26  TString hp = WF_HP;
27 
28  // ---------------------------------------------------
29  // read waveform and fill x wavearray
30  // ---------------------------------------------------
32  wavearray<double> time;
33 
34  double X,T;
35  ifstream in;
36  in.open(hp.Data());
37  for(int i=0;i<NZEROES;i++) x.append(0); // padding with zeroes
38  while(1) {
39  in >> X >> T;
40  if(!in.good()) break;
41  time.append(X);
42  x.append(T);
43  }
44  for(int i=0;i<NZEROES;i++) x.append(0); // padding with zeroes
45 
46  double rate = (time.size()-1)/(time[time.size()-1]-time[0]);
47  x.rate(rate);
48  cout << "rate time series " << x.rate() << endl;
49  cout << "size time series " << x.size() << endl;
50  cout << "start time series " << x.start()<< endl;
51  cout << "stop time series " << x.stop() << endl;
52 
53  x.resample(8192); // resample data to 8192 Hz
54 
55  // ---------------------------------------------------
56  // compute energy in time domain
57  // ---------------------------------------------------
58  double E_TIME=0;
59  double dt=1./x.rate();
60  geT = new gwavearray<double>(x); // energy array
61  for(int i=0;i<x.size();i++) {geT->data[i]=x[i]*x[i]*dt;E_TIME+=geT->data[i];}
62  cout << "E_TIME = " << E_TIME << endl;
63 
64  // ---------------------------------------------------
65  // perform WDM transform level=512
66  // ---------------------------------------------------
67  double R = x.rate();
68  int level = 32;
69 
70  double dF = R/(2*level); // frequency resolution
71  double dT = level/R; // time resolution
72 
73  cout<< "R="<< R << " dF= " << dF << " dT= " << dT << endl;
74 
75  WDM<double> wdtf(level, level, 4, 10);
77  w.Forward(x, wdtf); // WDM Transform
78 
79  cout << "size WDM " << w.size() << endl;
80  cout << "start WDM " << w.start() << endl;
81  cout << "stop WDM " << w.stop() << endl;
82 
83  // ---------------------------------------------------
84  // compute energy in time-freq domain
85  // ---------------------------------------------------
86  int slices = w.sizeZero();
87  float df2 = w.resolution();
88  float dt2 = 1./(2*df2);
89  int layers = w.maxLayer()+1;
90  double E_WDM=0;
91  geWDM = new gwavearray<double>(slices); // energy array
92  geWDM->rate(w.rate()/level);geWDM->start(w.start());
93  for(int i=1;i<=slices;i++) {
94  geWDM->data[i-1] = 0;
95  for(int j=1;j<=layers;j++) { // loop over the frequencies with the same time
96 #ifdef USE_WDM_00_90
97  // en = (E00+E90)/2
98  double en = pow(w.getSample(i-1,j-1),2)+pow(w.getSample(i-1,-(j-1+0.01)),2);
99  en/=2*R;
100 #else
101  // en = E00
102  double en = pow(w.getSample(i-1,j-1),2)/R;
103 #endif
104  geWDM->data[i-1] += en;
105  E_WDM += en;
106  }
107  }
108  cout << "slices "<< slices << " layers"<< layers << endl;
109  cout << "E_WDM = " << E_WDM << endl;
110  cout << "E_WDM/E_TIME = " << E_WDM/E_TIME << endl;
111 
112  // ---------------------------------------------------
113  // Plot Waveform in time/fft domain
114  // ---------------------------------------------------
115  gx = new gwavearray<double>(&x);
116  gx->Draw(GWAT_TIME);
117  //gx->Draw(GWAT_TIME,"FULL");
118  //gx->Draw(GWAT_FFT);
119 
120 #ifdef SAVE_PLOTS
121  // save plot to file
122  watplot* plot = gx->GetWATPLOT(); // get pointer to watplot object
123  char gtitle[256]; sprintf(gtitle,"Wolfy waveform");
124  plot->gtitle(gtitle,"time(sec)","amplitude");
125  //plot->gtitle(gtitle,"frequency(hz)","amplitude");
126  TString gfile="Wolfy_Time.png";
127  (*plot) >> gfile;
128 #endif
129 
130  // ---------------------------------------------------
131  // Plot Waveform in time-freq domain
132  // Note : Energy in the plot is not normalize : (Norm is 1/R)
133  // ---------------------------------------------------
134  WTS = new watplot(const_cast<char*>("wtswrc"));
135  double start,stop;
136  gx->GetTimeRange(start,stop);
137  double flow = 64;
138  double fhigh = 4096;
139 #ifdef USE_WDM_00_90
140  WTS->plot(&w, 2, start, stop,const_cast<char*>("COLZ"));
141 #else
142  WTS->plot(&w, 4, start, stop,const_cast<char*>("COLZ"));
143 #endif
144  WTS->hist2D->GetYaxis()->SetRangeUser(flow, fhigh);
145 #ifdef SAVE_PLOTS
146  WTS->canvas->SaveAs("Wolfy_WDM.png"); // dump plot
147 #endif
148 
149 
150  // ---------------------------------------------------
151  // Compare energy eT vs eWDM
152  // ---------------------------------------------------
153  greT = new gwavearray<double>(geWDM->size()); // energy array
154  greT->start(geWDM->start());
155  greT->rate(geWDM->rate());
156  // energy of time domain signal is integrated over the WDM bin size (=level)
157  for(int i=0;i<greT->size()-1;i++) {
158  greT->data[i]=0;
159  for(int j=0;j<level;j++) greT->data[i]+=geT->data[i*level+j];
160  }
161 
162 // geT->Draw(GWAT_TIME);
163 // geWDM->Draw(GWAT_TIME);
164  greT->Draw(GWAT_TIME); // black is greT
165  // the left edge of the bin is shifted in time by +half_sample because
166  // the WDM pixel time is the central time of the pixel (the same is for the frequency)
167  geWDM->start(geWDM->start()-1./(2.*geWDM->rate()));
168  greT->Draw(geWDM,GWAT_TIME,"SAME",kRed); // red is geWDM
169 #ifdef SAVE_PLOTS
170  // save plot to file
171  plot = greT->GetWATPLOT(); // get pointer to watplot object
172  gfile="Wolfy_Time_Freq_Comp.png";
173  (*plot) >> gfile;
174 #endif
175 }
size_t append(const wavearray< DataType_t > &)
Definition: wavearray.cc:793
void gtitle(TString title="", TString xtitle="", TString ytitle="")
Definition: watplot.cc:1256
int slices
virtual void rate(double r)
Definition: wavearray.hh:141
TString("c")
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<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ 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
gwavearray< double > * geT
plot gtitle(gtitle,"frequency (Hz)","strain/#sqrt{Hz}")
int layers
wavearray< double > hp
Definition: DrawInspiral.C:43
virtual void start(double s)
Definition: wavearray.hh:137
int j
Definition: cwb_net.C:28
i drho i
void CompareWDMvsTime()
void plot(wavearray< double > &, char *=NULL, int=1, double=0., double=0., bool=false, float=0., float=0., bool=false, float=0., bool=false)
Definition: watplot.cc:150
gwavearray< double > * geWDM
TH2F * hist2D
Definition: watplot.hh:193
wavearray< double > w
Definition: Test1.C:27
virtual size_t size() const
Definition: wavearray.hh:145
TCanvas * canvas
Definition: watplot.hh:192
x plot
watplot * WTS
gwavearray< double > * gx
double fhigh
#define WF_HP
double GetTimeRange(double &tMin, double &tMax, double efraction=0.9999999)
Definition: gwavearray.cc:359
double dt
gwavearray< double > * greT
#define NZEROES
double flow
WDM< double > wdtf(lev, 2 *lev, 6, 10)
double T
Definition: testWDM_4.C:11
ifstream in
virtual void stop(double s)
Definition: wavearray.hh:139
double resolution(int=0)
Definition: wseries.hh:155
void resample(const wavearray< DataType_t > &, double, int=6)
Definition: wavearray.cc:503
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Definition: wseries.cc:246
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
TString gfile
watplot * GetWATPLOT()
Definition: gwavearray.hh:88
DataType_t * data
Definition: wavearray.hh:319
DataType_t getSample(int n, double m)
Definition: wseries.hh:185
double dT
Definition: testWDM_5.C:12
size_t sizeZero()
Definition: wseries.hh:144
int maxLayer()
Definition: wseries.hh:139
void Draw(GWAT_DRAW type=GWAT_TIME, TString options="ALP", Color_t color=kBlack)
Definition: gwavearray.cc:89