Correctly normalized signal energy

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

src/BayesWavePostPEC.c

@@ -22,6 +22,7 @@ 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);
 int energy(int imin, int imax, double *a, double deltaF, double *arg);
+double rhonorm(int imax, double *a, double deltaF);
 
 /* ============================  MAIN PROGRAM  ============================ */
 
@@ -142,8 +143,9 @@ int main(int argc, char *argv[])
   double deltaT = dataPtr->timeData->deltaT;
 
   // JBK - About to guess delta f !!
-  double deltaF = fmax / (2*N); 
-
+  double deltaF = fmax / (N/2); 
+  fprintf(stdout, "Max freq is %e", fmax);
+  
   int tsize = dataPtr->timeData->data->length/4;
 
   double TwoDeltaToverN = 2.0 * deltaT / ((double) dataPtr->timeData->data->length);
@@ -305,6 +307,7 @@ int main(int argc, char *argv[])
   //work-space for chain diagnostics (periodic time-domain waveforms)
   double *ty=dvector(0,data->N-1);
           
+  int test=0;
   while(!feof(signalParams))
   {
 
@@ -389,10 +392,25 @@ int main(int argc, char *argv[])
 
 
 	// Try calculating signal energy
-	// rhorec = malloc(imax*sizeof(double));
 	double rhorec[data->imax];
 	energy(data->imin, data->imax, hrec[ifo], deltaF, rhorec);
 
+	// Try plotting the signal energy
+	double freq;
+	if (test == 0) {
+	  FILE *rhoout = fopen("rho.txt","w");
+	  for (i=0;i<data->imax;i++) {
+	    freq = deltaF*i;
+	    fprintf(rhoout, "%e %e \n", freq, rhorec[i]);
+	  }
+	  test=1;
+	  fclose(rhoout);
+	}
+
+	// Sanity check that rho is normalized to 1
+	// double testrho;
+	// testrho = rhonorm(data->imax, rhorec, deltaF);
+	// fprintf(stdout, "%e\n", testrho);
 
         if(mc>=frame*frameInterval && mc<M)
         {
@@ -490,7 +508,7 @@ int energy(int imin, int imax, double *a, double deltaF, double *arg)
   {
     norm += (a[i*2]*a[i*2] + a[i*2+1]*a[i*2+1]);
   }
-  norm *= 2.0*deltaF;
+  norm *= deltaF;
 
   // Calculate energy
   for(i=0; i<imax; i++)
@@ -501,12 +519,27 @@ int energy(int imin, int imax, double *a, double deltaF, double *arg)
       } 
     else
       {
-	arg[i] = a[i*2]*a[i*2] + a[i*2+1]*a[i*2+1];
+	arg[i] = (a[i*2]*a[i*2] + a[i*2+1]*a[i*2+1])/norm;
       }
   }
   return(0);
 }
 
+
+double rhonorm(int imax, double *a, double deltaF)
+{
+  int i;
+  double norm;
+  // Calculate normalization
+  norm = 0.0;
+  for(i=0; i<imax; i++)
+  {
+    norm += (a[i]);
+  }
+  norm *= deltaF;
+  return(norm);
+}
+
 double innerproduct(int imin, int imax, double *a, double *b, double Sn, double *Snf)
 {
  	int i;

src/Makefile

@@ -6,7 +6,7 @@ CCFLAGS = -Wall -g -std=gnu99 -O3 -pedantic
 
 OBJS = Subroutines.o BayesLine.o
 
-all: $(OBJS) BayesWaveBurst BayesWavePost
+all: $(OBJS) BayesWaveBurst BayesWavePostPEC
 
 BayesLine.o : BayesLine.c BayesLine.h
 	gcc $(CCFLAGS) -c BayesLine.c $($INCDIR:%=-I%)
@@ -17,9 +17,9 @@ Subroutines.o : Subroutines.c BayesLine.c BayesLine.h Declarations.h Constants.h
 BayesWaveBurst: BayesWaveBurst.c Declarations.h Constants.h numrec.h  $(OBJS)
 	gcc $(CCFLAGS) -o BayesWaveBurst BayesWaveBurst.c $(OBJS) $(LIBDIR:%=-L%) $(INCDIR:%=-I%) $(LIBS:%=-l%)
 
-BayesWavePost: BayesWavePost.c Declarations.h Constants.h numrec.h  $(OBJS)
-	gcc $(CCFLAGS) -o BayesWavePost BayesWavePost.c $(OBJS) $(LIBDIR:%=-L%) $(INCDIR:%=-I%) $(LIBS:%=-l%)
+BayesWavePostPEC: BayesWavePostPEC.c Declarations.h Constants.h numrec.h  $(OBJS)
+	gcc $(CCFLAGS) -o BayesWavePostPEC BayesWavePostPEC.c $(OBJS) $(LIBDIR:%=-L%) $(INCDIR:%=-I%) $(LIBS:%=-l%)
 
 clean:
-	rm *.o BayesWaveBurst BayesWavePost
+	rm *.o BayesWaveBurst BayesWavePostPEC