There will be a short amount of downtime, for git.ligo.org, docs.ligo.org, and chat.ligo.org, starting around approximately 10am CDT on Tuesday 18th June 2019. This is to enable access controls for GitLab Pages. More information can be found here.

Commit b5359b23 authored by James Clark's avatar James Clark

Merge branch 'freq_at_max' into 'master'

Adding the frequency at maximum amplitude as a moment

See merge request !52
parents 3e85e9d3 9578a2ed
Pipeline #46721 failed with stages
in 30 seconds
......@@ -45,6 +45,7 @@ struct Moments
double networkOverlap;
double *hmax;
double *t_at_hmax;
double *f_at_max_amp;
double *medianfrequency;
double *mediantime;
double *bandwidth90percentinterval;
......@@ -107,6 +108,9 @@ void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char
int post_process_model(char model_name[], char model_type[], struct Data *data, struct Model *model, struct Chain *chain, double **hdata, double **hoft, double **psd, double **one, double **hinj, int **dimension, int injectionFlag, LIGOTimeGPS GPStrig, bool lite);
void get_time_domain_envelope(double *envelope, double *hdata_sample, int N, int imin, int imax);
double get_time_at_max(double *hdata_sample, int N, int imin, int imax, double deltaT);
void get_freq_at_max(double *hdata_sample, int N, int imin, int imax, double deltaT, double *f_at_max_amp);
/* ============================ MAIN PROGRAM ============================ */
......@@ -1199,6 +1203,15 @@ void ComputeWaveformMoments(struct Data *data, double **h, double **psd, struct
get_hmax_and_t(N, wt[ifo], deltaT, &m->hmax[ifo], &m->t_at_hmax[ifo]);
}
/*
Frquency at maximum amplitude
*/
if(ifo<NI)
{
get_freq_at_max(wf[ifo], N, imin, imax, deltaT, &m->f_at_max_amp[ifo]);
}
}
......@@ -1351,6 +1364,7 @@ void initialize_moments(struct Moments *m, int NI)
m->snr = malloc((NI+1)*sizeof(double));
m->hmax = malloc((NI+1)*sizeof(double));
m->t_at_hmax = malloc((NI+1)*sizeof(double));
m->f_at_max_amp = malloc((NI+1)*sizeof(double));
m->medianfrequency = malloc((NI+1)*sizeof(double));
m->bandwidth90percentinterval = malloc((NI+1)*sizeof(double));
m->mediantime = malloc((NI+1)*sizeof(double));
......@@ -1370,6 +1384,7 @@ void free_moments(struct Moments *m)
free(m->snr);
free(m->hmax);
free(m->t_at_hmax);
free(m->f_at_max_amp);
free(m->medianfrequency);
free(m->bandwidth90percentinterval);
free(m->mediantime );
......@@ -1411,6 +1426,7 @@ void setup_output_files(int NI, int injectionFlag, char *outdir, char *model, ch
}
fprintf(momentsOut[0], "h_max ");
fprintf(momentsOut[0], "t_at_h_max ");
fprintf(momentsOut[0], "f_at_max_amp");
fprintf(momentsOut[0], "\n");
}
......@@ -1446,6 +1462,7 @@ void setup_output_files(int NI, int injectionFlag, char *outdir, char *model, ch
}
fprintf(momentsOut[ifo], "h_max ");
fprintf(momentsOut[ifo], "t_at_h_max ");
fprintf(momentsOut[ifo], "f_at_max_amp");
fprintf(momentsOut[ifo], "\n");
}
}
......@@ -1492,6 +1509,7 @@ void print_moments(int ifo, int injectionFlag, FILE **outfile, struct Moments *m
}
fprintf(outfile[ifo], "%e ",m->hmax[ifo]);
fprintf(outfile[ifo], "%e ",m->t_at_hmax[ifo]);
fprintf(outfile[ifo], "%e ",m->f_at_max_amp[ifo]);
fprintf(outfile[ifo], "\n");
}
......@@ -2033,6 +2051,7 @@ void print_stats(struct Data *data, LIGOTimeGPS GPS, int injectionFlag, char *ou
double *bandwidth = NULL;
double *hmax = NULL;
double *t_at_hmax = NULL;
double *f_at_max_amp = NULL;
int seconds = GPS.gpsSeconds;
int nanoseconds = GPS.gpsNanoSeconds;
......@@ -2117,6 +2136,7 @@ void print_stats(struct Data *data, LIGOTimeGPS GPS, int injectionFlag, char *ou
fprintf(statsfile, "bandwidth ");
fprintf(statsfile, "h_max ");
fprintf(statsfile, "t_at_h_max ");
fprintf(statsfile, "f_at_max_amp ");
fprintf(statsfile, "\n");
count=0;
......@@ -2142,7 +2162,7 @@ void print_stats(struct Data *data, LIGOTimeGPS GPS, int injectionFlag, char *ou
bandwidth = malloc(count*sizeof(double));
hmax = malloc(count*sizeof(double));
t_at_hmax = malloc(count*sizeof(double));
f_at_max_amp = malloc(count*sizeof(double));
//ignore header again
fgets(filerow,10000,momentsfile);
......@@ -2163,6 +2183,7 @@ void print_stats(struct Data *data, LIGOTimeGPS GPS, int injectionFlag, char *ou
}
fscanf(momentsfile, "%lg",&hmax[i]);
fscanf(momentsfile, "%lg",&t_at_hmax[i]);
fscanf(momentsfile, "%lg",&f_at_max_amp[i]);
}
//sort each moment array
......@@ -2173,6 +2194,7 @@ void print_stats(struct Data *data, LIGOTimeGPS GPS, int injectionFlag, char *ou
gsl_sort(bandwidth,1,count);
gsl_sort(hmax,1,count);
gsl_sort(t_at_hmax,1,count);
gsl_sort(f_at_max_amp,1,count);
//get and print medians for stats file
fprintf(statsfile, "%i ", map );
......@@ -2184,6 +2206,7 @@ void print_stats(struct Data *data, LIGOTimeGPS GPS, int injectionFlag, char *ou
fprintf(statsfile, "%g ", gsl_stats_median_from_sorted_data (bandwidth,1,count) );
fprintf(statsfile, "%g ", gsl_stats_median_from_sorted_data (hmax,1,count));
fprintf(statsfile, "%g ", gsl_stats_median_from_sorted_data (t_at_hmax,1,count) + (GPStime + 2.0 - data->Tobs));
fprintf(statsfile, "%g ", gsl_stats_median_from_sorted_data (f_at_max_amp,1,count) );
fprintf(statsfile, "\n");
fclose(momentsfile);
......@@ -3464,3 +3487,69 @@ void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char
}
void get_time_domain_envelope(double *envelope, double *hdata_sample, int N, int imin, int imax)
{
//Allocate memory for hoft_original, hdata_shifted, and hoft_shifted
double *hoft_original = malloc(N*sizeof(double));
double *hoft_shifted = malloc(N*sizeof(double));
double *hdata_shifted = malloc(N*sizeof(double));
//Get time domain waveform
get_time_domain_waveforms(hoft_original, hdata_sample, N, imin, imax);
//Shift frequency domain waveform by pi/2
for(int i = 0; i<(int)(N/2.0); i++){
hdata_shifted[2*i] = -1*hdata_sample[2*i+1];
hdata_shifted[2*i+1] = hdata_sample[2*i];
}
//Get time domain waveform after shift
get_time_domain_waveforms(hoft_shifted, hdata_shifted, N, imin, imax);
//Calculate the envelope
for(int i=0; i<N; i++){
envelope[i]=sqrt(hoft_original[i]*hoft_original[i]+hoft_shifted[i]*hoft_shifted[i]);
}
free(hoft_original);
free(hoft_shifted);
free(hdata_shifted);
}
double get_time_at_max(double *hdata_sample, int N, int imin, int imax, double deltaT)
{
//grab the time domain envelope
double *amplitude_envelope = malloc(N*sizeof(double));
get_time_domain_envelope(amplitude_envelope, hdata_sample, N, imin, imax);
//find the iteration of the maximum amplitude
int max_iter = 0;
double max_time = 0;
double max_amp = 0;
for(int i=0; i<N; i++){
if(amplitude_envelope[i]>max_amp){
max_amp = amplitude_envelope[i];
max_iter = i;
}
}
//find the time of the maximum amplitude
max_time = (double)max_iter*deltaT;
free(amplitude_envelope);
return(max_time);
}
void get_freq_at_max(double *hdata_sample, int N, int imin, int imax, double deltaT, double *f_at_max_amp)
{
//grab the time at maximum amplitude
double max_time = get_time_at_max(hdata_sample, N, imin, imax, deltaT);
//grab the instantaneous frequency at that time
double *instantaneous_freq = malloc(N*sizeof(double));
double *hoft_sample = malloc(N*sizeof(double));
get_time_domain_waveforms(hoft_sample, hdata_sample, N, imin, imax);
tf_tracks(hoft_sample, instantaneous_freq, deltaT, N);
*f_at_max_amp = instantaneous_freq[(int)(max_time/deltaT)];
free(instantaneous_freq);
free(hoft_sample);
}
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