Commit 30f62edd authored by James Clark's avatar James Clark

Merge branch 'post-no-clean' into 'master'

More stable application of signal+glitch model

See merge request !61
parents d011d35f 385b9e52
Pipeline #47201 passed with stages
in 2 minutes and 53 seconds
...@@ -195,6 +195,8 @@ int main(int argc, char *argv[]) ...@@ -195,6 +195,8 @@ int main(int argc, char *argv[])
start[j]=data->bw_model[j]; start[j]=data->bw_model[j];
j++; j++;
} }
//terminate the string
start[j]=0;
//burn off next two digits //burn off next two digits
j+=3; j+=3;
...@@ -248,14 +250,14 @@ int main(int argc, char *argv[]) ...@@ -248,14 +250,14 @@ int main(int argc, char *argv[])
t_rise = 0.4; // the standard LIGO choice t_rise = 0.4; // the standard LIGO choice
alpha = (2.0*t_rise/data->bw_seglen); 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); tukey(glitch_model, alpha, data->bw_N);
// FFT // FFT
fftw_wrapper(glitch_model, data->bw_N, 1); fftw_wrapper(glitch_model, data->bw_N, 1);
// FFT scaling // FFT scaling
norm = 2.0/(double)data->bw_N; norm = 1.0/(double)data->bw_N;
// zero pad the high frequencies // zero pad the high frequencies
for(i=0; i< data->bw_N; i++) for(i=0; i< data->bw_N; i++)
...@@ -298,7 +300,7 @@ int main(int argc, char *argv[]) ...@@ -298,7 +300,7 @@ int main(int argc, char *argv[])
//renormalize glitch model amplitude (undo windowing) //renormalize glitch model amplitude (undo windowing)
tukey_scale(&s1, &s2, alpha, N); 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 // Initialise and populate clean data
for(i=0; i< N; i++) for(i=0; i< N; i++)
...@@ -380,7 +382,11 @@ void fftw_wrapper(double *data, int N, int flag) ...@@ -380,7 +382,11 @@ void fftw_wrapper(double *data, int N, int flag)
int n; int n;
double *timeData = (double *)malloc(N*sizeof(double)); double *timeData = (double *)malloc(N*sizeof(double));
fftw_complex *freqData = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*N); 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) switch (flag)
{ {
//Reverse transform //Reverse transform
...@@ -397,31 +403,27 @@ void fftw_wrapper(double *data, int N, int flag) ...@@ -397,31 +403,27 @@ void fftw_wrapper(double *data, int N, int flag)
freqData[N/2][0] = freqData[0][1]; freqData[N/2][0] = freqData[0][1];
freqData[N/2][1] = freqData[0][1] = 0.0; 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 //The FFT
fftw_execute(reverse); fftw_execute(reverse);
fftw_destroy_plan(reverse);
//Copy output back into data //Copy output back into data
for(n=0; n<N; n++) data[n] = timeData[n]; for(n=0; n<N; n++) data[n] = timeData[n];
break; break;
//Forward transform //Forward transform
case 1: case 1:
//fill timeData with contents of data //fill timeData with contents of data
for(n=0; n<N; n++) timeData[n] = data[n]; for(n=0; n<N; n++) timeData[n] = data[n];
//setup FFTW plan //fill timeData with contents of data
fftw_plan forward = fftw_plan_dft_r2c_1d(N, timeData, freqData, FFTW_MEASURE); for(n=0; n<N; n++) timeData[n] = data[n];
//The FFT //The FFT
fftw_execute(forward); fftw_execute(forward);
fftw_destroy_plan(forward);
//Copy output back into data //Copy output back into data
for(n=0; n<N/2; n++) for(n=0; n<N/2; n++)
{ {
...@@ -432,7 +434,6 @@ void fftw_wrapper(double *data, int N, int flag) ...@@ -432,7 +434,6 @@ void fftw_wrapper(double *data, int N, int flag)
//get DC and Nyquist where FFTW wants them //get DC and Nyquist where FFTW wants them
data[1] = freqData[N/2][1]; data[1] = freqData[N/2][1];
break; break;
default: default:
...@@ -440,7 +441,10 @@ void fftw_wrapper(double *data, int N, int flag) ...@@ -440,7 +441,10 @@ void fftw_wrapper(double *data, int N, int flag)
exit(1); exit(1);
break; break;
} }
fftw_destroy_plan(reverse);
fftw_destroy_plan(forward);
free(timeData); free(timeData);
fftw_free(freqData); fftw_free(freqData);
} }
......
...@@ -1096,10 +1096,7 @@ void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf, ...@@ -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_complex *hf = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*N);
fftw_plan reverse = fftw_plan_dft_c2r_1d(N, hf, ht, FFTW_MEASURE); fftw_plan reverse = fftw_plan_dft_c2r_1d(N, hf, ht, FFTW_MEASURE);
for (i=0; i < N; i++) for (i=0; i < N; i++) ht[i] = 0.0;
{
ht[i] = 0.0;
}
hf[0][0] = 0.0; hf[0][0] = 0.0;
hf[0][1] = 0.0; hf[0][1] = 0.0;
...@@ -1146,10 +1143,13 @@ void print_colored_time_domain_waveforms(char filename[], double *h, int N, doub ...@@ -1146,10 +1143,13 @@ void print_colored_time_domain_waveforms(char filename[], double *h, int N, doub
int i; int i;
double t; double t;
FILE *waveout = fopen(filename,"w"); FILE *waveout = fopen(filename,"w");
fftw_complex *hf = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*N);
double *ht = malloc(N*sizeof(double)); 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; for(i=0; i<N; i++) ht[i] = 0.0;
hf[0][0] = 0.0; hf[0][0] = 0.0;
...@@ -1158,8 +1158,13 @@ void print_colored_time_domain_waveforms(char filename[], double *h, int N, doub ...@@ -1158,8 +1158,13 @@ void print_colored_time_domain_waveforms(char filename[], double *h, int N, doub
{ {
if(i>imin && i<imax) if(i>imin && i<imax)
{ {
hf[i][0] = h[2*i+1]; hf[i][0] = h[2*i];
hf[i][1] = 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 ...@@ -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][0] = hf[0][1];
hf[N/2][1] = hf[0][1] = 0.0; 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); fftw_execute(reverse);
double norm = sqrt((double)N); double norm = sqrt((double)N);
......
...@@ -842,18 +842,20 @@ int main(int argc, char *argv[]) ...@@ -842,18 +842,20 @@ int main(int argc, char *argv[])
/* */ /* */
/* Post process signal+glitch model */ /* Post process signal+glitch model */
/* */ /* */
/***************************************/ /***************************************/
post_process_model("full", "signal", data, model, chain, hdata, hoft, psd, one, hinj, dimension, injectionFlag, GPStrig, lite); if(data->fullModelFlag)
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_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 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[]) ...@@ -861,7 +863,7 @@ int main(int argc, char *argv[])
/* Post process glitch model */ /* 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[]) ...@@ -869,7 +871,7 @@ int main(int argc, char *argv[])
/* Post process cleaning phase */ /* 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[]) ...@@ -877,7 +879,7 @@ int main(int argc, char *argv[])
/* Post process noise phase */ /* 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);
/*****************************************/ /*****************************************/
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment