Now computing both freq and time domain moments

git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@119 c56465c9-8126-4a4f-9d7d-ac845eff4865

postprocess/s6/plot_pec.py

@@ -241,6 +241,19 @@ for ifo in ifoList:
     plt.savefig('sig_energy.png')
 
 
+    # --------------------------
+    # For debugging, plot signal energy (time domain)
+    # ----------------------
+    plt.figure()
+    t, e = np.loadtxt('t_rho.txt', unpack=True)
+    plt.plot(t,e)
+    plt.grid()
+    plt.ylabel('Signal energy')
+    plt.xlabel('Time')
+    plt.xlim(1.98, 2.02)
+    plt.savefig('t_sig_energy.png')
+
+
     # --------------------------
     # For debugging, plot time domain waveform
     # ----------------------

src/BayesWavePostPEC.c

@@ -22,11 +22,12 @@ void print_help_message(void);
 double innerproduct(int imin, int imax, double *a, double *b, double Sn, double *Snf);
 void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf, double eta, double Tobs, double fix, int imin, int imax, double tmin, double tmax);
 double get_energy(int imin, int imax, double *a, double deltaF);
-int get_rhof(int imin, int imax, double *a, double deltaF, double *arg);
+double get_rhof(int imin, int imax, double *a, double deltaF, double *arg);
 double rhonorm(int imax, double *a, double deltaF);
 double get_f0(int imax, double *a, double deltaF);
 double get_band(int imax, double *a, double deltaF, double f0);
 void get_time_domain_waveforms(double *hoft, double *h, int N, double *Snf, double eta, double Tobs, double fix, int imin, int imax, double tmin, double tmax);
+double get_rhot(int N, double *a, double deltaT, double *arg);
 
 /* ============================  MAIN PROGRAM  ============================ */
 
@@ -397,30 +398,38 @@ int main(int argc, char *argv[])
 
 	// Calculate signal energy frequency distribution (aka rhof)	
 	// Calculate frequency moments for the reconstructed signal
-	double f0_rec;
-	double band_rec;
+	double f0_rec, band_rec, energy_rec;
 	double rhof_rec[data->imax];
-	get_rhof(data->imin, data->imax, hrec[ifo], deltaF, rhof_rec);
+	energy_rec = get_rhof(data->imin, data->imax, hrec[ifo], deltaF, rhof_rec);
 	f0_rec = get_f0(data->imax, rhof_rec, deltaF);
 	band_rec = get_band(data->imax, rhof_rec, deltaF, f0_rec);
 	
 	// Calculate frequency moments for the injected signal
-	double f0_inj;
-	double band_inj;
+	double f0_inj, band_inj, energy_inj;
 	double rhof_inj[data->imax];
-	get_rhof(data->imin, data->imax, hinj[ifo], deltaF, rhof_inj);
+	energy_inj = get_rhof(data->imin, data->imax, hinj[ifo], deltaF, rhof_inj);
 	f0_inj = get_f0(data->imax, rhof_inj, deltaF);
 	band_inj = get_band(data->imax, rhof_inj, deltaF, f0_inj);
 
-	// ------ Sanity Checks! ------------- 
-	if (test == 1) {
-
 	// Calculate the time domain waveform
 	double hoft_rec[N];
+	double hoft_inj[N];
 	get_time_domain_waveforms(hoft_rec, hrec[ifo], N, data->Snf[ifo], 1.0, data->Tobs, fix, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
+	get_time_domain_waveforms(hoft_inj, hinj[ifo], N, data->Snf[ifo], 1.0, data->Tobs, fix, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
 
+	// Calculate the time moments for the reconstructed signal
+	double t0_rec, dur_rec, t_energy_rec;
+	double rhot_rec[N];
+	t_energy_rec = get_rhot(N, hoft_rec, deltaT, rhot_rec);
+	t0_rec =  get_f0(N, rhot_rec, deltaT);
+	dur_rec = get_band(N, rhot_rec, deltaT, t0_rec);
+		
+
+	// ------ Sanity Checks! ------------- 
+	if (test == 1) {
 	  double freq, t;
 	  FILE *rhoout = fopen("rho.txt","w");
+	  FILE *t_rhoout = fopen("t_rho.txt","w");
 	  FILE *hoftout = fopen("sanity_time.txt","w");
 	  double testrho;
 	  double energy;
@@ -431,6 +440,13 @@ int main(int argc, char *argv[])
 	  }	  
 	  fclose(rhoout);
 
+	  // Plot the time distribution of the singal energy
+	  for (i=0;i<N;i++) {
+	    t = i*deltaT;
+	    fprintf(t_rhoout, "%e %e \n", t, rhot_rec[i]);
+	  }	  
+	  fclose(t_rhoout);
+
 	  // Plot the time domain waveform
 	  for (i=0;i<N;i++) {
 	    t = i*deltaT;
@@ -441,12 +457,17 @@ int main(int argc, char *argv[])
 	  // Sanity check that rho is normalized to 1	  
 	  testrho = rhonorm(data->imax, rhof_rec, deltaF);
 	  fprintf(stdout, "The normalization of rhof is: %e\n", testrho);
-
+	  testrho = rhonorm(N, rhot_rec, deltaT);
+	  fprintf(stdout, "The normalization of rhot is: %e \n", testrho);
 	  energy = get_energy(data->imin, data->imax, hrec[ifo], deltaF);
 	  fprintf(stdout, "The total signal energy is: %e \n", energy);
 	  test=0;
+	  fprintf(stdout, "The time domain signal energy is: %e \n", t_energy_rec);
+	  fprintf(stdout, "The injected/recovered signal energy is: %e / %e  \n", energy_inj, energy_rec);
 	  fprintf(stdout, "The injected/recovered central frequency is: %e / %e Hz \n", f0_inj, f0_rec);
 	  fprintf(stdout, "The injected/recovered bandwidth is: %e / %e Hz \n", band_inj, band_rec);
+	  fprintf(stdout, "The injected/recovered central time is: %e seconds \n", t0_rec);
+	  fprintf(stdout, "The injected/recovered duration is: %e seconds \n", dur_rec);
 	}
 	// ------- End Sanity checks ---------
 
@@ -548,9 +569,10 @@ double get_energy(int imin, int imax, double *a, double deltaF)
   return(result);
 }
   
-int get_rhof(int imin, int imax, double *a, double deltaF, double *arg)
+double get_rhof(int imin, int imax, double *a, double deltaF, double *arg)
 {
   // Calculate rho_f (Normalized signal energy frequency distribution.  See moments.pdf eqn 1.5)
+  
   int i;
   double norm;
 
@@ -563,9 +585,31 @@ int get_rhof(int imin, int imax, double *a, double deltaF, double *arg)
       arg[i] = 2*(a[i*2]*a[i*2] + a[i*2+1]*a[i*2+1])/norm;
     }
   }
-  return(0);
+  // Return value is the total signal energy
+  return(norm);
 }
 
+double get_rhot(int N, double *a, double deltaT, double *arg)
+{
+  // Calculate rho_t (Normalized signal energy frequency distribution.  See moments.pdf eqn 1.5)
+  int i;
+  double norm;
+
+  // We'll normalize by the total signal energy
+  norm = 0;
+  for(i=0;i<N; i++) {
+    norm += a[i]*a[i];
+  }
+  norm *= deltaT;
+
+  for(i=0; i<N; i++) {
+    arg[i] = a[i]*a[i]/norm;
+  }
+  // Return value is the total signal energy
+  return(norm);
+}
+
+
 double rhonorm(int imax, double *a, double deltaF)
 {
   // Calculate normalization.  This is a sanity check, since rho is supposed to have a normalization of 1.
@@ -615,12 +659,13 @@ double innerproduct(int imin, int imax, double *a, double *b, double Sn, double
   double arg;
 
   arg = 0.0;
-	for(i=imin; i<imax; i++)
+  for(i=imin; i<imax; i++)
   {
-		arg += (a[i*2]*b[i*2] + a[i*2+1]*b[i*2+1])/Snf[i];
+    //arg += (a[i*2]*b[i*2] + a[i*2+1]*b[i*2+1])/Snf[i];
+    arg += (a[i*2]*b[i*2] + a[i*2+1]*b[i*2+1]);
   }
-
-  arg *= 4.0/Sn;
+  // arg *= 4.0/Sn;
+  arg *= 4.0;
 
 	return(arg);
 }
@@ -680,30 +725,30 @@ void get_time_domain_waveforms(double *hoft, double *h, int N, double *Snf, doub
   double x;
   // double x,t;
 
-  h[0] = 0.0;
-  h[1] = 0.0;
+  hoft[0] = 0.0;
+  hoft[1] = 0.0;
   for(i=1; i< N/2; i++)
   {
     if(i>imin && i<imax)
     {
       // x = sqrt(Snf[i]*eta); // Dividing by the Snf should whiten the waveform
       x = sqrt(eta);        //JBK - try making an UNWHITENED time domain waveform
-      h[2*i] /= x;
-      h[2*i+1] /= -x;
+      hoft[2*i] = h[2*i] / x;
+      hoft[2*i+1] = h[2*i+1]/(-x);
     }
     else
     {
-      h[2*i]=h[2*i+1]=0.0;
+      hoft[2*i]=hoft[2*i+1]=0.0;
     }
   }
 
-  drealft(h-1,N,-1);
+  drealft(hoft-1,N,-1);
 
   // double t0 = Tobs-2.0;
 
   for(i=0; i<N; i++) {
     // t = (double)(i)/(double)(N)*Tobs;
     //    if(t>=tmin && t<tmax) fprintf(waveout,"%e %e\n", t-t0, h[i]/fix);
-    hoft[i] = h[i]/fix;
+    hoft[i] = hoft[i]/fix;
   }
 }