diff --git a/src/BayesWaveCleanFrame.c b/src/BayesWaveCleanFrame.c index b9dd82dea053dd1690665beda7c26cc471e9a322..4fbc7d88a4755f31c68b3cce7cf3236f1dd26762 100644 --- a/src/BayesWaveCleanFrame.c +++ b/src/BayesWaveCleanFrame.c @@ -195,6 +195,8 @@ int main(int argc, char *argv[]) start[j]=data->bw_model[j]; j++; } + //terminate the string + start[j]=0; //burn off next two digits j+=3; @@ -248,14 +250,14 @@ int main(int argc, char *argv[]) t_rise = 0.4; // the standard LIGO choice alpha = (2.0*t_rise/data->bw_seglen); - // Tukey filer the glitch reconstruction to avoid edge effects + // Tukey filter the glitch reconstruction to avoid edge effects tukey(glitch_model, alpha, data->bw_N); // FFT fftw_wrapper(glitch_model, data->bw_N, 1); // FFT scaling - norm = 2.0/(double)data->bw_N; + norm = 1.0/(double)data->bw_N; // zero pad the high frequencies for(i=0; i< data->bw_N; i++) @@ -298,7 +300,7 @@ int main(int argc, char *argv[]) //renormalize glitch model amplitude (undo windowing) tukey_scale(&s1, &s2, alpha, N); - norm = sqrt(1.0/(data->bw_dt*data->bw_seglen))*s2; + norm = sqrt(0.5/(data->bw_dt*data->bw_seglen))*s2; // Initialise and populate clean data for(i=0; i< N; i++) @@ -380,7 +382,11 @@ void fftw_wrapper(double *data, int N, int flag) int n; double *timeData = (double *)malloc(N*sizeof(double)); fftw_complex *freqData = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*N); - + + //setup FFTW plans + fftw_plan reverse = fftw_plan_dft_c2r_1d(N, freqData, timeData, FFTW_MEASURE); + fftw_plan forward = fftw_plan_dft_r2c_1d(N, timeData, freqData, FFTW_MEASURE); + switch (flag) { //Reverse transform @@ -397,31 +403,27 @@ void fftw_wrapper(double *data, int N, int flag) freqData[N/2][0] = freqData[0][1]; freqData[N/2][1] = freqData[0][1] = 0.0; - //setup FFTW plan - fftw_plan reverse = fftw_plan_dft_c2r_1d(N, freqData, timeData, FFTW_MEASURE); - + //The FFT fftw_execute(reverse); - fftw_destroy_plan(reverse); - + //Copy output back into data for(n=0; n<N; n++) data[n] = timeData[n]; - + break; //Forward transform case 1: - + //fill timeData with contents of data for(n=0; n<N; n++) timeData[n] = data[n]; - //setup FFTW plan - fftw_plan forward = fftw_plan_dft_r2c_1d(N, timeData, freqData, FFTW_MEASURE); + //fill timeData with contents of data + for(n=0; n<N; n++) timeData[n] = data[n]; //The FFT fftw_execute(forward); - fftw_destroy_plan(forward); - + //Copy output back into data for(n=0; n<N/2; n++) { @@ -432,7 +434,6 @@ void fftw_wrapper(double *data, int N, int flag) //get DC and Nyquist where FFTW wants them data[1] = freqData[N/2][1]; - break; default: @@ -440,7 +441,10 @@ void fftw_wrapper(double *data, int N, int flag) exit(1); break; } - + + fftw_destroy_plan(reverse); + fftw_destroy_plan(forward); + free(timeData); fftw_free(freqData); } diff --git a/src/BayesWaveIO.c b/src/BayesWaveIO.c index 9210255bf22891d6627bc5458cf330f90c97c764..d8b0d66e2c75c7e5fae5224e0566a7a27695a7f2 100644 --- a/src/BayesWaveIO.c +++ b/src/BayesWaveIO.c @@ -1096,10 +1096,7 @@ void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf, fftw_complex *hf = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*N); fftw_plan reverse = fftw_plan_dft_c2r_1d(N, hf, ht, FFTW_MEASURE); - for (i=0; i < N; i++) - { - ht[i] = 0.0; - } + for (i=0; i < N; i++) ht[i] = 0.0; hf[0][0] = 0.0; hf[0][1] = 0.0; @@ -1146,10 +1143,13 @@ void print_colored_time_domain_waveforms(char filename[], double *h, int N, doub int i; double t; FILE *waveout = fopen(filename,"w"); - - fftw_complex *hf = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*N); + double *ht = malloc(N*sizeof(double)); - + + fftw_complex *hf = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*N); + fftw_plan reverse = fftw_plan_dft_c2r_1d(N, hf, ht, FFTW_MEASURE); + + for(i=0; i<N; i++) ht[i] = 0.0; hf[0][0] = 0.0; @@ -1158,8 +1158,13 @@ void print_colored_time_domain_waveforms(char filename[], double *h, int N, doub { if(i>imin && i<imax) { - hf[i][0] = h[2*i+1]; - hf[i][1] = h[2*i]; + hf[i][0] = h[2*i]; + hf[i][1] = h[2*i+1]; + } + else + { + hf[i][0] = 0.0; + hf[i][1] = 0.0; } } @@ -1167,8 +1172,6 @@ void print_colored_time_domain_waveforms(char filename[], double *h, int N, doub hf[N/2][0] = hf[0][1]; hf[N/2][1] = hf[0][1] = 0.0; - - fftw_plan reverse = fftw_plan_dft_c2r_1d(N, hf, ht, FFTW_MEASURE); fftw_execute(reverse); double norm = sqrt((double)N); diff --git a/src/BayesWavePost.c b/src/BayesWavePost.c index 6e4ec4f363f2fcde1a996b968decc954d6f0dbcc..7d6043434b7ec7cc9a11b01afad1fc4b0a9e0d9d 100644 --- a/src/BayesWavePost.c +++ b/src/BayesWavePost.c @@ -842,18 +842,20 @@ int main(int argc, char *argv[]) /* */ /* Post process signal+glitch model */ /* */ - /***************************************/ - post_process_model("full", "signal", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); - post_process_model("full", "glitch", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); - post_process_model("full", "full", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); - + /***************************************/ + if(data->fullModelFlag) + { + post_process_model("full", "signal", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); + post_process_model("full", "glitch", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); + post_process_model("full", "full", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); + } /***************************************/ /* */ /* Post process signal model */ /* */ /***************************************/ - post_process_model("signal", "signal", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); + if(data->signalModelFlag) post_process_model("signal", "signal", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); /***************************************/ @@ -861,7 +863,7 @@ int main(int argc, char *argv[]) /* Post process glitch model */ /* */ /***************************************/ - post_process_model("glitch", "glitch", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); + if(data->glitchModelFlag) post_process_model("glitch", "glitch", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); /*****************************************/ @@ -869,7 +871,7 @@ int main(int argc, char *argv[]) /* Post process cleaning phase */ /* */ /*****************************************/ - post_process_model("clean", "glitch", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); + if(data->cleanModelFlag) post_process_model("clean", "glitch", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); /*****************************************/ @@ -877,7 +879,7 @@ int main(int argc, char *argv[]) /* Post process noise phase */ /* */ /*****************************************/ - post_process_model("noise", "noise", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); + if(data->noiseModelFlag) post_process_model("noise", "noise", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); /*****************************************/