Skip to content
Snippets Groups Projects
Commit 593918b1 authored by Jonah Kanner's avatar Jonah Kanner :nerd:
Browse files

Now output time domain waveform ... plot_pec waveform plotting is broken 

git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@133 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 1563744f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 134 additions and 112 deletions
src/BayesWavePostPEC.c
+ 134
− 112
View file @ 593918b1
......@@ -299,7 +299,7 @@ int main(int argc, char *argv[])
fprintf(momentsOut[ifo], "# overlap energy_inj energy_rec time_inj time_rec duration_inj duration_rec f_inj f_rec band_inj ban_rec \n");
}
// opening files for moments in whitened data
// Open output files for moments in whitened data
FILE **w_momentsOut = malloc(NI*sizeof(FILE*));
for(ifo=0; ifo<NI; ifo++){
sprintf(filename,"w_moments.dat.%i", ifo);
......@@ -307,6 +307,14 @@ int main(int argc, char *argv[])
fprintf(w_momentsOut[ifo], "# overlap energy_inj energy_rec time_inj time_rec duration_inj duration_rec f_inj f_rec band_inj ban_rec \n");
}
// Open output files for reconstructed time domain waveform (colored noise)
FILE **t_recOut = malloc(NI*sizeof(FILE*));
for(ifo=0; ifo<NI; ifo++){
sprintf(filename,"t_recOut.dat.%i", ifo);
t_recOut[ifo] = fopen(filename,"w");
fprintf(t_recOut[ifo], "# Each line is h(t) for one sample in the chain");
}
mkdir("injected_waveforms",S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
int mc=0;
......@@ -324,9 +332,10 @@ int main(int argc, char *argv[])
double *ty=dvector(0,data->N-1);
int test=1;
int sampcount=0;
while(!feof(signalParams))
{
sampcount++;
/******************************************************************************/
/* */
/* Get current state of MCMC output into format for BayesWave code library */
......@@ -398,118 +407,130 @@ int main(int argc, char *argv[])
double white;
for(ifo=0; ifo<NI; ifo++) {
for(white=0; white<2; white++) {
hihi = innerproduct(imin, N/2, hinj[ifo], hinj[ifo], 1.0, psd[ifo]);
hrhr = innerproduct(imin, N/2, hrec[ifo], hrec[ifo], 1.0, psd[ifo]);
hrhi = innerproduct(imin, N/2, hrec[ifo], hinj[ifo], 1.0, psd[ifo]);
double overlap;
overlap = hrhi/sqrt(hihi*hrhr);
// Calculate signal energy frequency distribution (aka rhof)
// Calculate frequency moments for the reconstructed signal
double f0_rec, band_rec, energy_rec;
double rhof_rec[data->imax];
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, band_inj, energy_inj;
double rhof_inj[data->imax];
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);
// 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);
// Calculate the time moments for the injected signal
double t0_inj, dur_inj, t_energy_inj;
double rhot_inj[N];
t_energy_inj = get_rhot(N, hoft_inj, deltaT, rhot_inj);
t0_inj = get_f0(N, rhot_inj, deltaT);
dur_inj = get_band(N, rhot_inj, deltaT, t0_inj);
// fprintf(overlapsOut[ifo], "# energy_inj energy_rec time_inj time_rec duration_inj duration_rec f_inj f_rec band_inj ban_rec \n");
if (white == 0) {
fprintf(momentsOut[ifo], "%e %e %e %e %e %e %e %e %e %e %e \n", overlap, energy_inj, energy_rec, t0_inj, t0_rec, dur_inj, dur_rec, f0_inj, f0_rec, band_inj, band_rec);
} else {
fprintf(w_momentsOut[ifo], "%e %e %e %e %e %e %e %e %e %e %e \n", overlap, energy_inj, energy_rec, t0_inj, t0_rec, dur_inj, dur_rec, f0_inj, f0_rec, band_inj, band_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;
// Plot the freq distribution of the singal energy
for (i=0;i<data->imax;i++) {
freq = deltaF*i;
fprintf(rhoout, "%e %e \n", freq, rhof_rec[i]);
}
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;
fprintf(hoftout, "%e %e \n", t, hoft_rec[i]);
for(white=0; white<2; white++) {
hihi = innerproduct(imin, N/2, hinj[ifo], hinj[ifo], 1.0, psd[ifo]);
hrhr = innerproduct(imin, N/2, hrec[ifo], hrec[ifo], 1.0, psd[ifo]);
hrhi = innerproduct(imin, N/2, hrec[ifo], hinj[ifo], 1.0, psd[ifo]);
double overlap;
overlap = hrhi/sqrt(hihi*hrhr);
// Calculate signal energy frequency distribution (aka rhof)
// Calculate frequency moments for the reconstructed signal
double f0_rec, band_rec, energy_rec;
double rhof_rec[data->imax];
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, band_inj, energy_inj;
double rhof_inj[data->imax];
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);
// 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);
// Calculate the time moments for the injected signal
double t0_inj, dur_inj, t_energy_inj;
double rhot_inj[N];
t_energy_inj = get_rhot(N, hoft_inj, deltaT, rhot_inj);
t0_inj = get_f0(N, rhot_inj, deltaT);
dur_inj = get_band(N, rhot_inj, deltaT, t0_inj);
// fprintf(overlapsOut[ifo], "# energy_inj energy_rec time_inj time_rec duration_inj duration_rec f_inj f_rec band_inj ban_rec \n");
if (white == 0) {
fprintf(momentsOut[ifo], "%e %e %e %e %e %e %e %e %e %e %e \n", overlap, energy_inj, energy_rec, t0_inj, t0_rec, dur_inj, dur_rec, f0_inj, f0_rec, band_inj, band_rec);
} else {
fprintf(w_momentsOut[ifo], "%e %e %e %e %e %e %e %e %e %e %e \n", overlap, energy_inj, energy_rec, t0_inj, t0_rec, dur_inj, dur_rec, f0_inj, f0_rec, band_inj, band_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;
// Plot the freq distribution of the singal energy
for (i=0;i<data->imax;i++) {
freq = deltaF*i;
fprintf(rhoout, "%e %e \n", freq, rhof_rec[i]);
}
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;
fprintf(hoftout, "%e %e \n", t, hoft_rec[i]);
}
fclose(hoftout);
// 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 / %e \n", t_energy_inj, 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 / %e seconds \n", t0_inj, t0_rec);
fprintf(stdout, "The injected/recovered duration is: %e / %e seconds \n", dur_inj, dur_rec);
}
// ------- End Sanity checks ---------
// Output some sweet waveforms
if((sampcount % 1000) == 0) {
if( white == 0 ){
for(i=0; i<N; i++) {
fprintf(t_recOut[ifo], "%e ", hoft_rec[i]);
}
fprintf(t_recOut[ifo], "\n");
}
}
// Francesco's code for outputting injected waveform in same way as gnu plot flag
// Possible replsced by new big output files
// if(mc>=frame*frameInterval && mc<M)
//{
// sprintf(filename,"injected_waveforms/signal_%i_wave_%d.dat.%i", (int)tempMin, frame, ifo);
// for(i=0; i< N; i++) ty[i] = hinj[ifo][i];
// print_time_domain_waveforms(filename, ty, N, data->Snf[ifo], 1.0, data->Tobs, fix, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
// if(ifo==NI-1) frame++;
//}
// Whiten the data before the next iteration of the whiten loop
whiten_data(imin, N/2, hrec[ifo], psd[ifo]);
whiten_data(imin, N/2, hinj[ifo], psd[ifo]);
}
}
fclose(hoftout);
// 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 / %e \n", t_energy_inj, 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 / %e seconds \n", t0_inj, t0_rec);
fprintf(stdout, "The injected/recovered duration is: %e / %e seconds \n", dur_inj, dur_rec);
}
// ------- End Sanity checks ---------
if(mc>=frame*frameInterval && mc<M)
{
sprintf(filename,"injected_waveforms/signal_%i_wave_%d.dat.%i", (int)tempMin, frame, ifo);
for(i=0; i< N; i++) ty[i] = hinj[ifo][i];
// print_time_domain_waveforms(filename, ty, N, data->Snf[ifo], 1.0, data->Tobs, fix, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
if(ifo==NI-1) frame++;
}
// Whiten the data before the next iteration of the whiten loop
whiten_data(imin, N/2, hrec[ifo], psd[ifo]);
whiten_data(imin, N/2, hinj[ifo], psd[ifo]);
}
}
mc+=chain->cycle;
}
}
......@@ -525,6 +546,7 @@ int main(int argc, char *argv[])
for(ifo=0;ifo<NI;ifo++){
fclose(momentsOut[ifo]);
fclose(w_momentsOut[ifo]);
fclose(t_recOut[ifo]);
}
return 0;
......
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