diff --git a/src/BayesWavePost.c b/src/BayesWavePost.c
index b7087b2622e97e85a9cc52dd78df59fba75e73c4..c50c9dfddbe934b5279e26b1dc6e4e1a7f06c499 100644
--- a/src/BayesWavePost.c
+++ b/src/BayesWavePost.c
@@ -3108,26 +3108,28 @@ int post_process_model(char model_name[], char model_type[], struct Data *data,
   // Calculate median
   for (ifo=0; ifo<NI; ifo++)
   {
+    //Compute the median hoff waveform                                                                                                                                  
+    median_and_CIs(waveform_mega[ifo], N, Nwaveform, deltaF, medianWaveformFile[ifo], hoft_median[ifo]);
     
-    // Print fourier domain residuals from a random draw from the chain
-    sprintf(filename,"%s/fourier_domain_%s_residual_%s.dat",outdir,model_type,data->ifos[ifo]);
+    // Print fourier domain residuals from the median waveform reconstruction                                                                                           
+    sprintf(filename,"%s/fourier_domain_%s_median_residual_%s.dat",outdir,model_type,data->ifos[ifo]);
     residualFile[ifo] = fopen(filename,"w");
-    
-    
+
     if(glitchFlag || signalFlag || fullFlag)
-    {
-      for (ii = data->imin; ii < data->N/2; ii++)
       {
-        residual[ifo][2*ii]   -= waveform_mega[ifo][random_draw][2*ii];
-        residual[ifo][2*ii+1] -= waveform_mega[ifo][random_draw][2*ii+1];
+	for (ii = data->imin; ii < data->N/2; ii++)
+	  {
+	    residual[ifo][2*ii]   -= hoft_median[ifo][2*ii];
+	    residual[ifo][2*ii+1] -= hoft_median[ifo][2*ii+1];
+	  }
       }
-    }
-    
+
     for (ii = 0; ii < data->N/2; ii++)
-    {
-      fprintf(residualFile[ifo],"%g %g %g\n", (double)ii*deltaF, residual[ifo][2*ii], residual[ifo][2*ii+1]);
-    }
+      {
+	fprintf(residualFile[ifo],"%g %g %g\n", (double)ii*deltaF, residual[ifo][2*ii], residual[ifo][2*ii+1]);
+      }
     fclose(residualFile[ifo]);
+
     
     // Calc and print median PSD
     median_and_CIs(psd_mega[ifo], N/2, Nwaveform, deltaF, medianPSDFile[ifo], psd_median[ifo]);
@@ -3154,9 +3156,8 @@ int post_process_model(char model_name[], char model_type[], struct Data *data,
         }
       }
       median_and_CIs(psd_mega[ifo], N/2, Nwaveform, deltaF, medianFrequencyWaveformFile[ifo], hoff_median[ifo]);
-      
-      
-      
+
+
       // Get (whitened) time domain waveforms, calc a print median
       for (n = 0; n < Nwaveform; n++) {
         whiten_data(data->imin, data->imax, waveform_mega[ifo][n], psd[ifo]);
@@ -3200,7 +3201,7 @@ int post_process_model(char model_name[], char model_type[], struct Data *data,
     fclose(whitenedGeocenterFile[0]);
     fclose(coloredGeocenterFile[0]);
   }
-  
+
   for(ifo=0; ifo<NI; ifo++)
   {
     if(data->bayesLineFlag)
@@ -3348,7 +3349,7 @@ void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char
   Nsamp=0;
 
   //get residual file
-  sprintf(filename,"%s/fourier_domain_%s_residual_%s.dat",outdir,model,ifos[ifo]);
+  sprintf(filename,"%s/fourier_domain_%s_median_residual_%s.dat",outdir,model,ifos[ifo]);
   dFile = fopen(filename,"r");
 
   while(!feof(dFile))
@@ -3367,7 +3368,7 @@ void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char
   
   // combine Fourier coefficints
   //get residual file
-  sprintf(filename,"%s/fourier_domain_%s_residual_%s.dat",outdir,model,ifos[ifo]);
+  sprintf(filename,"%s/fourier_domain_%s_median_residual_%s.dat",outdir,model,ifos[ifo]);
   dFile = fopen(filename,"r");
   
   //get PSD file
@@ -3379,7 +3380,7 @@ void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char
   {
     fscanf(dFile,"%lg %lg %lg",&f,&dre,&dim);
     fscanf(nFile,"%lg %lg %lg %lg %lg %lg",&f,&Sn,&Sn25,&Sn75,&Sn05,&Sn95);
-    if(f>fmin)
+    if(f>=fmin)
     {
       samples[2*nn]   = dre/sqrt(Sn/(Tobs/2.));
       samples[2*nn+1] = dim/sqrt(Sn/(Tobs/2.));