Logo coherent WaveBurst  
Library Reference Guide
Logo
Test2.C
Go to the documentation of this file.
1 //Test 2: Similar, but witness channel includes some uncorrelated noise:
2 //Generate white noise time series w(t), x(t), y(t)
3 //Witness channel: A(t) = 0.6*y(t) + w(t)
4 //Target channel: H(t) = x(t) + 0.8*w(t)
5 //Compare the ASD (or simply the RMS) of H(t) before and after regression.
6 //Since regression removes the w(t) correlated noise but adds in some y(t)
7 //noise from the witness channel, I think the regression should be able to
8 //reduce the noise amplitude by a factor of
9 // sqrt(1^2+(0.8*0.6)^2)/sqrt(1^2+0.8^2) = sqrt(1.23)/sqrt(1.64) ~ 0.87
10 
11 {
12  //Time
13  #define LENGHT 1200
14  #define SCRATCH 32
15  #define RATE 2048
16 
17  using namespace CWB;
18 
19  //define x channel properties -> gauss 1
20  int N = RATE*(LENGHT+2*SCRATCH);
22  x.rate(RATE);
23  x.resize(N);
24  x.start(0);
25  x.stop(LENGHT+2*SCRATCH);
26 
27  //define w channel properties -> gauss 2
30 
31  // time series is filled with white noise data:
32  TRandom3 rnd(0);
33  for(int i=0; i<N; i++) x[i] = rnd.Gaus();
34  for(int i=0; i<N; i++) w[i] = rnd.Gaus();
35  for(int i=0; i<N; i++) y[i] = rnd.Gaus();
36 
37  //Fill target and witness
38  wavearray<double> A = w; y*= 0.6; A += y;
39  wavearray<double> H = x; w*= 0.8; H += w;
40 
41  //Make WDM transform, resolution = 1Hz
42  int lev=H.rate()/2;
43  WDM<double> wdtf(lev, 2*lev, 6, 10);
45  tfmap.Forward(H, wdtf);
46 
47  //Adding target channel
48  regression r;
49  r.add(tfmap,const_cast<char*>("hchannel"));
50 
51  //Adding witness channel
52  r.add(A,const_cast<char*>("witness"));
53 
54  //Calculate prediction
55  r.setFilter(5);
56  r.setMatrix(SCRATCH);
57  r.solve(0, 20, 'h');
58  r.apply(0.2);
59 
60  wavearray<double> clean=r.getClean();
61  cout << "Ratio rms: (" << clean.rms() << "/" << H.rms() << ")= " << clean.rms()/H.rms() << endl;
62 
63  cout << "x : " << x.mean() << " " << x.rms() << endl;
64  cout << "clean : " << clean.mean() << " " << clean.rms() << endl;
65  clean -= x;
66  cout << "clean-x : " << clean.mean() << " " << clean.rms() << endl;
67 
68  wavearray<double> eigen=r.getVEIGEN(-1);
69  eigen.rate(11);
70  watplot plot;
71  TCanvas *c1 = plot.canvas;
72  c1->Divide(1,2);
73 
74  c1->cd(1);
75  plot.plot(eigen,const_cast<char*>("alp"),1);
76  plot.graph[0]->SetTitle("Eigen-values of all layers");
77 
80  c1->cd(2);
81  TGraph* g=new TGraph(freq.size(),freq.data,rank.data);
82  g->SetLineColor(1);
83  g->Draw("alp");
84  g->SetTitle("Ranking for all layers");
85 
86  c1->SaveAs("Test2.png");
87 
88 
89  exit(0);
90 
91 }
wavearray< double > x
Definition: Test2.C:21
static double g(double e)
Definition: GNGen.cc:116
void setMatrix(double edge=0., double f=1.)
Definition: regression.cc:425
wavearray< double > y
Definition: Test2.C:29
Definition: ced.hh:42
virtual void rate(double r)
Definition: wavearray.hh:141
wavearray< double > w
Definition: Test2.C:28
#define RATE
wavearray< double > rank
Definition: Regression_H1.C:80
std::vector< TGraph * > graph
Definition: watplot.hh:194
size_t add(WSeries< double > &target, char *name, double fL=0., double fH=0.)
Definition: regression.cc:91
virtual void start(double s)
Definition: wavearray.hh:137
i drho i
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
virtual double rms()
Definition: wavearray.cc:1206
virtual size_t size() const
Definition: wavearray.hh:145
TCanvas * canvas
Definition: watplot.hh:192
#define SCRATCH
TRandom3 rnd(0)
wavearray< double > freq
Definition: Regression_H1.C:79
x plot
TCanvas * c1
void apply(double threshold=0., char c='a')
Definition: regression.cc:709
wavearray< double > getRank(int n)
Definition: regression.hh:152
WSeries< double > tfmap
#define LENGHT
void solve(double th, int nE=0, char c='s')
Definition: regression.cc:610
wavearray< double > vfreq
Definition: regression.hh:189
static double A
Definition: geodesics.cc:26
size_t setFilter(size_t)
Definition: regression.cc:276
wavearray< double > getVEIGEN(int n=-1)
Definition: regression.cc:357
wavearray< double > getClean()
Definition: regression.hh:135
regression r
Definition: Regression_H1.C:44
WDM< double > wdtf(lev, 2 *lev, 6, 10)
int N
Definition: Test2.C:20
virtual double mean() const
Definition: wavearray.cc:1071
virtual void stop(double s)
Definition: wavearray.hh:139
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Definition: wseries.cc:246
DataType_t * data
Definition: wavearray.hh:319
int lev
Definition: Regression_H1.C:38
virtual void resize(unsigned int)
Definition: wavearray.cc:463
exit(0)