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