Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
STFT.cc
Go to the documentation of this file.
1 #include "STFT.hh"
2 #include "TPaletteAxis.h"
3 
4 CWB::STFT::STFT(wavearray<double> x, int nfft, int noverlap, TString ztype,
5  TString fwindow, double fparam, TString name) :
6  isLogz(false), title("") {
7 
8  this->name = name;
9  this->canvas = NULL;
10  this->h2 = NULL;
11  this->ztype = ztype;
12  int iztype=0;
13  if(ztype.CompareTo("amplitude")==0) iztype=1;
14  if(ztype.CompareTo("energy")==0) iztype=2;
15  if(iztype==0) {cout << "CWB::STFT::STFT: Error wrong ztype (amplitude/energy)" << endl;exit(1);}
16 
17  if (nfft<=0) {cout << "CWB::STFT::STFT: Error nfft must be positive" << endl;exit(1);}
18  if (noverlap<=0||noverlap>=nfft) {cout << "CWB::STFT::STFT: Error noverlap must be > 0 & < nfft" << endl;exit(1);}
19  if (x.rate()<=0) {cout << "CWB::STFT::STFT: Error sample rate must be > 0" << endl;exit(1);}
20  if (x.start()<0) {cout << "CWB::STFT::STFT: Error start time must be positive" << endl;exit(1);}
21  if (int(x.size())<nfft) {cout << "CWB::STFT::STFT: Error size must be > nfft " << endl;exit(1);}
22 
23  //cout << "nfft : " << nfft << endl;
24  //cout << "sample rate : " << x.rate() << endl;
25  //cout << "start time : " << x.start() << endl;
26 
27  // Compute Window
28  CWB::Window wnd(const_cast<char*>(fwindow.Data()),nfft,fparam);
29  window = new double[nfft];
30  for (int i=0;i<nfft;i++) window[i] = wnd.GetValue(i);
31 
32  double df=(double)x.rate()/(double)(nfft);
33  int loops = x.size()/nfft;
34 
35  int nshift=nfft-noverlap;
36  int N=int(nfft/nshift);
37  double Tmin=x.start();
38  double Tmax=x.size()/x.rate()+x.start();
39  int nT=int(x.size()/nshift);
40  double Fmin=0.;
41  double Fmax=x.rate()/2.;
42  int nF=nfft/2;
43  h2 = new TH2D(name.Data(),"Spectrogram", nT, Tmin, Tmax, nF, Fmin, Fmax);
44  wavearray<double> y(nfft);
45 //double tEN=0;
46 //double fEN=0;
47  double dt=1./x.rate();
48  double norm=sqrt(x.rate());
49  for (int n=0;n<N*loops;n++) {
50  int shift=n*nshift;
51  if(shift+nfft>int(x.size())) break;
52  for (int i=0;i<nfft;i++) y.data[i]=norm*x.data[i+shift]*window[i];
53 //for (int i=0;i<nfft;i++) tEN+=pow(y.data[i],2);
54  y.FFT(1);
55  double time = (shift+nfft/2)*dt;
56  for (int i=0;i<nfft;i+=2) {
57  double psd=sqrt(pow(y.data[i],2)+pow(y.data[i+1],2))*sqrt(1/df);
58 //fEN+=pow(psd,2);
59  double frequency=(double)i*df/2.;
60  if(iztype==2) psd*=psd; // energy
61  h2->SetBinContent(int(nT*time/(Tmax-x.start())),int(nF*frequency/Fmax),psd);
62  }
63 //cout << tEN/x.rate() << " " << fEN*df << " " << df << endl;exit(0);
64  }
65  y.resize(0);
66 }
67 
68 CWB::STFT::~STFT() {
69 
70  if(h2!=NULL) delete h2;
71  if(canvas!=NULL) delete canvas;
72  if(window!=NULL) delete [] window;
73 }
74 
75 void
76 CWB::STFT::Draw(double t1, double t2, double f1, double f2, double z1, double z2,
77  int dpaletteId, Option_t* goption) {
78 
79 // TCanvas* tcanvas = (TCanvas*) gROOT->FindObject(name);
80 // if (tcanvas!=NULL) {cout << "CWB::STFT::STFT: Error Canvas " << name.Data() << " already exist" << endl;exit(1);}
81  if(canvas!=NULL) delete canvas;
82  //canvas = new TCanvas(name, "LVC experiment", 300,40, 1000, 600);
83  canvas = new TCanvas(name, "LVC experiment", 300,40, 800, 600);
84  canvas->Clear();
85  canvas->ToggleEventStatus();
86  canvas->SetGridx(false);
87  canvas->SetGridy(false);
88  canvas->SetLogz(isLogz);
89  canvas->SetFillColor(kWhite);
90  canvas->SetRightMargin(0.10);
91  canvas->SetLeftMargin(0.10);
92  canvas->SetBottomMargin(0.13);
93  canvas->SetBorderMode(0);
94  //canvas->SetWindowSize(1200,600);
95  //canvas->SetGrayscale();
96 
97  // remove the red box around canvas
98  gStyle->SetFrameBorderMode(0);
99  gROOT->ForceStyle();
100 
101  gStyle->SetTitleH(0.050);
102  gStyle->SetTitleW(0.95);
103  gStyle->SetTitleY(0.98);
104  gStyle->SetTitleFont(12,"D");
105  gStyle->SetTitleColor(kBlue,"D");
106  gStyle->SetTextFont(12);
107  gStyle->SetTitleFillColor(kWhite);
108  gStyle->SetLineColor(kWhite);
109  gStyle->SetNumberContours(256);
110 // gStyle->SetMarkerStyle(7);
111 // gStyle->SetMarkerSize(2);
112  gStyle->SetCanvasColor(kWhite);
113  gStyle->SetStatBorderSize(1);
114 
115  if (dpaletteId==DUMMY_PALETTE_ID) {
116  if (paletteId!=0) {
117  SetPlotStyle(paletteId);
118  } else {
119  gStyle->SetPalette(1,0);
120  }
121  } else {
122  if (dpaletteId!=0) {
123  SetPlotStyle(dpaletteId);
124  } else {
125  gStyle->SetPalette(1,0);
126  }
127  }
128 
129  canvas->cd();
130  canvas->SetLogz(isLogz);
131 
132  h2->SetStats(kFALSE);
133  h2->SetTitleFont(12);
134  h2->SetTitle(title);
135  h2->SetFillColor(kWhite);
136 
137  h2->GetXaxis()->SetNdivisions(506);
138  h2->GetXaxis()->SetLabelFont(42);
139  h2->GetXaxis()->SetLabelOffset(0.014);
140  h2->GetXaxis()->SetTitleOffset(1.4);
141  h2->GetYaxis()->SetTitleOffset(1.2);
142  h2->GetYaxis()->SetNdivisions(506);
143  h2->GetYaxis()->SetLabelFont(42);
144  h2->GetYaxis()->SetLabelOffset(0.01);
145  h2->GetZaxis()->SetLabelFont(42);
146  h2->GetZaxis()->SetNoExponent(false);
147  h2->GetZaxis()->SetNdivisions(506);
148 
149  h2->GetXaxis()->SetTitleFont(42);
150  h2->GetXaxis()->SetTitle("Time (sec)");
151  h2->GetXaxis()->CenterTitle(true);
152 
153  h2->GetYaxis()->SetTitleFont(42);
154  h2->GetYaxis()->SetTitle("Frequency (Hz)");
155  h2->GetYaxis()->CenterTitle(true);
156 
157  //char ztitle[256];
158  //sprintf(ztitle,"Normalized tile %s",ztype.Data());
159 
160  h2->GetZaxis()->SetTitleOffset(0.6);
161  h2->GetZaxis()->SetTitleFont(42);
162  //h2->GetZaxis()->SetTitle(ztitle);
163  h2->GetZaxis()->CenterTitle(true);
164 
165  h2->GetXaxis()->SetLabelSize(0.03);
166  h2->GetYaxis()->SetLabelSize(0.03);
167  h2->GetZaxis()->SetLabelSize(0.03);
168 
169  if(title.Sizeof()==1) {
170  char stitle[256];
171  sprintf(stitle,"Spectrogram (Normalized tile %s)",ztype.Data());
172  h2->SetTitle(stitle);
173  }
174  double dt=h2->GetXaxis()->GetBinWidth(0);
175  double df=h2->GetYaxis()->GetBinWidth(0);
176 
177  int nt1=0;
178  int nt2=h2->GetNbinsX();
179  int nf1=0;
180  int nf2=h2->GetNbinsY();
181 
182  if(t2>t1) nt1=int((t1-h2->GetXaxis()->GetXmin())/dt);
183  if(t2>t1) nt2=int((t2-h2->GetXaxis()->GetXmin())/dt);
184  if(f2>f1) nf1=int(f1/df);
185  if(f2>f1) nf2=int(f2/df);
186 
187  double h2max=0.0;
188  for (int i=nt1;i<nt2;i++) {
189  for (int j=nf1;j<nf2;j++) {
190  double binc=h2->GetBinContent(i,j);
191  if(binc>h2max) h2max=binc;
192  }
193  }
194  //h2->GetZaxis()->SetRangeUser(0,int(h2max)+1);
195  h2->GetZaxis()->SetRangeUser(0,h2max);
196 
197  if(t2>t1) h2->GetXaxis()->SetRangeUser(t1,t2);
198  if(f2>f1) h2->GetYaxis()->SetRangeUser(f1,f2);
199  if(z2>z1) h2->GetZaxis()->SetRangeUser(z1,z2);
200 
201  h2->Draw(goption);
202 
203  // change palette's width
204  canvas->Update();
205  TPaletteAxis *palette = (TPaletteAxis*)h2->GetListOfFunctions()->FindObject("palette");
206  palette->SetX1NDC(0.91);
207  palette->SetX2NDC(0.933);
208  palette->SetTitleOffset(0.92);
209  palette->GetAxis()->SetTickSize(0.01);
210  canvas->Modified();
211 
212 }
213 
214 // -----------------------------------------------------------------
215 // http://ultrahigh.org/2007/08/20/making-pretty-root-color-palettes/
216 // -----------------------------------------------------------------
217 void
218 CWB::STFT::SetPlotStyle(int paletteId) {
219 
220  const Int_t NRGBs = 5;
221  const Int_t NCont = 255;
222 
223  if (fabs(paletteId)==1) {
224  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
225  Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
226  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
227  Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
228  if (paletteId<0) {
229  TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
230  } else {
231  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
232  }
233  } else
234  if (fabs(paletteId)==2) {
235  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
236  Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
237  //Double_t red[NRGBs] = { 0.00, 0.00, 0.00, 1.00, 1.00 };
238  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
239  //Double_t green[NRGBs] = { 0.00, 1.00, 1.00, 1.00, 0.00 };
240  Double_t blue[NRGBs] = { 1.00, 1.00, 0.00, 0.00, 0.00 };
241  if (paletteId<0) {
242  TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
243  } else {
244  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
245  }
246  } else
247  if (fabs(paletteId)==3) {
248  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
249  Double_t red[NRGBs] = { 0.00, 0.09, 0.18, 0.09, 0.00 };
250  Double_t green[NRGBs] = { 0.01, 0.02, 0.39, 0.68, 0.97 };
251  Double_t blue[NRGBs] = { 0.17, 0.39, 0.62, 0.79, 0.97 };
252  if (paletteId<0) {
253  TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
254  } else {
255  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
256  }
257  } else
258  if (fabs(paletteId)==4) {
259  Double_t stops[NRGBs] = { 0.00, 0.50, 0.75, 0.875, 1.00 };
260  Double_t red[NRGBs] = { 1.00, 1.00, 1.00, 1.00, 1.00 };
261  Double_t green[NRGBs] = { 1.00, 0.75, 0.50, 0.25, 0.00 };
262  Double_t blue[NRGBs] = { 0.00, 0.00, 0.00, 0.00, 0.00 };
263  if (paletteId<0) {
264  TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
265  } else {
266  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
267  }
268  } else
269  if (fabs(paletteId)==5) { // Greyscale palette
270  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
271  Double_t red[NRGBs] = { 1.00, 0.84, 0.61, 0.34, 0.00 };
272  Double_t green[NRGBs] = { 1.00, 0.84, 0.61, 0.34, 0.00 };
273  Double_t blue[NRGBs] = { 1.00, 0.84, 0.61, 0.34, 0.00 };
274  if (paletteId<0) {
275  TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
276  } else {
277  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
278  }
279  } else
280  if (fabs(paletteId)==57) { // kBird palette (is defined only in ROOT6)
281  Double_t stops[9] = { 0., 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0};
282  Double_t red[9] = { 0.2082, 0.0592, 0.0780, 0.0232, 0.1802, 0.5301, 0.8186, 0.9956, 0.9764};
283  Double_t green[9] = { 0.1664, 0.3599, 0.5041, 0.6419, 0.7178, 0.7492, 0.7328, 0.7862, 0.9832};
284  Double_t blue[9] = { 0.5293, 0.8684, 0.8385, 0.7914, 0.6425, 0.4662, 0.3499, 0.1968, 0.0539};
285  if (paletteId<0) {
286  TColor::CreateGradientColorTable(9, stops, blue, green, red, NCont);
287  } else {
288  TColor::CreateGradientColorTable(9, stops, red, green, blue, NCont);
289  }
290  }
291  gStyle->SetNumberContours(NCont);
292 
293  return;
294 }
295 
296 void
297 CWB::STFT::Print(TString pname) {
298 
299  canvas->Print(pname);
300 
301 /*
302  if(TString(pname).Contains(".png")!=0) { // fix gray background for png plots
303  TString gname=pname;
304  gname.ReplaceAll(".png",".gif");
305  canvas->Print(gname);
306  char cmd[1024];
307  sprintf(cmd,"convert %s %s",gname.Data(),pname.Data());
308  //cout << cmd << endl;
309  gSystem->Exec(cmd);
310  sprintf(cmd,"rm %s",gname.Data());
311  //cout << cmd << endl;
312  gSystem->Exec(cmd);
313  } else {
314  canvas->Print(pname);
315  }
316 */
317 }
NRGBs
#define NRGBs
wavearray::size
virtual size_t size() const
Definition: wavearray.hh:127
noverlap
int noverlap
Definition: TestDelta.C:20
green
Double_t green[nRGBs]
Definition: cwb_report_prod_1.C:78
CWB::STFT::ztype
TString ztype
Definition: STFT.hh:78
wavearray::rate
virtual void rate(double r)
Definition: wavearray.hh:123
CWB::STFT::~STFT
~STFT()
Definition: STFT.cc:68
CWB::STFT::Print
void Print(TString pname)
Definition: STFT.cc:297
name
par[0] name
Definition: CWB_Plugin_HEN_SIM_Config.C:87
n
int n
Definition: cwb_net.C:10
const_cast< char * >
cout<< "skymap size : "<< L<< endl;for(int l=0;l< L;l++) sm.set(l, l);sm > const_cast< char * >("skymap.dat")
TString
TString("c")
Definition: cwb_report_skymap.C:111
palette
int palette
Definition: DrawGnetwork2.C:17
x
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
Definition: BaseFunctionWDM.C:51
frequency
double frequency
Definition: DrawClusterWithSkyplot.C:20
psd
wavearray< double > psd(33)
nfft
int nfft
Definition: TestDelta.C:19
canvas
return wmap canvas
Definition: DrawClusterWithSkyplot.C:80
blue
Double_t blue[nRGBs]
Definition: cwb_report_prod_1.C:79
Fmin
double Fmin
Definition: TestDelta.C:27
CWB::STFT::h2
TH2D * h2
Definition: STFT.hh:72
wavearray::start
virtual void start(double s)
Definition: wavearray.hh:119
j
int j
Definition: cwb_net.C:10
i
i drho i
Definition: cwb_epparameters.C:85
N
#define N
CWB::Window::GetValue
double GetValue(unsigned i)
Definition: Window.cc:59
CWB::STFT::Draw
void Draw(double t1=0.0, double t2=0.0, double f1=0.0, double f2=0.0, double z1=0.0, double z2=0.0, int dpaletteId=DUMMY_PALETTE_ID, Option_t *option="colfz")
Definition: STFT.cc:76
nT
nT
Definition: cbc_plots.C:659
NCont
#define NCont
DUMMY_PALETTE_ID
#define DUMMY_PALETTE_ID
Definition: gskymap.hh:47
stops
Double_t stops[nRGBs]
Definition: cwb_report_prod_1.C:76
time
double time[6]
Definition: cbc_plots.C:435
CWB::STFT::canvas
TCanvas * canvas
Definition: STFT.hh:71
Tmax
double Tmax
Definition: TestDelta.C:26
CWB::STFT::SetPlotStyle
void SetPlotStyle(int paletteId=1)
Definition: STFT.cc:218
Tmin
double Tmin
Definition: TestDelta.C:25
int
i() int(T_cor *100))
f2
TF1 * f2
Definition: cbc_plots.C:1710
dt
double dt
Definition: DrawPhaseShift.C:13
title
char title[256]
Definition: SSeriesExample.C:1
CWB::STFT::window
double * window
Definition: STFT.hh:80
fabs
double fabs(const Complex &x)
Definition: numpy.cc:37
Fmax
double Fmax
Definition: TestDelta.C:28
df
double df
Definition: HowToAccessToWSeries.C:36
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
wavearray::data
DataType_t * data
Definition: wavearray.hh:301
CWB::STFT::STFT
STFT(wavearray< double > x, int nfft, int noverlap, TString ztype="amplitude", TString fwindow="hann", double fparam=0.0, TString name="stft")
Definition: STFT.cc:4
t1
TH1 * t1
Definition: cbc_plots_sim4_ifar.C:1896
CWB::STFT::name
TString name
Definition: STFT.hh:76
shift
double shift[NIFO_MAX]
Definition: cwb2G_parameters.C:358
y
wavearray< double > y
Definition: Test10.C:31
red
Double_t red[nRGBs]
Definition: cwb_report_prod_1.C:77
exit
exit(0)