Logo coherent WaveBurst  
Library Reference Guide
Logo
testWDM_8.C
Go to the documentation of this file.
1 { // accuracy of time delay filters:
2  // compare time shifts done in the TF domain using WDM TD filters vs the correct values
3  // WDM + TD Filters + WDM Inverse vs the trivial time shift of original data
4 
5  // TO REMOVE BOUNDARY ARTIFACTS, ONE SECOND OF DATA IS DISCARDED FOR EACH OF THE OPERATIONS:
6  // (i) APPLICATION OF TD FILTERS
7  // (ii) INVERSE TRANSFORMATION
8 
9  int Rate = 1024; // sampling rate (Hz)
10  int Duration = 20; // duration (seconds)
11  int N = Rate*Duration; // number of samples
12  wavearray<double> ts(N); //time series container
13  ts.rate(Rate);
14 
15  // time series is filled with white noise data:
16  TRandom3 rnd(0);
17  for(int i=0; i<N; i++) ts[i] = rnd.Gaus();
18 
19  WDM<double> wdm(32, 64, 4, 8); // define a WDM transform (32 bands)
20  WSeries<double> tfmap; // TF map container
21  tfmap.Forward(ts, wdm);
22 
23 
25  pWDM->setTDFilter(18); // computes TD filters
26 
27  double* tmp = new double[pWDM->nWWS/2]; // will store time delayed coefficients
28  for(int j=1*32; j<19*32; ++j) // time indices between 1-19s (1 sec discarded)
29  for(int i=0;i<=32; ++i) // frequency bands
30  tmp[j*33+i] = (double)pWDM->getPixelAmplitude(i, j, 13, false); //delay by 13 samples
31 
32  // overwrite original tmap:
33  double* TFMap = pWDM->pWWS; // all TF coefficients are stored in TFMap
34  for(int j=32; j<19*32; ++j) // time indices between 1-19s
35  for(int i=0;i<=32; ++i) // frequency bands
36  TFMap[j*33+i] = tmp[j*33+i];
37 
38  // inverse to time domain and plot shifted signal on top of original (ZOOM around 10s to see)
39  tfmap.Inverse();
40 
41  wavearray<double> ts2(N);
42  tfmap.getLayer(ts2, 0);
43 
44 
45  TH1F* h = new TH1F("h", "Time delay filter relative error", 100, -0.01, 0.01);
46  for(int j=1024*2; j<18*1024; ++j) // time indices between 2-18s (another second discarded)
47  h->Fill(ts2[j]/ts[j-13] - 1);
48  h->Draw();
49 
50  delete [] tmp;
51 
52  // NOTE: TD FILTERS TRUNCATION DEFINES THEIR PRECISION;
53  // FOR EXAMPLE, USING
54  // pWDM->setTDFilter(18)
55  // INSTEAD IMPROVES THE RELATIVE ERROR BY MORE THAN AN ORDER OF MAGNITUDE.
56 
57 }
virtual void rate(double r)
Definition: wavearray.hh:141
int Duration
Definition: testWDM_8.C:10
WDM< double > wdm(nLAYERS, nLAYERS, 6, 10)
DataType_t * pWWS
Definition: WaveDWT.hh:141
int j
Definition: cwb_net.C:28
i drho i
void setTDFilter(int nCoeffs, int L=1)
Definition: WDM.cc:639
double getPixelAmplitude(int, int, int, bool=false)
Definition: WDM.cc:748
wavearray< double > ts(N)
wavearray< double > h
Definition: Regression_H1.C:25
int N
Definition: testWDM_8.C:11
double * tmp
Definition: testWDM_5.C:31
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
Definition: wseries.cc:193
WSeries< double > tfmap
unsigned long nWWS
pointer to wavelet work space
Definition: WaveDWT.hh:142
WDM< double > * pWDM
Definition: testWDM_5.C:28
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Definition: wseries.cc:246
WaveDWT< DataType_t > * pWavelet
Definition: wseries.hh:456
TRandom3 rnd(0)
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
Definition: wseries.cc:291