diff --git a/cmddefs.cpp b/cmddefs.cpp index df4d7b5d..93bf7839 100644 --- a/cmddefs.cpp +++ b/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) // diff --git a/eval.cpp b/eval.cpp index 27932417..54ddd145 100644 --- a/eval.cpp +++ b/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'" ); } diff --git a/fftw/fftwrap.cpp b/fftw/fftwrap.cpp index 34d58873..381bfc0b 100644 --- a/fftw/fftwrap.cpp +++ b/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 & 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 * 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= 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; } } diff --git a/fftw/fftwrap.h b/fftw/fftwrap.h index 64d7b3cb..fb91f2e8 100644 --- a/fftw/fftwrap.h +++ b/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; + }; diff --git a/spindles/spindles.cpp b/spindles/spindles.cpp index e17cb6ea..03923d2c 100644 --- a/spindles/spindles.cpp +++ b/spindles/spindles.cpp @@ -223,6 +223,10 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param ) // (optional) upper bound for core amplitude threshold, e.g. 4.5 < x < 10 if th-max=10 double maximal_threshold = param.has( "th-max" ) ? param.requires_dbl( "th-max" ) : -9 ; + + // winsorize cwt-amplitude + const bool winsorize = param.has( "winsor" ); + const double winsorize_th = winsorize ? param.requires_dbl( "winsor" ) : 0; // minimum spindle core duration (relates to 'th') const double min0_dur_sec = param.has( "min0" ) ? param.requires_dbl( "min0" ) : 0.3; @@ -403,6 +407,12 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param ) const bool cache_peaks_sec = param.has( "cache-peaks-sec" ); const std::string cache_peaks_sec_name = cache_peaks_sec ? param.value( "cache-peaks-sec" ) : ""; + + // + // Spindle-level QC filters + // + + spindle_qc_t q( param ); // // Spindle propagation @@ -623,11 +633,10 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param ) // // Run CWT // - + CWT cwt; cwt.set_sampling_rate( Fs[s] ); - for (int fi=0;fi baseline_fft; - do_fft( d , Fs[s] , &baseline_fft ); + // + // Baseline CWT on signal + // + + q.init( d , Fs[s] ); // @@ -788,10 +795,26 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param ) // Get results for this F_C // - const std::vector & results = cwt.results(fi); - + std::vector results = cwt.results(fi); + + + // + // Add to Q-filter (before any winsorizing or smoothing, etc) + // + q.init_spindle( results ); + + + // + // Optionally winsorize + // + if ( winsorize ) + { + logger << " winsorizing at " << winsorize_th << "\n"; + MiscMath::winsorize( &results , winsorize_th ); + } + // // Get a moving average of the result, 0.1 windows, and get mean // @@ -928,36 +951,7 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param ) threshold_max[p] = maximal_threshold * reaveraged[p]; } - - - // - // Verbose signal display with thresholds - // - - if ( show_cwt_coeff ) - { - - writer.var( "RAWCWT" , "Raw CWT coefficient" ); - writer.var( "CWT" , "CWT coefficient" ); - writer.var( "CWT_TH" , "CWT primary threshold" ); - writer.var( "CWT_TH2" , "CWT secondary threshold" ); - writer.var( "CWT_THMAX" , "CWT maximum threshold" ); - - int np = cwt.points(); - int nf = cwt.freqs(); - if ( np != np0 ) Helper::halt( "internal problem in cwt()" ); - - for (int ti=0;ti averaged_corr = averaged; - for ( int i=0;i locked; - + if ( characterize && some_data ) { - + characterize_spindles( edf , param , signals(s) , bandpass_filtered_status , frq[fi] , @@ -1322,22 +1316,104 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param ) d , // original EEG signal &spindles , // this will be annotated/reduced ( hms ? &starttime : NULL) , - &baseline_fft, + &q, tlocking ? &locked : NULL , // mean signal around spindle troughs p_in_pos_hw , p_in_neg_hw , p_in_pos_slope , p_in_neg_slope - ); - + ); } - + // // Get mean spindle parameters for this channel/frequency // - + if ( characterize ) spindle_stats( spindles , means ); + // + // Verbose signal display with thresholds + // + + if ( show_cwt_coeff ) + { + + int np = cwt.points(); + int nf = cwt.freqs(); + if ( np != np0 ) Helper::halt( "internal problem in cwt()" ); + + // add spindle track for initial/putative spindles (pre-merge, QC) + std::vector sputative(np, false); + if ( spindles1.size() != 0 ) + { + + interval_t curr = spindles1[0]; + int si = 0; + int ti = 0; + int ns = spindles1.size(); + while ( ti < np ) + { + // skip to next event? + if ( (*tp)[ti] >= curr.stop ) + { + ++si; + if ( si == ns ) break; + curr = spindles1[si]; + } + + if ( (*tp)[ti] >= curr.start && (*tp)[ti] < curr.stop ) + sputative[ti] = true; + + ++ti; + } + } + + std::vector sfinal(np, false); + if ( spindles.size() != 0 ) + { + + interval_t curr = spindles[0].tp; + int si = 0; + int ti = 0; + int ns = spindles.size(); + while ( ti < np ) + { + // skip to next event? + if ( (*tp)[ti] >= curr.stop ) + { + ++si; + if ( si == ns ) break; + curr = spindles[si].tp; + } + + if ( (*tp)[ti] >= curr.start && (*tp)[ti] < curr.stop ) + sfinal[ti] = true; + + ++ti; + } + } + + + for (int ti=0;ti 0 ) + if ( characterize && tlocking && spindles.size() > 0 ) { double tlock_min = locked.begin()->first; - double tlock_max = (--locked.end())->first; - - writer.var( "TLOCK" , "Average EEG amplitude time-locked to spindle peak" ); - + double tlock_max = (--locked.end())->first; + std::map::const_iterator ll = locked.begin(); while ( ll != locked.end() ) { @@ -2348,10 +2425,8 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param ) ++ll; } writer.unlevel( "MSEC" ); - } - // // Estiamte of spindle density to console // @@ -2499,6 +2574,8 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param ) avg.so_nearest /= (double)sp_epoch ; writer.value( "AMP" , avg.amp ); + writer.value( "ACT_MX" , avg.norm_amp_max ); + writer.value( "ACT_MN" , avg.norm_amp_mean ); writer.value( "DUR" , avg.dur ); writer.value( "FRQ" , avg.frq ); writer.value( "NOSC" , avg.nosc ); @@ -2520,6 +2597,8 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param ) const int e = edf.timeline.display_epoch( epoch ) - 1; qd.add( writer.faclvl_notime() , "AMP" , e , avg.amp ); + qd.add( writer.faclvl_notime() , "ACT_MX" , e , avg.norm_amp_max ); + qd.add( writer.faclvl_notime() , "ACT_MN" , e , avg.norm_amp_mean ); qd.add( writer.faclvl_notime() , "DUR" , e , avg.dur ); qd.add( writer.faclvl_notime() , "FRQ" , e , avg.frq ); qd.add( writer.faclvl_notime() , "NOSC" , e , avg.nosc ); @@ -2595,6 +2674,7 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param ) writer.var( "MINS" , "Number of minutes for spindle detection" ); writer.var( "DENS" , "Spindle density (per minute)" ); writer.var( "AMP" , "Mean spindle amplitude" ); + writer.var( "AMP" , "Mean spindle amplitude" ); writer.var( "DUR" , "Mean spindle duration" ); writer.var( "FWHM" , "Mean spindle FWHM" ); writer.var( "NOSC" , "Mean spindle number of oscillations" ); @@ -2607,6 +2687,10 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param ) writer.value( "N01" , nspindles_premerge ); // original writer.value( "N02" , nspindles_postmerge ); // post merging writer.value( "N" , (int)spindles.size() ) ; // post merging and QC + + writer.value( "P01" , (int)spindles.size() / (double)nspindles_premerge ); // relative final vs orig (impact of merge) + writer.value( "P02" , (int)spindles.size() / (double)nspindles_postmerge ); // reltive final vs post-merge (impact of Q) + writer.value( "MINS" , t_minutes ); writer.value( "DENS" , spindles.size() / t_minutes ); @@ -2616,6 +2700,8 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param ) writer.value( "Q" , means[ "Q" ] ); writer.value( "AMP" , means["AMP"] ); + writer.value( "ACT_MX" , means["ACT_MX"] ); + writer.value( "ACT_MN" , means["ACT_MN"] ); writer.value( "DUR" , means["DUR"] ); writer.value( "FWHM" , means["FWHM"] ); writer.value( "NOSC" , means["NOSC"] ); @@ -2995,8 +3081,8 @@ void characterize_spindles( edf_t & edf , const std::vector * original_signal , std::vector * spindles , clocktime_t * starttime , - std::map * baseline , - std::map * locked , + spindle_qc_t * q , + std::map * locked , std::vector * p_in_pos_hw , std::vector * p_in_neg_hw , std::vector * p_in_pos_slope , @@ -3106,46 +3192,15 @@ void characterize_spindles( edf_t & edf , } - const int s = edf.header.signal( new_label ); + + // + // Output + // + + const int n = spindles->size(); - // - // Output - // - - const int n = spindles->size(); - - // - // Spindle-level QC filters (set default at 0, i.e spindle-activity must be more likely) - // - - bool qc_q = param.has( "noq" ) ? false : true; - double qc_qmin = 0 , qc_qmax = -1; - if ( qc_q ) - { - if ( param.has( "q" ) ) { qc_q = true; qc_qmin = param.requires_dbl( "q" ); } - if ( param.has( "q-max" ) ) { qc_q = true; qc_qmax = param.requires_dbl( "q-max" ); } - if ( qc_qmin < 0 ) qc_q = false; - } - - if ( qc_q && ( target_f < 10 || target_f > 16 ) ) - { - Helper::halt( "spindle fc is outside of 10-16 Hz range, but 'q' is set\n add q=-9 to remove freq-based QC filter" ); - } - - // if ( 1 ) - // { - // std::cout << " INT = " << edf.timeline.wholetrace() << "\n"; - // slice_t slice( edf , s , edf.timeline.wholetrace() ); - // const std::vector * d = slice.pdata(); - // std::cout << "N = " << d->size() << "\n"; - // std::cout << std::setprecision(12); - // for ( int ii=0; ii < d->size(); ii++ ) - // { - // std::cout << "FXX\t" << ii << "\t" << (*d)[ii] << "\n"; - // } - // } - + // // track if we QC any spindles out at this step // @@ -3179,7 +3234,7 @@ void characterize_spindles( edf_t & edf , slice_t slice( edf , s , spindle->tp ); // nb. we are going to mean-centre these data too - // for this particular spindle... avoids someo weird cases + // for this particular spindle... avoids some weird cases // where we have huge offsets in original signal amplitudes around // this spindle, meaning that even after BP filter, we can // have spindle intervals with no zero-crossings @@ -3202,6 +3257,7 @@ void characterize_spindles( edf_t & edf , // spindle->isa = 0; + spindle->norm_amp_max = 0; if ( averaged != NULL ) { @@ -3210,10 +3266,22 @@ void characterize_spindles( edf_t & edf , int stop = spindle->stop_sp; for (int s=start;s<=stop;s++) - spindle->isa += (*averaged)[s] ; + { + spindle->isa += (*averaged)[s] ; + + // track max (ACT_MX) + if ( (*averaged)[s] > spindle->norm_amp_max ) + spindle->norm_amp_max = (*averaged)[s]; + } + + // ACT_MN (i.e. simple mean, no duration info) + spindle->norm_amp_mean = spindle->isa / (double)(stop-start+1); + + // ISA_S (i.e. includes duration info as norm by Fs only + spindle->isa /= (double)Fs; + } - - spindle->isa /= (double)Fs; + // @@ -4029,89 +4097,27 @@ void characterize_spindles( edf_t & edf , spindle->frq = spindle->nosc / (double)spindle->dur; // - // FFT on original data, compared to baseline + // Q stats? // - - - if ( baseline ) - { - - slice_t slice0( edf , s0 , spindle->tp ); - - // copy over for ranges - std::map spindle_fft = *baseline; - - // fixed at: - // 0.5..4 - // 4..8 - - // 10..13.5 - // 13.5..16 - - // 20..30 - - do_fft( slice0.nonconst_pdata() , Fs , &spindle_fft ); - - // calculate enrichment (log10-scale), so set min to v. low... - double q_spindle = -999 , q_baseline = -999; - - std::map::const_iterator ff = spindle_fft.begin(); - while ( ff != spindle_fft.end() ) - { - - const double & baseline_band_power = (*baseline)[ ff->first ] ; - const double & spindle_band_power = ff->second; - const freq_range_t & band = ff->first; - - // relative enrichment (to baseline) [ log scale ] - double re = spindle_band_power - baseline_band_power; - - // relative enrichment (to baseline) - //double re = spindle_band_power / baseline_band_power; - - // store - spindle->enrich[ ff->first ] = re; - - // calculate overall q score - // take 'spindle' as the two middle categories - - // quality score: 10..16 is spindle range - // - if ( band.first <= 16 && band.second >= 10 ) - { - // i.e. get largest of slow and fast bands - if ( re > q_spindle ) q_spindle = re; - } - else - { - // i.e. get largest of non-spindle bands - if ( re > q_baseline) q_baseline = re; - } - ++ff; - } + + if ( q->qc_q ) + { - // relative relative enrichment [ log scale ] - spindle->qual = q_spindle - q_baseline ; - - // relative relative enrichment [ log scale ] - //spindle->qual = q_spindle / q_baseline ; - - // i.e. max( B_S / B_S_0 ) / max( B_NS / B_NS_0 ) - // or logs + // get Q-score + spindle->qual = q->qual( spindle ); - // QUAL filter? + // failed? + if ( spindle->qual < q->qc_qmin ) spindle->include = false; + if ( q->qc_qmax > 0 && spindle->qual > q->qc_qmax ) spindle->include = false; - if ( qc_q ) - { - if ( spindle->qual < qc_qmin ) spindle->include = false; - if ( qc_qmax > 0 && spindle->qual > qc_qmax ) spindle->include = false; - } - + if ( q->verbose_all || ( q->verbose_failed && ! spindle->include ) ) + q->output( spindle , i ); } - + + // - // [ REMOVE ] Optional, time-locked analysis? [ for QC+ spindles only ] + // [ REMOVE - i.e. do generically ] Optional, time-locked analysis? [ for QC+ spindles only ] // if ( locked && spindle->include ) @@ -4123,10 +4129,7 @@ void characterize_spindles( edf_t & edf , const double window_left = - ( window_sec / 2.0 ); // start point (left of window) for peak minus half window - int orig_sp = spindle->start_sp + lowest_idx - ( window_sec / 2.0 ) * Fs ; - - // logger << "os = " <start_sp + lowest_idx - ( window_sec / 2.0 ) * Fs ; uint64_t centre = tp[ lowest_idx ]; interval_t i0; @@ -4148,19 +4151,6 @@ void characterize_spindles( edf_t & edf , const int bin = pos / step_tp ; const double fbin = window_left + bin * step_sec ; - // weight by CWT for spindle... - - // if ( 0 ) - // std::cout << "TL\t" - // << edf.id << "\t" - // << target_f << "\t" - // << orig_sp << "\t" - // << i << "\t" - // << l << "\t" - // << (*averaged)[orig_sp] << "\t" - // << (*d0)[l] << "\t" - // << (*averaged)[orig_sp] * (*d0)[l] << "\n"; - ++orig_sp; (*locked)[ fbin ] += (*d0)[l]; @@ -4192,9 +4182,10 @@ void characterize_spindles( edf_t & edf , spindles->clear(); for (int i=0;ipush_back( copy_spindles[i] ); - logger << " QC'ed spindle list from " << copy_spindles.size() << " to " << spindles->size() << "\n"; + logger << " Q-score pruned spindle list from " << copy_spindles.size() << " to " << spindles->size() << "\n"; } - + else + logger << " no spindles removed based on Q-scores\n"; // @@ -4211,7 +4202,6 @@ void characterize_spindles( edf_t & edf , } } - // // Remove tmp channel we created // @@ -4221,20 +4211,16 @@ void characterize_spindles( edf_t & edf , int s = edf.header.signal( new_label ); edf.drop_signal( s ); } - - - + } void per_spindle_output( std::vector * spindles , param_t & param , - clocktime_t * starttime , - std::map * baseline ) + clocktime_t * starttime , + spindle_qc_t & q ) { - const bool enrich_output = param.has( "enrich" ); - const int n = spindles->size(); const bool has_coupling = param.has( "so" ); @@ -4284,6 +4270,8 @@ void per_spindle_output( std::vector * spindles , } writer.value( "AMP" , spindle->amp ); + writer.value( "ACT_MX" , spindle->norm_amp_max ); + writer.value( "ACT_MN" , spindle->norm_amp_mean ); writer.value( "DUR" , spindle->dur ); writer.value( "FWHM" , spindle->fwhm ); writer.value( "NOSC" , spindle->nosc ); @@ -4357,30 +4345,13 @@ void per_spindle_output( std::vector * spindles , if ( has_if ) writer.value( "IF" , spindle->if_spindle ); - - // - // Enrichment relative to the baseline - // - - if ( baseline ) - { - + // Q-score + if ( q.qc_q ) + { writer.value( "Q" , spindle->qual ); writer.value( "PASS" , spindle->include ); - - if ( enrich_output ) - { - std::map::const_iterator bb = spindle->enrich.begin(); - while ( bb != spindle->enrich.end() ) - { - writer.level( globals::print( bb->first ) , globals::band_strat ); - writer.value( "ENRICH" , bb->second ); - ++bb; - } - writer.unlevel( globals::band_strat ); - } } - + } // end of per-spindle output @@ -4390,69 +4361,160 @@ void per_spindle_output( std::vector * spindles , - - - -void do_fft( const std::vector * d , const int Fs , std::map * freqs ) +void spindle_qc_t::proc( param_t & param ) { - - // Fixed parameters:: use 4-sec segments with 2-second - // overlaps and Hanning window - double overlap_sec = 2; - double segment_sec = 4; - double length_sec = d->size() / (double)Fs; + if ( param.has( "noq" ) ) + { + qc_q = false; + return; + } - // check length - if ( length_sec <= ( segment_sec + overlap_sec ) ) + qc_q = true; + + // thresholds + + if ( param.has( "q" ) ) + qc_qmin = param.requires_dbl( "q" ); + else + qc_qmin = 0; + + if ( qc_qmin <= -9 ) { - overlap_sec = 0; - segment_sec = length_sec; + qc_q = false; + return; } - const int total_points = d->size(); - const int segment_points = segment_sec * Fs; - const int noverlap_points = overlap_sec * Fs; - - int noverlap_segments = floor( ( total_points - noverlap_points) - / (double)( segment_points - noverlap_points ) ); + if ( param.has( "q-max" ) ) + qc_qmax = param.requires_dbl( "q-max" ); + + // outputs + verbose_failed = param.has( "q-verbose" ); + verbose_all = param.has( "q-verbose-all" ); + + // defaults : 6 Hz (5 cycles) + // : 30 Hz (7 cycles) + // OR?? : 35 Hz ( 7 cycles) + + ofc.clear(); + ofc.push_back( 6 ); + ofc.push_back( 35 ); + ocycles.clear(); + ocycles.push_back( 5 ); + ocycles.push_back( 7 ); - // std::cout << "total_points " << total_points << " " << segment_points << " " << noverlap_points << " " << noverlap_segments << "\n"; + // target fcs + if ( param.has( "q-frq" ) ) + { + ofc.clear(); + ocycles.clear(); + std::vector pp = param.dblvector( "q-frq" ); + if ( pp.size() % 2 != 0 ) + Helper::halt( "expecting q-frq=fc,cycles(,fc,cycles,...}" ); + for (int i=0; iclear(); +} + +void spindle_qc_t::init( const std::vector * _ts , int _Fs ) +{ + + // keep a copy of the original time-series + // (this is *assumed* to still be in scope whenever q.output() is called) - // 1 .. 25 - // for (double f=0.5 ; f <= 25.5 ; f++ ) - // (*freqs)[ freq_range_t( f , f+1 ) ] = 0; + ts = _ts; + Fs = _Fs; - (*freqs)[ freq_range_t( 0.5 , 4 ) ] = 0 ; - (*freqs)[ freq_range_t( 4 , 8 ) ] = 0 ; - (*freqs)[ freq_range_t( 10 , 13.5 ) ] = 0 ; - (*freqs)[ freq_range_t( 13.5 , 16 ) ] = 0 ; - (*freqs)[ freq_range_t( 20 , 30 ) ] = 0 ; + ocwt.resize( nob ); - // populate - pwelch.psdmean( freqs ); + CWT cwt; + cwt.set_sampling_rate( Fs ); + for (int fi=0; fi::iterator ff = freqs->begin(); - while ( ff != freqs->end() ) + for (int fi=0; fisecond = log10( ff->second ); - //ff->second = ff->second ; - ++ff; + ocwt[fi] = cwt.results(fi); + const double om = MiscMath::mean( ocwt[fi] ); + for (int i=0; i & s ) +{ + scwt = s; + double sm = MiscMath::mean( scwt ); + for (int i=0; i o( nob , 0 ); + const int np = spindle->stop_sp - spindle->start_sp + 1; + + for (int i=spindle->start_sp; i<= spindle->stop_sp; i++) + { + s += scwt[i]; + for (int j=0; j omx ) omx = o[j]; + if ( omx <= 0 ) omx = 1e-12; + + return log2( s ) - log2( omx ); } +void spindle_qc_t::output( spindle_t * spindle , const int n ) +{ + // verbose outputs + const int np = spindle->stop_sp - spindle->start_sp + 1; + + int s1 = spindle->start_sp - 2 * Fs; + int s2 = spindle->stop_sp + 2 * Fs; + if ( s1 < 0 ) s1 = 0; + if ( s2 >= ts->size() ) s2 = ts->size() - 1 ; + + for (int i=s1; i<=s2; i++) + { + const bool in_spindle = i >= spindle->start_sp && i <= spindle->stop_sp; + + std::cout << "QS\t" + << n << "\t" + << in_spindle << "\t" + << (*ts)[i] << "\t" + << scwt[i] ; + for (int j=0; j & spindles , std::map & results ) @@ -4460,6 +4522,8 @@ void spindle_stats( const std::vector & spindles , std::map & spindles , std::map enrich; // versus baseline - std::vector::const_iterator ii = spindles.begin(); while ( ii != spindles.end() ) { @@ -4538,25 +4600,24 @@ void spindle_stats( const std::vector & spindles , std::mapisa; posisa += ii->posisa; negisa += ii->negisa; + + norm_amp_max += ii->norm_amp_max; + norm_amp_mean += ii->norm_amp_mean; possp += ii->possp; negsp += ii->negsp; qual += ii->qual; - // relative enrichment compared to baseline - std::map::const_iterator ss = ii->enrich.begin(); - while ( ss != ii->enrich.end() ) - { - enrich[ ss->first ] += ss->second; - ++ss; - } - ++ii; } results[ "AMP" ] = amp /(double)denom; + + results[ "ACT_MX" ] = norm_amp_max /(double)denom; + results[ "ACT_MN" ] = norm_amp_mean /(double)denom; + results[ "TOTDUR" ] = dur; results[ "DUR" ] = dur / (double)denom; results[ "FWHM" ] = fwhm / (double)denom; @@ -4608,307 +4669,7 @@ void spindle_stats( const std::vector & spindles , std::map::iterator ee = enrich.begin(); - while ( ee != enrich.end() ) - { - results[ "E" + globals::print( ee->first ) ] = ee->second / (double)denom; - ++ee; - } - - -} - - - -annot_t * spindle_bandpass( edf_t & edf , param_t & param ) -{ - - - // - // Attach signals - // - - std::string signal_label = param.requires( "sig" ); - - signal_list_t signals = edf.header.signal_list( signal_label ); - - const int ns = signals.size(); - - // - // Obtain sampling freqs (Hz) - // - - std::vector Fs = edf.header.sampling_freq( signals ); - - - // - // Annotations to save - // - - annot_t * a = edf.timeline.annotations.add( "spindles-v2" ); - a->description = "Martin et al. spindles" ; - - - // - // For each signal - // - - for (int s = 0 ; s < ns ; s++ ) - { - - if ( edf.header.is_annotation_channel( signals(s) ) ) continue; - - // - // Based on: Martin et al. "Topography of age-related changes in - // sleep spindles", Neurobio Aging 34(2), 2013, pp 468-476 - // - // Method 'A4' from Warby et al. - // - - // [# Band-pass filter EEG, calculate RMS in sliding windows and apply a - // constant threshold. Detect a spindle if the RMS exceeds a constant - // threshold for 0.3-3 s.] - - // Parameters - - const double p_resolution = 0.25; - const double p_percentile = 95; - const double p_window_length = Fs[s] * p_resolution; - - // 1. Bandpass filter signal from C3-M2 in the 11-15 Hz band - // 2. Calculate the RMS of the bandpass filtered signal with a time - // resolution of 25 ms using a time window of 25 ms [# no overlap] - // 3. threshold <- 95th percentile of RMS signal [# only S2+S3+S4] - // 4. if ( RMS > threshold && 0.3s <= duration above threshold <= 3s ) - // then [Detect spindle] - - - // - // Filter entire signal - // - - std::vector fripple( 1, 0.01 ); - std::vector ftw( 1, 0.5 ); - - // ripple = 0.005 , transition width (Hz) = 0.5 Hz - dsptools::apply_fir( edf , signals(s) , fir_t::BAND_PASS , - 1 , // Kaiser window - fripple, ftw, // ripple , TW - 10 , 16 - ); - - // - // Get windows of 0.25seconds, no overlap (i.e. advance by 0.25) - // - - int ne = edf.timeline.set_epoch( p_resolution , p_resolution ); - - // - // Aggregate RMS per window - // - - std::vector rms; - - - // - // Get data - // - - while ( 1 ) - { - - int epoch = edf.timeline.next_epoch(); - - if ( epoch == -1 ) break; - - interval_t interval = edf.timeline.epoch( epoch ); - - // - // Get data - // - - slice_t slice( edf , signals(s) , interval ); - - const std::vector * signal = slice.pdata(); - - // - // Calculate RMS for window - // - - double t = MiscMath::rms( *signal ); - - rms.push_back(t); - - } // next 0.25 window - - - // - // Get threshold (95th percentile) - // - - const int n_bins = rms.size(); - - const int t95 = n_bins * ( p_percentile / 100.0 ); - - const double threshold = MiscMath::kth_smallest_preserve( rms , t95 ); - - const uint64_t bin_ms = p_resolution * globals::tp_1sec; - - std::vector spindles; - - uint64_t start = 0 , stop = 0; - - // count of current spindle 'length' (in 0.25s windows - int scnt = 0; - - // for 0.3 to 3s duration, means at least 2 bins, - // but not more than 12 - - for (int i=0;i= threshold ) - { - if ( scnt == 0 ) // start of putative spindle - { - scnt = 1; - start = i * bin_ms; - } - else // continue a window - { - ++scnt; - stop = (i+1) * bin_ms - 1; - } - } - else - { - - if ( scnt ) // close a window? - { - if ( scnt >= 2 && scnt <= 12 ) - { - spindles.push_back( spindle_t( start , stop , 0 , 0 ) ); - } - } - scnt = 0; // reset scnt in any case - } - } - // - // Characterisation of each spindle - // - - bool bandpass_filtered = true; - - characterize_spindles( edf , param , signals(s) , bandpass_filtered , 13 , - "bandpass" , - NULL , NULL , &spindles , NULL , NULL , NULL , NULL , NULL , NULL ); - - std::map means; - - spindle_stats( spindles , means ); - - // - // Save in the annotation class - // - - std::vector::const_iterator ii = spindles.begin(); - while ( ii != spindles.end() ) - { - const interval_t & spindle = ii->tp; - instance_t * instance = a->add( signals.label(s) , spindle , signals.label(s) ); - - //instance->set( "mid", (double)(tp_mid * globals::tp_duration ) ); - - ++ii; - } - - - // - // Per-spindle level output - // - - if ( false ) - { - - std::vector::const_iterator ii = spindles.begin(); - int cnt = 0; - - while ( ii != spindles.end() ) - { - - const spindle_t & spindle = *ii; - - writer.level( ++cnt , "SPINDLE" ); - writer.var( "SINGLE_SP_START" , "Single spindle start time-point" ); - writer.var( "SINGLE_SP_STOP" , "Single spindle stop time-point" ); - writer.var( "SINGLE_SP_DUR" , "Single spindle stop time-point" ); - - writer.value( "SINGLE_SP_START" , spindle.tp.start * globals::tp_duration ); - writer.value( "SINGLE_SP_STOP" , spindle.tp.stop * globals::tp_duration ); - writer.value( "SINGLE_SP_DUR" , (spindle.tp.stop-spindle.tp.start+1)/(double)globals::tp_1sec ); - - ++ii; - } - writer.unlevel( "SPINDLE" ); - } - - const double t_minutes = ( n_bins * p_resolution ) / 60.0 ; - - bool empty = spindles.size() == 0; - - if ( empty ) - std::cout << "INDIV" << "\t" - << edf.id << "\t" - << "[" << globals::current_tag << "]\t" - << signals.label(s) << "\t" - << 0 << "\t" - << t_minutes << "\t" - << 0 << "\t" - << 0 << "\t" - << "NA\t" - << "NA\t" - << "NA\t" - << "NA\t" - << "NA\t" - << "NA\n"; - else - std::cout << "INDIV" << "\t" - << edf.id << "\t" - << "[" << globals::current_tag << "]\t" - << signals.label(s) << "\t" - << spindles.size() << "\t" - << t_minutes << "\t" - << spindles.size() / t_minutes << "\t" - << means["TOTDUR"] << "\t" - << means["AMP"] << "\t" - << means["DUR"] << "\t" - << means["NOSC"] << "\t" - << means["FRQ"] << "\t" - << means["FFT"] << "\t" -// << means["TREND"] << "\t" -// << means["ABSTREND"] << "\t" - << means["SYMM"] << "\t" - << means["SYMM2"] << "\n"; - - - } // Next signal - - - // a->save( "spindle.annot" ); - - return a; - } -// helper function -void write_if_exists( const std::string & s , const std::map & means ) -{ - std::map::const_iterator ss = means.find( s ); - if ( ss != means.end() ) writer.value( s , ss->second ); -} - - diff --git a/spindles/spindles.h b/spindles/spindles.h index 85a337e1..a6967db9 100644 --- a/spindles/spindles.h +++ b/spindles/spindles.h @@ -51,6 +51,7 @@ struct spindle_t // spindle properties double amp, dur, fwhm, nosc, frq, fft, symm, symm2, isa; + double norm_amp_max, norm_amp_mean; double chirp, frq_h1, frq_h2, frq_range; // neg/pos defined by HW @@ -95,19 +96,68 @@ struct spindle_t }; -annot_t * spindle_bandpass( edf_t & , param_t & ); -annot_t * spindle_wavelet( edf_t & , param_t & ); +struct spindle_qc_t { + spindle_qc_t( param_t & param ) + { -// helper function for FFT -void do_fft( const std::vector * d , const int Fs , std::map * fft ); + qc_qmin = 0; + qc_qmax = -1; + + // cwt + ofc.clear(); + + // get params + proc( param ); + + } + + bool qc_q; + + double qc_qmin, qc_qmax; + + // original time-series (for verbose outputs) + const std::vector * ts; + int Fs; + + // CWT : non-spindle bands + int nob; // number of other-Fcs [ ofc.size() ] + std::vector ofc; + std::vector ocycles; + std::vector > ocwt; + + // spindle stats + std::vector scwt; // varies for each Fc + + // verbose output + bool verbose_failed; + bool verbose_all; + + // functions + void proc( param_t & ); + void init( const std::vector * ts , int Fs ); + void init_spindle( const std::vector & scwt ); + double qual( spindle_t * spindle ); + + // verbose output + void output( spindle_t * spindle , const int n ); + +}; + + +annot_t * spindle_wavelet( edf_t & , param_t & ); + +// helper function for FFT +void do_fft( const std::vector * d , const int Fs , spindle_qc_t * q , bool baseline ); +void do_fft( const std::vector * d , const int Fs , const int np, spindle_qc_t * q , bool baseline ); + // helper to get spindle stats void spindle_stats( const std::vector & spindles , std::map & ) ; // helper -void write_if_exists( const std::string & s , const std::map & means ) ; +//void write_if_exists( const std::string & s , const std::map & means ) ; void characterize_spindles( edf_t & edf , param_t & param , @@ -115,11 +165,11 @@ void characterize_spindles( edf_t & edf , bool bandpass_filtered , const double target_f, const std::string & label, - const std::vector * averaged , + const std::vector * averaged , const std::vector * original_signal , std::vector * spindles , - clocktime_t * starttime , - std::map * baseline = NULL , + clocktime_t * starttime , + spindle_qc_t * q = NULL , std::map * locked = NULL , std::vector * in_pos_hw = NULL , std::vector * in_neg_hw = NULL , @@ -130,13 +180,8 @@ void characterize_spindles( edf_t & edf , void per_spindle_output( std::vector * spindles , param_t & param , - clocktime_t * starttime , - std::map * baseline ); - - - - - + clocktime_t * starttime , + spindle_qc_t & q ); //