reimpl of q-filtering

cmddefs.cpp

@@ -2141,6 +2141,8 @@ void cmddefs_t::init()
   add_table( "SPINDLES" , "CH,F" , "Individual-level output" );
   add_var( "SPINDLES" , "CH,F" , "DENS" , "Spindle density (count per minute)" );
   add_var( "SPINDLES" , "CH,F" , "AMP" , "Mean spindle amplitude (uV or mV units)" );
+  add_var( "SPINDLES" , "CH,F" , "ACT_MX" , "Mean max spindle activity (normed CWT)" );
+  add_var( "SPINDLES" , "CH,F" , "ACT_MN" , "Mean average spindle activity (normed CWT)" );
   add_var( "SPINDLES" , "CH,F" , "DUR" , "Mean spindle duration (core+flanking region)" );
   add_var( "SPINDLES" , "CH,F" , "NOSC" , "Mean number of oscillations per spindle" );
   add_var( "SPINDLES" , "CH,F" , "FWHM" , "Mean spindle FWHM (full width at half maximum)" );
@@ -2172,6 +2174,8 @@ void cmddefs_t::init()
 
   add_table( "SPINDLES" , "CH,F,SPINDLE" , "Spindle-level output [per-spindle]" ); 
   add_var( "SPINDLES" , "CH,F,SPINDLE" , "AMP" , "Spindle amplitude (uV or mV units)" );
+  add_var( "SPINDLES" , "CH,F,SPINDLE" , "ACT_MX" , "Max spindle activity (normed CWT)" );
+  add_var( "SPINDLES" , "CH,F,SPINDLE" , "ACT_MN" , "Average spindle activity (normed CWT)" );
   add_var( "SPINDLES" , "CH,F,SPINDLE" , "CHIRP" , "Spindle chirp (-1 to +1)" );
   add_var( "SPINDLES" , "CH,F,SPINDLE" , "DUR" , "Spindle duration (seconds)" );
   add_var( "SPINDLES" , "CH,F,SPINDLE" , "FWHM" , "Spindle FWHM (seconds)" );
@@ -2192,6 +2196,9 @@ void cmddefs_t::init()
   add_var( "SPINDLES" , "CH,F,SPINDLE" , "SYMM" , "Symmetry index (relative position of peak)" );
   add_var( "SPINDLES" , "CH,F,SPINDLE" , "SYMM2" , "Folded symmetry index (0=symmetrical, 1=asymmetrical)" );
   hide_var( "SPINDLES" , "CH,F,SPINDLE" , "IF" , "Mean frequency per spindle over duration [if]" );
+
+  add_table( "SPINDLES" , "CH,F,B,SPINDLE" , "Band enrichment (per-spindle)" );
+  add_var( "SPINDLES" , "CH,F,B,SPINDLE" , "ENRICH" , "Spindle enrichment" );
 
   hide_table( "SPINDLES" , "CH,F,RELLOC" , "Mean IF stratified by relative location in spindle [if]" );
   hide_var( "SPINDLES" , "CH,F,RELLOC" , "IF" , "Mean frequency of all spindles, per relative position within the spindle (five bins)" );
@@ -2274,14 +2281,18 @@ void cmddefs_t::init()
   // show-coef verbose output
 
   add_table( "SPINDLES" , "F,CH,T" , "Verbose threshold/coefficient output [show-coeff]" );
+  add_var( "SPINDLES" , "F,CH,T" , "SEC" , "Time (sec)" );
   add_var( "SPINDLES" , "F,CH,T" , "RAWCWT" , "Raw CWT coefficient" );
   add_var( "SPINDLES" , "F,CH,T" ,"CWT" , "CWT coefficient" );
+  add_var( "SPINDLES" , "F,CH,T" ,"AVG" , "Averaged CWT coefficient" );
+  add_var( "SPINDLES" , "F,CH,T" ,"AVG_CORR" , "Averaged baseline-corrected CWT coefficient" );
   add_var( "SPINDLES" , "F,CH,T" ,"CWT_TH" , "CWT primary threshold" );
   add_var( "SPINDLES" , "F,CH,T" ,"CWT_TH2" , "CWT secondary threshold" );
   add_var( "SPINDLES" , "F,CH,T" ,"CWT_THMAX" , "CWT maximum threshold" );
+  add_var( "SPINDLES" , "F,CH,T" ,"PUTATIVE" , "Pre-QC spindle" );
+  add_var( "SPINDLES" , "F,CH,T" ,"SPINDLE" , "Post-QC spindle" );
 
 
-
   //
   // SO (duplicated from SO command below)
   //

eval.cpp

@@ -2686,9 +2686,8 @@ void proc_spindles( edf_t & edf , param_t & param )
   std::string method = param.has( "method" ) ? param.value( "method" ) : "wavelet" ; 
 
   annot_t * a = NULL;
-
-  if      ( method == "bandpass" ) a = spindle_bandpass( edf , param );
-  else if ( method == "wavelet" ) a = spindle_wavelet( edf , param );
+
+  if ( method == "wavelet" ) a = spindle_wavelet( edf , param );
   else Helper::halt( "SPINDLE method not recognized; should be 'bandpass' or 'wavelet'" );
 
 }

fftw/fftwrap.cpp

@@ -105,7 +105,7 @@ void FFT::init( int Ndata_, int Nfft_, int Fs_ , fft_t type_ , window_function_t
   //
 
   w.resize( Ndata , 1 ); // i.e. default of no window
-  
+
   normalisation_factor = 0;  
   if      ( window == WINDOW_TUKEY50 ) w = MiscMath::tukey_window(Ndata,0.5);
   else if ( window == WINDOW_HANN )    w = MiscMath::hann_window(Ndata);
@@ -114,7 +114,8 @@ void FFT::init( int Ndata_, int Nfft_, int Fs_ , fft_t type_ , window_function_t
   for (int i=0;i<Ndata;i++) normalisation_factor += w[i] * w[i];
   normalisation_factor *= Fs;  
   normalisation_factor = 1.0/normalisation_factor;
-  //  std::cout << "normalisation_factor = " << normalisation_factor << "\n";
+  //std::cout << "normalisation_factor = " << normalisation_factor << "\n";
+
 } 
 
 bool FFT::apply( const std::vector<double> & x )
@@ -833,6 +834,9 @@ void PWELCH::process()
   if ( average_adj ) 
     fft0.average_adjacent();
 
+  if ( ! do_normalization )
+    fft0.norm_fac( 1.0 );
+
   psd.resize( fft0.cutoff , 0 );
   N = fft0.cutoff;
 
@@ -1003,17 +1007,20 @@ void PWELCH::psdmean( std::map<freq_range_t,double> * f )
 
       const double & lwr = ii->first.first;
       const double & upr = ii->first.second; 
-
+      //      std::cout << " searching " << lwr << " " << upr << "\n";
+
       // add is l <= x < y 
       double r = 0;
       int c = 0;
       for (int i=0;i<N;i++) 
 	{
 	  if ( freq[i] >= upr ) break;
-	  if ( freq[i] >= lwr ) { ++c; r += psd[i]; }	  
+	  if ( freq[i] >= lwr ) { ++c; r += psd[i]; }
 	}
 
-      ii->second = r / (double) c;
+      if ( c != 0 ) 
+	ii->second = r / (double) c;
+
       ++ii;
     }
 }

fftw/fftwrap.h

@@ -208,7 +208,9 @@ class real_FFT
     }
 
   void init( int Ndata , int Nfft , int Fs , window_function_t window = WINDOW_NONE );
-
+
+  void norm_fac( const double f ) { normalisation_factor = f; } 
+
   void reset() ;
 
   ~real_FFT();
@@ -420,11 +422,12 @@ class PWELCH
 	 bool use_median = false ,
 	 bool calc_seg_sd = false , 
 	 bool average_adj = false ,
-	 bool use_nextpow2 = false ) 
+	 bool use_nextpow2 = false ,
+	 bool do_normalization = true ) 
    : data(data) , Fs(Fs) , M(M) , noverlap_segments(noverlap_segments) , 
      window(W),
      use_median(use_median), calc_seg_sd(calc_seg_sd),
-     average_adj(average_adj) , use_nextpow2( use_nextpow2 )
+     average_adj(average_adj) , use_nextpow2( use_nextpow2)  , do_normalization(do_normalization) 
   {
 
     // calculate implied overlap in actual data-points the above
@@ -534,6 +537,9 @@ class PWELCH
   // always set NFFT to the next power of 2
   bool use_nextpow2;
 
+  // do 1/N^2 norm (incl. window)
+  bool do_normalization;
+
 };