Logo coherent WaveBurst  
Library Reference Guide
Logo
RegressionL1.C
Go to the documentation of this file.
1 #define FRLIST_TARGET "input/L1_LDAS_C02_L2.frl"
2 #define FRLIST_WITNESS "input/L-R-9425.frl"
3 
4 #define nWITNESS 21
5 //#define nWITNESS 9
6 #define SCRATCH 8.
7 #define dFREQ 0.5
8 
9 #define LINE_60HZ
10 //#define LINE_120HZ
11 //#define LINE_180HZ
12 
13 #define GPS 942539008
14 #define LENGTH 616
15 
16 #define nFILTER 2
17 
18 using namespace CWB;
20 
21 void RegressionL1() {
22 
23  TString target = "L1:LDAS-STRAIN";
24 
25  TString witness[nWITNESS] = {
26  "L0:PEM-EX_MAGX",
27  "L1:SUS-MMT1_SENSOR_LL",
28  "L1:SUS-MMT1_SENSOR_LR",
29  "L1:SUS-MMT1_SENSOR_UR",
30  "L1:SUS-MMT1_SENSOR_UL",
31  "L1:SUS-MMT2_SENSOR_LL",
32  "L1:SUS-MMT2_SENSOR_LR",
33  "L1:SUS-MMT2_SENSOR_UR",
34  "L1:SUS-MMT2_SENSOR_UL",
35  "L1:SUS-MMT3_SENSOR_LL",
36  "L1:SUS-MMT3_SENSOR_LR",
37  "L1:SUS-MMT3_SENSOR_UR",
38  "L1:SUS-MMT3_SENSOR_UL",
39  "L1:SUS-ITMX_SENSOR_LL",
40  "L1:SUS-ITMX_SENSOR_LR",
41  "L1:SUS-ITMX_SENSOR_UR",
42  "L1:SUS-ITMX_SENSOR_UL",
43  "L1:SUS-ITMY_SENSOR_LL",
44  "L1:SUS-ITMY_SENSOR_LR",
45  "L1:SUS-ITMY_SENSOR_UR",
46  "L1:SUS-ITMY_SENSOR_UL"
47  };
48 
49 #ifdef LINE_60HZ
50  double fLow=50; double fHigh=70;
51 #endif
52 #ifdef LINE_120HZ
53  double fLow=110; double fHigh=130;
54 #endif
55 #ifdef LINE_180HZ
56  double fLow=170; double fHigh=190;
57 #endif
58 
59  // declare watplot
60  plot = new watplot(const_cast<char*>("plot"),200,20,800,500);
61  plot->gtitle(target);
62  plot->goptions(const_cast<char*>("alp logy"), 1, 0., 0., true, fLow, fHigh, true, 64);
63  //plot->goptions(const_cast<char*>("alp logy"), 1, 0., 0.);
64  //plot->gtitle("rms vs witness channel","witness","rms");
65 
66  // read target data
68  xx.start(GPS); xx.stop(GPS+LENGTH);
69 
70  frame frt(FRLIST_TARGET,target);
71  frt >> xx;
72 
73  // read witness data
76  for(int n=0;n<nWITNESS;n++) {
77  yy[n] = new wavearray<double>;
78  yy[n]->start(GPS); yy[n]->stop(GPS+LENGTH);
79  frw.setChName(witness[n]);
80  frw >> *yy[n];
81  cout << "yy : " << n << " " << yy[n]->size() << " " << yy[n]->rate() << endl;
82  }
83 
84  Biorthogonal<double> Bio(512);
85  WSeries<double> wT(Bio);
86  // resample target 16384Hz -> 512Hz
87  wT.Forward(xx,6);
88  wT.getLayer(xx,0);
89  // resample witness 0 2048Hz -> 512Hz
90  wT.Forward(*yy[0],3);
91  wT.getLayer(*yy[0],0);
92 
93  int level=xx.rate()/(2*dFREQ);
94  WDM<double> WD(level, level, 4, 12);
95 
96  wT.Forward(xx,WD);
97  wT.setlow(fLow);
98  wT.sethigh(fHigh);
99  regression rr(wT,const_cast<char*>("target"));
100  rr.add(*yy[0],const_cast<char*>("witness0"),59.,61.);
101  for(int n=1; n<nWITNESS; n++) rr.add(*yy[n],const_cast<char*>("witnessn"),0,3.);
102  for(int n=1; n<nWITNESS; n++) rr.add(1,n+1,const_cast<char*>("bi-witness"));
103  for(int n=1; n<nWITNESS; n++) rr.mask(n+1);
104 /*
105  for(int n=0;n<nWITNESS;n++) rr.add(*yy[n],"witness");
106  for(int n=2;n<=nWITNESS;n++) rr.add(1,n,"bi-witness");
107  for(int n=2;n<=nWITNESS;n++) rr.mask(n);
108 */
109  rr.setFilter(nFILTER);
110  rr.setMatrix(SCRATCH,1.);
111  rr.solve(0.,0,'h');
112  //rr.solve(0.,nWITNESS,'h');
113  rr.apply();
114 
115 /*
116  // extract rms vs witness channel
117  wavearray<double> rmsH = rr.getRRMS();
118  wavearray<double> rmsL = rr.getRRMS();
119  wavearray<double> rmsA = rr.getRRMS();
120 
121  rmsH.rate(1); rmsL.rate(1); rmsA.rate(1);
122  rmsH >> *plot; rmsL >> *plot; rmsA >> *plot;
123  return;
124 */
125 
126  // extract noise and clean data
127  wavearray<double> nn = rr.getNoise();
128  wavearray<double> cc = rr.getClean();
129 
130  xx >> *plot; cc >> *plot;
131 
132 #ifdef LINE_60HZ
133  plot->gtitle("L1 regression @ 60Hz");
134  *plot >> const_cast<char*>("L1_regression_60Hz.png");
135 #endif
136 #ifdef LINE_120HZ
137  plot->gtitle("L1 regression @ 120Hz");
138  *plot >> const_cast<char*>("L1_regression_1200Hz.png");
139 #endif
140 #ifdef LINE_180HZ
141  plot->gtitle("L1 regression @ 180Hz");
142  *plot >> const_cast<char*>("L1_regression_180Hz.png");
143 #endif
144 
145 }
void sethigh(double f)
Definition: wseries.hh:132
void gtitle(TString title="", TString xtitle="", TString ytitle="")
Definition: watplot.cc:1256
void setMatrix(double edge=0., double f=1.)
Definition: regression.cc:425
Definition: ced.hh:42
double fHigh
wavearray< double > target
void RegressionL1()
Definition: RegressionL1.C:21
int n
Definition: cwb_net.C:28
wavearray< double > getNoise()
Definition: regression.hh:141
#define nFILTER
Definition: RegressionL1.C:16
TString("c")
frame frt(FRLIST_NAME, HCHANNEL)
void setlow(double f)
Definition: wseries.hh:125
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
#define GPS
Definition: RegressionL1.C:13
void setChName(TString chName)
Definition: frame.hh:120
virtual size_t size() const
Definition: wavearray.hh:145
wavearray< double > xx
Definition: TestFrame1.C:11
#define FRLIST_TARGET
Definition: RegressionL1.C:1
void apply(double threshold=0., char c='a')
Definition: regression.cc:709
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
Definition: wseries.cc:193
void goptions(char *opt=NULL, int col=1, double t1=0., double t2=0., bool fft=false, float f1=0., float f2=0., bool psd=false, float t3=0., bool oneside=false)
Definition: watplot.cc:1221
void solve(double th, int nE=0, char c='s')
Definition: regression.cc:610
#define dFREQ
Definition: RegressionL1.C:7
size_t setFilter(size_t)
Definition: regression.cc:276
#define SCRATCH
Definition: RegressionL1.C:6
wavearray< double > getClean()
Definition: regression.hh:135
wavearray< double > yy
Definition: TestFrame5.C:12
watplot * plot
Definition: RegressionL1.C:19
double fLow
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
#define FRLIST_WITNESS
Definition: RegressionL1.C:2
void mask(int n, double flow=0., double fhigh=0.)
Definition: regression.cc:321
frame frw(FRLIST_NAME2)
#define nWITNESS
Definition: RegressionL1.C:4
#define LENGTH
Definition: RegressionL1.C:14