Added median PSD estimation method
authorTito Dal Canton <tito@dalcanton.it>
Mon, 17 Dec 2012 10:09:57 +0000 (11:09 +0100)
committerTito Dal Canton <tito@dalcanton.it>
Mon, 17 Dec 2012 13:41:46 +0000 (14:41 +0100)
Fixed --spectrum-type options in cbc/findchirp and cbc/spectrum
Added --sample-rate option to cbc/findchirp
Other minor fixes

cbc/findchirp/include/findchirp-cli.hpp
cbc/findchirp/src/findchirp-cli.cpp
cbc/findchirp/src/findchirp.cpp
cbc/spectrum/src/spectrum-cli.cpp
common/include/filter/GWFilter.hpp

index fb767c4..15ba412 100644 (file)
@@ -52,7 +52,7 @@ namespace GW {
     const char* get_bank_file() const {return m_var_map["bank-file"].as<std::string>().c_str();}
     float       get_snr_threshold() const {return m_var_map["snr-threshold"].as<float>();}
     int         get_perform_chisq2() const {return m_var_map["perform-chisq2"].as<int>();}
-    float       get_chisq2_threshold() const {return m_var_map["chisq2-threshold"].as<float>();}
+    float       get_chisq_threshold() const {return m_var_map["chisq-threshold"].as<float>();}
     int         get_chisq_bins() const {return m_var_map["chisq-bins"].as<int>();}
     const char* get_clustering_algo() const {return m_var_map["clustering-algo"].as<std::string>().c_str();}
     const char* get_coincidence_algo() const {return m_var_map["coincidence-algo"].as<std::string>().c_str();}
@@ -63,6 +63,8 @@ namespace GW {
     int         get_dynamic_exp() const {return m_var_map["dynamic-range-exponent"].as<int>();}
     int         get_segment_size() const {return m_var_map["segment-length"].as<int>();}
     int         get_segment_overlap() const {return m_var_map["segment-overlap"].as<int>();}
+    int         get_sample_rate() const {return m_var_map["sample-rate"].as<int>();}
+    std::string get_spectrum_type() const {return m_var_map["spectrum-type"].as<std::string>();}
 
   }; // CLIFindchirp
 
index 0036c05..6471e2e 100644 (file)
@@ -29,7 +29,7 @@ GW::CLIFindchirp::CLIFindchirp(int ac, char* av[]) :
     ("bank-file", po::value<string>()->default_value(""), "the template bank file name")
     ("snr-threshold", po::value<float>()->default_value(5.5f), "snr threshold for a detection candidate")
     ("perform-chisq2", po::value<int>()->default_value(1), "perform the chisq2 test")
-    ("chisq2-threshold", po::value<float>()->default_value(0.0f), "threshold for chisq test")
+    ("chisq-threshold", po::value<float>()->default_value(0.0f), "threshold for chisq test")
     ("chisq-bins", po::value<int>()->default_value(16), "number of bins for chisq test")
     ("clustering-algo", po::value<string>()->default_value("window"), "clustering algorithm")
     ("coincidence-algo", po::value<string>()->default_value("simple"), "coincidence algorithm")
@@ -40,6 +40,9 @@ GW::CLIFindchirp::CLIFindchirp(int ac, char* av[]) :
     ("dynamic-range-exponent", po::value<int>()->default_value(32), "dynamic range exponent")
     ("segment-length", po::value<int>()->default_value(1048576), "data segment length in samples")
     ("segment-overlap", po::value<int>()->default_value(524288), "data segment overlap in samples")
+    ("sample-rate", po::value<int>()->default_value(4096), "filter data at given sampling rate")
+    ("spectrum-type", po::value<string>()->default_value("median"), "PSD estimation method (median, mean)")
+
     ;
 
   store_notify_opts(ac, av);
index 0898de7..dd420ea 100644 (file)
@@ -48,18 +48,17 @@ int main(int ac, char* av[]) {
     string     strain_outfile           = myCli.get_write_strain();
     string     psd_outfile              = myCli.get_write_psd();
     string     snr_outfile              = myCli.get_write_snr_series();
-    string     psd_method               = "average";
     double     start_time               = myCli.get_gps_start_time();
     double     end_time                        = myCli.get_gps_end_time();
     double     data_length                     = end_time - start_time;
     cl_uint    data_size                       = 0;
     double     segment_length                  = 0;
     cl_uint       segment_size                         = 0;
-    cl_uint    sampling_rate                   = 4096;
+    cl_uint    sampling_rate                   = myCli.get_sample_rate();
     double     longest_template_length         = 0;
     cl_uint    longest_template_idx            = 0;
     cl_uint    number_of_chi2_bands     = myCli.get_chisq_bins();
-    double     chi2_threshold                  = myCli.get_chisq2_threshold();
+    double     chi2_threshold                  = myCli.get_chisq_threshold();
     cl_uint    chi2_bands[32];
     cl_float   snr_threshold                   = myCli.get_snr_threshold();
     double     flow                     = myCli.get_flow();
@@ -167,7 +166,7 @@ int main(int ac, char* av[]) {
         segment_size, 1. / segment_length, &myCli, &myLogger, &myRTE);
 
     // Calculates the PSD
-    myFilter.PSD(strain, psd, psd_method);
+    myFilter.PSD(strain, psd, myCli.get_spectrum_type());
 
     // Truncates the inverse PSD in time domain
     myFilter.SmoothPSD(psd, sampling_rate * 0.5, sampling_rate * 0.125);
index caa974d..351a637 100644 (file)
@@ -22,7 +22,7 @@ GW::CLISpectrum::CLISpectrum(int ac, char* av[]) :
   m_opt_desc.add_options()
     ("frame-cache", po::value<string>()->default_value(""), "path to input frame cache file")
     ("spectrum-file",  po::value<string>()->default_value(""), "the output file name for the psd")
-    ("spectrum-type",  po::value<string>()->default_value(""), "the method to use for psd calculation")
+    ("spectrum-type",  po::value<string>()->default_value("median"), "the method to use for psd calculation")
     ("channel-name",   po::value<string>()->default_value(""), "the name of the channel to read from the cache file")
     ("write-strain",   po::value<string>()->default_value(""), "write the frame file data to ascii file")
     ("dynamic-range-exponent", po::value<uint>()->default_value(20), "the dynamic exponent to apply (power of 10)")
index ac22f9d..a4f3ca7 100644 (file)
@@ -80,6 +80,8 @@ template <typename T> class Filter : public GWToolsBase
        cl_int HannWindow(      Array <T> * myArray
     );
 
+       T MedianPSDBias(unsigned int nn);
+       
        cl_int PSD(TimeSeries <T> * myArray,    //<! the source TimeSeries
                FrequencySeries <T> * myPSD, //<! the dest FrequencySeries for the PSD
                std::string method_name      //<! method to use
@@ -520,6 +522,22 @@ template <typename T> int Filter<T>::SmoothPSD(
 
 };
 
+/* Compute the bias in the median PSD estimation */
+template <typename T> T Filter<T>:: MedianPSDBias(unsigned int nn)
+{
+    const unsigned int nmax = 1000;
+    T ans = 1.;
+    unsigned int n = (nn - 1) / 2;
+
+    if (nn >= nmax)
+        return 0.69314718055994529;    // ln(2)
+    for (int i = 1; i <= n; ++i) {
+        ans -= 1.0 / (2*i);
+        ans += 1.0 / (2*i + 1);
+    }
+    return ans;
+}
+
 /******************************************************************************/
 /*!
  * The Filter::PSD
@@ -537,43 +555,80 @@ template <typename T> int Filter<T>::PSD(TimeSeries <T> * myArray,
 
     // Determining the segment size and the maximum number of segments
     // with 0.5 overlap
-
-       cl_int segment_size       = myPSD->length;
-       cl_int number_of_segments = myArray->length / (segment_size / 2) - 1;
-
-       // check consistency of input and output sampling
-       // FIXME would it be more useful to just reset it to the correct value,
-       // maybe also printing a warning?
-       if ((myPSD->delta_f * myArray->delta_t * segment_size) != 1.) {
-               m_logger->Log(GW::ERROR, __func__, "Incorrect FrequencySeries->delta_f");
-               return -1;
-       }
+    // FIXME get this parameters from application's CLI
+    cl_int segment_size          = myPSD->length;
+    cl_int number_of_segments = myArray->length / (segment_size / 2) - 1;
+
+    // check consistency of input and output sampling
+    // FIXME would it be more useful to just reset it to the correct value,
+    // maybe also printing a warning?
+    if ((myPSD->delta_f * myArray->delta_t * segment_size) != 1.) {
+        m_logger->Log(GW::ERROR, __func__, "Incorrect FrequencySeries->delta_f");
+        return -1;
+    }
 
     m_logger->Log(GW::DEBUG,__func__,
     "PSD is calculated using %d segments", number_of_segments);
 
-       Array <T> tmpArray  (segment_size, m_commif, m_logger, m_runtime);
-       Array <T> tmpArray2 (segment_size, m_commif, m_logger, m_runtime);
+    Array <T> tmpArray  (segment_size, m_commif, m_logger, m_runtime);
+    Array <T> tmpArray2 (segment_size, m_commif, m_logger, m_runtime);
+    std::vector< Array<T> > arrays;
+
+    myFFT->CreatePlan(segment_size);
+    for (unsigned int i = 0; i < number_of_segments; i++ ) {
+        m_logger->Log(GW::INSANE, __func__, "Calculating FFT of segment %d", i);
+        myArray->Copy((cl_int) i * segment_size / 2, (cl_int)0, (cl_int) segment_size, &tmpArray);
+        HannWindow(&tmpArray);
+        myFFT->Forward(&tmpArray, &tmpArray2);
+        tmpArray2.abs();
+        tmpArray2.sqr();
+        arrays.push_back(tmpArray2);
+    }
+    if (method_name == "mean") {
+        // calculate mean of individual FFTs
+        *myPSD = 0;
+        for (unsigned int i = 0; i < number_of_segments; i++ ) {
+            *myPSD = *myPSD + arrays[i];
+        }
+        *myPSD /= number_of_segments;
+    } else if (method_name == "median") {
+        // calculate median of individual FFTs
+        std::vector<double> values;
+        T                   *seg_buf = (T *)malloc(sizeof(T) * segment_size * number_of_segments),
+                            *median_buf = (T *)malloc(sizeof(T) * segment_size);
+
+        for (unsigned int i = 0; i < number_of_segments; i++ ) {
+            arrays[i].Fetch(seg_buf + i * segment_size);
+        }
+        for (unsigned int j = 0; j < segment_size; j++) {
+            values.clear();
+            for (unsigned int i = 0; i < number_of_segments; i++ ) {
+                values.push_back(seg_buf[i * segment_size + j].real());
+            }
+            std::sort(values.begin(), values.end());
+            if (number_of_segments & 1 == 0)
+                median_buf[j] = 0.5 * (values[number_of_segments / 2 - 1] + values[number_of_segments / 2]);
+            else
+                median_buf[j] = values[number_of_segments / 2];
+        }
+        myPSD->Dispatch(median_buf);
+        free(seg_buf);
+        free(median_buf);
+        // correct for median bias in Gaussian noise
+        *myPSD /= MedianPSDBias(number_of_segments);
+    }
 
-    *myPSD = 0;
+    // calculate mean power of Hann window
+    tmpArray = 1.;
+    HannWindow(&tmpArray);
+    tmpArray.abs();
+    tmpArray.sqr();
+    double win_mean_power = tmpArray.sum().real() / segment_size;
 
-       myFFT->CreatePlan(segment_size);
-       for (int i = 0; i < number_of_segments; i++ ) {
-               m_logger->Log(GW::INSANE,__func__,"Calculating FFT of segment %d", i);
-               myArray->Copy((cl_int) i * segment_size / 2, (cl_int)0, (cl_int) segment_size, &tmpArray);
-               HannWindow(&tmpArray);
-               myFFT->Forward(&tmpArray, &tmpArray2);
-               tmpArray2.abs();
-               tmpArray2.sqr();
-               *myPSD = *myPSD + tmpArray2;
-       }
+    // normalize the PSD
+    *myPSD *= myArray->delta_t / segment_size / win_mean_power;
 
-       // Normalizing
-       // Note: normalizing for two-sided FrequencySeries
-       // FIXME calculate the Hann average power instead of hardcoding it
-       *myPSD *= 1. * myArray->delta_t / segment_size / 0.375 / number_of_segments;
-       
-       return CL_SUCCESS;
+    return CL_SUCCESS;
 };