Skip to content

Spindles and slow oscillations (SO)

Command Description
SPINDLES Wavelet-based spindle detection
SO Slow oscillation detection


Wavelet-based sleep spindle detection

This command detects spindles using a wavelet-based approach. Optionally, it also detects slow oscillations (SO) and the temporal coupling between spindles and SO. A target frequency is specified, which corresponds to the center frequency of the Morlet wavelet. A single SPINDLES command can detect spindles at different frequencies, and on different channels; there are also options for collating putative spindles that are coincident in time but observed on different channels and/or at neighbouring frequencies.

The SPINDLES command has quite a large number of options (although most of them are not necessary for its basic operation). For clarity, this section is split into a number of sub-sections that describe different aspects of this command:

After running the wavelet transform, spindles are detected based on the following rules:

  • the magnitude of the wavelet transform is smoothed (default sliding window of 0.1, win) and scaled by the whole-signal mean (or median, if the median flag is used)

  • here, whole-signal means all unmasked epochs available to this command (i.e. may be only all N2 epochs rather than all epochs)

  • a spindle is called when a core interval of at least 0.3 seconds (min0) exhibits a signal at least 4.5 times (th) the mean (or median)

  • furthermore, all spindles must also exhibit an extended interval of at least 0.5 seconds (min) with a signal at least 2 times (th2) the mean (or median)

  • putative spindle intervals cannot be longer than 3 seconds (max); further, spindles within 0.5 seconds (merge) are merged (unless the resulting spindle is greater than max seconds)

  • putative spindles are also filtered based on a quality metric q, the threshold of which is given by the q parameter

Following this heuristic procedure, in addition to the count and density (count per minute) of spindles, the following quantities are estimated to characterize the morphology of each spindle, after bandpass-filtering the signal for FC ± 2 Hz:

  • spindle duration (of core plus flanking interval) in seconds, and number of oscillations

  • spindle amplitude, based on the maximum peak-to-peak amplitude

  • mean and maximum of the transformed wavelet statistic

  • spindle frequency, based on counting zero-crossings, and secondarily based on the modal frequency from a FFT

  • spindle symmetry, based on the relative location of the spindle's central "peak" (the point of maximum peak-to-peak amplitude), and a second "folded" symmetry index, calculated as 2|S-0.5|

  • spindle chirp, based on the log ratio of mean peak-to-peak time intervals in the first versus the second half of the spindle

  • a measure of integrated spindle activity (ISA), based on the sum of normalized wavelet coefficients in the spindle interval

  • a quality measure q based on the ratio of 1) the relative enrichment of power in the sigma band (for the spindle compared to baseline) versus 2) the comparable relative enrichment in power across all bands (i.e. again for the spindle compared to baseline)

Averaging over all spindles in a given run, Luna reports the mean spindle density, amplitude, duration, frequency, number of oscillations, chirp, symmetry indices, etc. Luna also reports per-spindle statistics, as well as per-epoch counts of spindles.

Basic usage


The most basic parameter is fc, which specifies the target frequency (or frequencies) for the wavelet(s). Combined with the cycles parameter, this defines the distribution of spindle frequencies that are targeted by the wavelet. See also the CWT-DESIGN command, which gives further information on wavelet properties.

The primary output variables for an individual are spindle density (DENS), mean duration (DUR) and frequency (FRQ or FFT, both are typically very similar). An additional metric of interest is integrated spindle activity (ISA), which, for a given spindle, is the sum of the normalized wavelet coefficients in the spindle interval. Higher values indicate more evidence in favour of a spindle, and longer and higher amplitude spindles. Averaging across spindles detected, for an individual, ISAS is the mean ISA, which reflects the average amplitude and duration of that individual's spindles; ISAM is the mean ISA per minute, which reflects the combination of spindle amplitude, duration and density; ISAT is the total ISA, which reflects the combination of spindle amplitude, duration and count. Because it is normalized by the individual's baseline, ISAS may be a better metric on which to base inter-individual comparisons, compared to AMP which is more likely to reflect technical differences in recordings.

Primary options
Parameter Example Description
sig sig=C3,F3 Restrict analysis to these channels (otherwise, all channels are included)
fc fc=11,15 FC value(s), e.g. to find 11 and 15 Hz spindles (default 13.5 Hz)
cycles cycles=12 Number of cycles, where more cycles gives better frequency but poorer temporal resolution (default 7)
th th=6 Multiplicative threshold for core spindle detection (default 4.5)
th2 th2=3 Multiplicative threshold for non-core spindle detection (default=2)
median median Flag to indicate that the median, not mean, is used for thresholding
q q=0.3 Quality metric criterion for individual spindles (default 0)
Secondary parameters

Most users will not need to alter these.

Parameter Example Description
fc-lower fc-lower=9 Lower limit if iterating over multiple FC values
fc-upper fc-upper=15 Upper limit if iterating over multiple FC values
fc-step fc-step=2 Increment step if iterating over multiple FC values
th-max th-max=10 Maximum threshold for spindle core (default: none)
min min=1 Minimum duration for an entire spindle (default 0.5 seconds)
min0 min0=0.3 Minimum duration for a spindle core (default 0.3 seconds)
max max=2 Maximum duration for an entire spindle (default 3 seconds)
win win=0.2 Smoothing window for wavelet coefficients (default 0.1 seconds)
local local=120 Use local window (in seconds) to define baseline for spindle detection


Individual-level output (strata: F x CH)

Variable Description
DENS Spindle density (count per minute)
N Total number of spindles
AMP Mean spindle amplitude (uV or mV units)
DUR Mean spindle duration (core+flanking region)
NOSC Mean number of oscillations per spindle
FWHM Mean spindle FWHM (full width at half maximum)
ISA_S Mean integrated spindle activity (ISA) per spindle
ISA_M Mean integrated spindle activity (ISA) per minute
ISA_T Total integrated spindle activity (ISA)
FRQ Mean spindle frequency (from counting zero-crossings)
FFT Mean spindle frequency (from FFT)
CHIRP Mean chirp metric per spindle
SYMM Mean spindle symmetry metric
SYMM2 Mean spindle folded-symmetry metric
Q Mean spindle quality metric
DISPERSION Mean dispersion index of epoch spindle count
DISPERSION_P P-value for test of over-dispersion
MINS Total duration of signal entered into the analysis (minutes)
NE Number of epochs
N01 Number of spindles prior to merging
N02 Number of spindles post merging, prior to QC

Epoch-level output (strata: E x F x CH)

Variable Description
N Number of spindles observed in that epoch (for that target frequency/channel)

Spindle-level output (strata: SPINDLE x F x CH)

Variable Description
AMP Spindle amplitude (uV or mV units)
CHIRP Spindle chirp (-1 to +1)
DUR Spindle duration (seconds)
FWHM Spindle FWHM (seconds)
NOSC Number of oscillations
FRQ Spindle frequency based on counting zero-crossings in bandpass filtered signal
FFT Spindle frequency based on FFT
ISA Integrated spindle activity
MAXTSTAT Maximum wavelet statistic
MEANSTAT Mean wavelet statistic
Q Quality metric
PASS Flag (0/1) for whether this spindle passes the quality metric criterion
START Start position of the spindle (seconds elapsed since start of EDF)
STOP Stop position of the spindle (seconds elapsed since start of EDF)
START_SP Start position of the spindle (in sample-units relative to current in-memory EDF)
STOP_SP Stop position of the spindle (in sample-units relative to the current in-memory EDF)
SYMM Symmetry index (relative position of peak)
SYMM2 Folded symmetry index (0=symmetrical, 1=asymmetrical)


Here we estimate spindles for all NREM2 sleep for the tutorial individual nsrr02:

luna s.lst nsrr02 -o out.db -s "MASK ifnot=NREM2 & RE & SPINDLES fc=11,15 sig=EEG"

We see some output is sent to the console describing the process. For slow (i.e. target frequency of 11 Hz) spindles:


 detecting spindles around F_C 11Hz
 wavelet with 7 cycles
 smoothing window = 0.1s
 detection thresholds (core, flank, max)  = 4.5, 2x
 core duration threshold (core, min, max) = 0.3, 0.5, 3s
 pruned spindle count from 296 to 275
 filtering at 9 to 13
 filtering channel EEG_BP_11_4, bandpass FIR order 71
 estimated spindle density is 1.37845

The relevant strata in out.db are given by destrat, namely per-individual, per-epoch and per-spindle respectively:

destrat out.db
  [SPINDLES]    : F CH              : 2 level(s)    : AMP CHIRP DENS DISPERSION DISPERSION_P
                :                   :               : DUR FFT FRQ FWHM ISA_M ISA_S ISA_T
                :                   :               : MINS N N01 N02 NE NOSC Q SYMM SYMM2
                :                   :               :
                :                   :               : 
  [SPINDLES]    : E F CH            : (...)         : N
                :                   :               : 
  [SPINDLES]    : F CH SPINDLE      : 692 level(s)  : AMP CHIRP DUR FFT FRQ FWHM ISA
                :                   :               : MAXSTAT MEANSTAT NOSC PASS Q START
                :                   :               : START_SP STOP STOP_SP SYMM SYMM2
                :                   :               :

To get the spindle density and mean spindle duration and frequency for this individual:

destrat out.db +SPINDLES -r F CH -v DENS DUR FFT -p 3 
ID      CH    F    DENS     DUR     FFT
nsrr02  EEG   11   1.378    0.847   11.280
nsrr02  EEG   15   2.090    0.847   13.475

We will further query the output using lunaR and ldb(). In R:


k <- ldb( "out.db" )

We can use lx() to see the strata (factor combinations) that correspond to the destrat output above:

RE : BL 

Looking at the per-epoch data:

d <- lx( k , "SPINDLES" , "CH_E_F" ) 
> head(d)
      ID  CH  E  F N
1 nsrr02 EEG 92 11 2
2 nsrr02 EEG 92 15 0
3 nsrr02 EEG 93 11 4
4 nsrr02 EEG 93 15 0
5 nsrr02 EEG 98 11 2
6 nsrr02 EEG 98 15 1

Not that it is particularly informative in itself, but we can plot the number of spindles per epoch for fast and slow spindles:

par(mfcol=c(2,1), mar=c(4,4,1,1))
plot( d$E[ d$F == 11 ]  , d$N[ d$F == 11 ]  , pch= 20 , col="orange" , 
      xlab="Epoch" , ylab="N spindles" , main="11 Hz spindles" ) 
plot( d$E[ d$F == 15 ]  , d$N[ d$F == 15 ]  , pch= 20 , col="skyblue" , 
      xlab="Epoch" , ylab="N spindles" , main="15 Hz spindles" )


Luna provides a rough sanity-check to test for unusual distributions of spindles per epoch, as might occur if a movement artifact, for example, causes multiple spindles to be scored in one epoch. To a first approximation, we can treat the distribution of spindles per epoch as Poisson-distributed, and test for over-dispersion (i.e. if the variance is significantly greater than the mean). This is provided by the DISPERSION and DISPERSION_P variables in the individual-level output:

destrat out.db +SPINDLES -r F CH -v DISPERSION DISPERSION_P -p 3 
nsrr02  EEG  11   1.143        0.606
nsrr02  EEG  15   1.243        0.038

In practice, we take values of DISPERSION above 2.0 to reflect potentially problematic recordings/patterns of spindle detection, so these don't look too bad. Potentially the epochs with more than five spindles may questionable (especially as we did not apply any type of artifact correction prior to detecting spindles):

table( d$F , d$N ) 
       0   1   2   3   4   5   6   7
  11 206 131  50   7   3   1   1   0
  15 154 139  62  27  14   2   0   1

Finally, here is an example spindle detected in this individual, shown alongside the bandpass filtered EEG for delta, theta, alpha and beta power as well as sigma power. As expected (i.e. by definition), the spindle corresponds to a brief increase in sigma power:



Primarily in the context of lunaR, in the near future we'll be adding tutorials and helper functions to facilitate automatically making these types of plots quickly for all spindles in a dataset.


The core of the spindle detection algorithm involves imposing discrete thresholds on the quantitative metrics derived from the continuously-varying EEG. As such, there are a lot of magic numbers inherent in the spindle detection process. One of the key parameters is the threshold which the normalized wavelet coefficient must exceed for a putative spindle to be scored (the default of which is 4.5 times the mean). But is 4.5 a good magic number? What about 2 instead? Or what about 20? And so on...

To get some insight into whether this may be an appropriate choice for the particular individuals and analyses under consideration, Luna has a set of functions that evaluate different thresholds based on Otsu's method, a technique from computer vision that attempts to set a threshold so as to minimize the within-class variance (and so correspondingly, maximize the between class variance) in the two classes that are formed by placing a dichotomizing threshold on a continuous measure. In our scenario, the measure is the wavelet coefficient, and we are trying to place a threshold that roughly corresponds to "spindle" and "non-spindle".

If the continuous metric exhibits clear bimodality, this approach should easily find a threshold that effectively splits the two peaks. In the context of spindles, where things are likely to be a bit messier, we consider the distribution of wavelet coefficients to be a mixture of two classes: the null distribution of noise, i.e. where there is no true spindle, and a second, less frequency class representing points where a spindle is present (and where we will expect higher values of the wavelet coefficient).

Parameter Example Description
empirical empirical Empirically determine thresholds
set-empirical set-empirical Use empirically determined thresholds for spindle detection
verbose-empirical verbose-empirical Output extensive information on threshold estimation

Individual-level output (option: empirical, strata: F x CH)

Variable Description
EMPTH Empirically-determined threshold
EMPF Relative frequency of above-thresholds points based on EMPTH
MEAN_OVER_MEDIAN Ratio of mean to median, to index skewness of the wavelet coefficients

Between-class variance over range of thresholds (option: empirical, strata: TH x F x CH)

Variable Description
SIGMAB Between-class variance for given threshold

Here we use this approach on the three tutorial individuals: we run a basic command to estimate spindles (for all NREM2 sleep, with no other artifact detection in place), adding the empirical threshold:

luna s.lst sig=EEG -o out.db -s "MASK ifnot=NREM2 & RE & SPINDLES fc=11,15 empirical" 

In lunaR, we use ldb() to load the resulting out.db file:

k <- ldb("out.db")
We extract the relevant table of results by different thresholds (TH):

d <- lx( k , "SPINDLES" , "CH_F_TH" )
As this table contains all three individuals, we'll use lid() and basic R functions to create three plots separately, of the between-class variance as a function of threshold, split by fast and slow spindles. First, we get a list of the individual IDs:

ids <- unique( d$ID ) 
Then we loop to create three plots:
for (id in ids) { 
t <- lid( d , id )  
plot( t$TH[ t$F==11 ] ,  t$SIGMAB[ t$F == 11 ] , 
      type="l" , lwd=2 , col="orange" , 
      main=id, xlab="Threshold" , ylab="Between-class variance" , ylim=c(0,1)) 
lines( t$TH[ t$F==15 ] ,  t$SIGMAB[ t$F == 15 ] , type="l" , lwd=2 , col="skyblue" ) 
legend( 5,0.2,c("11 Hz","15 Hz"),fill=c("orange","skyblue"))


For the first two individuals, based on Otsu's criterion at least, this suggests that thresholds around 4 or 5 are optimal, in the sense of maximizing the between class variance (i.e. in wavelet coefficients between "spindle" and "non-spindle" classes). For the third individual, the results look strange, suggesting that extreme threshold such as 20 times the mean is optimal. This is mirrored in the individual-level output for EMPTH and EMPF: the optimal thresholds and the corresponding area under the curve:

destrat out.db +SPINDLES -r F CH -v EMPTH EMPF -p 3
ID       CH   F   EMPF   EMPTH
nsrr01   EEG  11  0.034  4.000
nsrr01   EEG  15  0.067  3.250
nsrr02   EEG  11  0.041  4.000
nsrr02   EEG  15  0.024  6.250
nsrr03   EEG  11  0.001  20.000
nsrr03   EEG  15  0.002  20.000

This suggests there may be something off with the EEG channel for nsrr03. As well as looking at the raw data, one approach may be to impose some automated artifact detection prior to spindle detection and threshold estimation (i.e. as we would typically recommend). Re-running with aberrant epochs removed prior to spindle detection (i.e. by including ARTIFACTS and SIGSTATS commands:

luna s.lst sig=EEG -o out2.db -s "MASK ifnot=NREM2 & RE \
                                  & ARTIFACTS mask & SIGSTATS mask th=3,3,3 & RE \
                                  & SPINDLES fc=11,15 empirical" 

This appears to normalize the "optimal" thresholds somewhat, in that the area under the curve (i.e. roughly corresponding to the portion of time in "spindle" versus "non-spindle" states) are ball-park what we may expect to observe in a typical individual who has around 2 or so spindles per minute detected at a central scalp EEG channel, with each spindle lasting approximately 1 second (i.e. 2/60 = 0.033):

destrat out2.db +SPINDLES -r F CH -v EMPTH EMPF -p 3
ID      CH   F    EMPF    EMPTH
nsrr01  EEG  11   0.046   3.500
nsrr01  EEG  15   0.092   2.750
nsrr02  EEG  11   0.041   4.000
nsrr02  EEG  15   0.029   5.500
nsrr03  EEG  11   0.033   4.000
nsrr03  EEG  15   0.066   3.000

Returning to the plots, things look a bit better:


What to make of all this? We do not advise selecting different "optimal" thresholds for each individual based on these analyses. Rather, we suggest that looking at these (and other) metrics can be useful in basic sanity-checking for whether a given analysis for a given individual is likely to yield meaningful results.

Quality metrics

Luna calculates a simple QC metric for each spindle, by considering the relative increase of non-spindle activity within the spindle interval, relative to spindle-activity. We use five fixed bands, which unlike the PSD command, are not altered by changing special variables. The three non-spindle bands are: delta (0.5-4 Hz), theta (4-8 Hz), and beta (defined here as 20-30 Hz). The two spindle bands are slow sigma (10-13.5 Hz) and fast sigma (13.5-16 Hz).

For each channel, Luna calculates the baseline power for each band, Pband. For the interval spanning the ith spindle, we calculate band power Piband. We then calculate relative enrichment for each band (on a log10 scale), as:

    Eiband = log10(Piband) - log10(Pband).

The Q score for each spindle is based on the difference between the largest spindle-band relative enrichment versus the largest non-spindle band relative enrichment:

    Qi = max( Eislow-sigma , Eifast-sigma ) - max( Eidelta , Eitheta , Eibeta )

That is, numbers less than zero indicate that a non-spindle band is showing greater relative enrichment (relative to the baseline log-scaled power for that band) than either spindle band. Such spindles are more likely to reflect artifact which had a non-specific effect on the power spectrum.

To illustrate this QC metric, we'll detect spindles for the second tutorial individual:

luna s.lst 2 sig=EEG -o out.db \
  -s "MASK ifnot=NREM2 & RE & \
      SPINDLES fc=11,15 enrich q=-9 ftr=sp1 cycles=12"

Here, we detect both fast and slow spindles, targeting 11 and 15 Hz respectively; we've set cycles to 12 instead of the default 7, in order to increase frequency-specificity. We've added the enrich option to output Eislow-sigma and Eifast-sigma along with Eidelta Eitheta and Eibeta for each spindle/channel/target-frequency/band (and so stratified by CHxFxSPINDLExB). We also set q equal to -9, which effectively turns off any QC filtering, so that we can see what some spindles that would not otherwise have been called with the default threshold (q=0) look like. Finally, the ftr option writes out two files that detail where spindles occur.

In lunaR:

lattach( lsl( "s.lst" ) , 2 )
We'll attach as annotations the two newly-created files:
ladd.annot.file( "" ) 
ladd.annot.file( "" ) 
We'll then load our previously-generated results, saved in out.db:
k <- ldb( "out.db" ) 
and extract the table which contains Qi, the per-spindle quality metric.
d <- lx( k , "SPINDLES" , "F" , "CH" , "SPINDLE" )
Plotting the distribution of this separately for slow and fast spindles:
hist( d$Q[ d$F == 11 ] , breaks=30 , main="11 Hz spindles" , xlab="Q" , xlim=c(-2,2) ) 
hist( d$Q[ d$F == 15 ] , breaks=30 , main="15 Hz spindles" , xlab="Q" , xlim=c(-2,2) )
We see that most spindles (as would be expected) have quality scores greater than 0, meaning that the interval showed greater spindle-band enrichment than non-spindle-band enrichment.


Of the full set (i.e. running without the default q=0 threshold in place), we see that approximately 10% of spindles have q less than 0:

table( d$Q > 0 , d$F ) 
         11  15
  FALSE  68  36
  TRUE  405 370

Let's visualize some q<0 and q>0 spindles (we'll actually select qc>1 spindles, to ensure that we don't randomly select spindles with q very near 0). For simplicity, we'll do this just for fast spindles.

First, we'll save the labels of the two new annotation tracks we appended (i.e. you can see the names from running lannots()):

a15 = "SPINDLES-spindles-EEG-wavelet-15-sp1"
We'll then extract the intervals of each annotation instance, i.e. each spindle:
s15 <- lannots( a15 )
Confirming that we get 406 spindles, i.e. as output above, we can see length(s15) is indeed 406.

For visualization, here we'll just select five low Q and five high Q spindles:

lo15 <- sample(d$SPINDLE[ d$Q < 0 & d$F == 15 ])[1:5]
hi15 <- sample(d$SPINDLE[ d$Q > 1 & d$F == 15 ])[1:5]

Using lunaR functions, we can then write a function to plot the raw EEG signal for each spindle, with a 4-second flanking window on either side:

f1 <- function(i,a) { 
x <- ldata.intervals( i , chs = "EEG" , annots = a , w = 4 )
y <- x$EEG
plot( x$SEC , y , type="l" , axes=F,xlab="",ylab="" , ylim=c(-125,125) )
t <- x$SEC[ x[,lsanitize(a)] == 1 ]
points( t , rep( -125,length(t)) , col=rgb(0,0,255,150,max=255) , pch="|" )

Finally, we'll generate a matrix of plots, with the five low Q spindles in the left column, and the five high Q spindles in the right column:

par(mfrow=c(5,2),mar=c(0,0,0,0) )
for (i in 1:5) {
f1( i = s15[ lo15[ i ] ] , a = a15 )
f1( i = s15[ hi15[ i ] ] , a = a15 )

As you can see, the low Q spindles (that would not have been detected using Luna's default q=0 setting) appear to reflect some sources of artifact, whereas, as a whole, the high Q spindles look more like canonical spindles.


Be aware that setting q too high may unnecessarily miss many true spindles though. A value of q=0.3 (i.e. log10(2)) requires that the spindle-band enrichment is at least twice non-spindle band enrichment; a value of q=1 requires a 10-fold enrichment. Considering the distribution of mean Q over all individuals in the study may provide a good starting point, if you want to select non-default values.


Naturally, neither this particular automated spindle detection algorithm, nor this particular quality metric, are guaranteed to always get it right... We fully expect that this algorithm will flag as a "spindle" events that most scorers would not, and vice versa. When first approaching this tool, we'd advise experimenting with the different thresholds and parameters, and spot-check studies in the way we have above. (The lunaR tutorial shows how to automatically iterate through all spindles detected.)

In general, if visual review seems to indicate likely artifact, the spindle is likely to be a false positive. In contrast, however, otherwise high Q spindles may appear relatively modest on visual review, and may not have been detected by some human scorers. Human scoring does not constitute a real gold-standard, however, and both human and automatic approaches impose inherently arbitrary thresholds on what is likely a continuously-varying and heterogeneous set of physiological processes. As such, we'd advise against automatically assuming the automatic detector is wrong, at least for this latter class of events ("low amplitude spindles").

Overall, the usual principles of a) trying to ensure that high-quality data go into the analysis, and b) working with large samples and appropriate statistical techniques downstream, are likely to be just as (if not more) important than spending too much tweaking the detection procedure, or agonizing over what is or is not "a true spindle".

Collating spindles

When detecting spindles on multiple channels or at multiple, nearby frequencies, what presumably is the same underlying physiological event may often be detected multiple times. These options assist in collating the results of spindle detection across channels and/or frequencies:

  • you can control whether two temporally close spindles detected on the same channel/target-frequency should be merged and considered as one event or not (merge parameter, the default is 0.5 seconds)

  • second, within channel, you can collate temporally overlapping spindles detected at nearby frequencies (i.e. we would expect a great deal of overlap in the results of an analysis with FC=10 Hz and FC=11 Hz).

  • third, across channels, you can collate temporally overlapping spindles detected at nearby target frequencies; these may or may not represent physiologically distinct events depending on the montage and channels being collated.

These commands essentially take a list of spindles and generate another list of so-called m-spindles: groups of spindles that, on some definition, "overlap", i.e. represent the same underlying event.

The primary output from this analysis is an estimate of m-spindle density (MSP_DENS), i.e. which is taken to reflect the total number of unique spindles (i.e. allowing for potentially overlapping definitions). For each m-spindle (set a spindles), Luna also calculates an mean frequency (a weighted mean of the constituent spindle frequencies, with weights given by the ISA of each spindle). This is output for each spindle and also used to generate a frequency-conditioned estimate of m-spindle density.

Parameter Example Description
merge merge=0.2 Merge two putative spindles if within this interval (default 0.5 seconds)
collate collate Within each channel, collate overlapping spindles with similar frequencies (to MSPINDLES)
collate-channels collate-channels As above, except merge across channels also
th-frq th-frq=1 Criterion for merging two spindles if difference in frequency is within this range (default 2 Hz)
list-all-spindles list-all-spindles List all spindles that comprise each m-spindle

Secondary parameters:

Parameter Example Description
th-interval th-interval=0.5 Merge two overlapping spindles if the ratio of intersection to union is at least this value (default 0, i.e. any overlap)
window window=0.5 Set window around each spindle when defining temporal overlap
hms hms Show clock-time of each m-spindle

Individual-level summaries of m-spindles (option: collate, strata: none)

Variable Description
MSP_DENS m-spindle density
MSP_N m-spindle count
MSP_MINS Denominator for density, i.e. minutes of signal analyzed

m-spindle density stratified by m-spindle frequency (option: collate, strata: F )

Variable Description
MSP_DENS m-spindle density conditional on m-spindle frequency

Merged-spindle output (option: collate, strata: MSPINDLE)

or Merged-spindle output (option: collate-within-channel, strata: CH x MSPINDLE)

Variable Description
MSP_DUR Duration of this m-spindle
MSP_F Estimated frequency of this m-spindle
MSP_FL Lower frequency of this m-spindle
MSP_FU Upper frequency of this m-spindle
MSP_SIZE Number of spindles in this m-spindle
MSP_STAT Statistic for m-spindle
MSP_START Start time (seconds elapsed from EDF start) of m-spindle
MSP_STOP Stop time (seconds elapsed from EDF start) of m-spindle

Spindle to m-spindle mappings (option: list-all-spindles, strata: SPINDLE x MSPINDLE)

Variable Description
SCH Spindle label (channel:target frequency)
FFT Spindle estimated frequency (via FFT)
START Spindle start time (elapsed seconds from EDF start)
STOP Spindle stop time (elapsed seconds from EDF start)

Additional output (option: hms, strata: MSPINDLE or CH x MSPINDLE))

Variable Description
MSP_START_HMS Merged spindle start clock-time
MSP_STOP_HMS Merged spindle stop clock-time
luna s.lst 2 sig=EEG -o out.db -s "MASK ifnot=NREM2 & RE & \
                                   SPINDLES fc-lower=10 fc-upper=16 fc-step=0.25 \
                                            collate list-all-spindles"

The total m-spindle density is 4.89 spindles per second, i.e. this estimates the total unique number of spindles over this broad frequency range ( with target frequencies of 10-16 Hz and each wavelet spanning approximately plus/minus 2 Hz):

destrat out.db +SPINDLES 
nsrr02   4.8972     199.5       977

We can also view the m-spindle density as a function of the estimated m-spindle frequency: i.e. these numbers will sum to the total MSP_DENS estimate (or close to it, as we only consider the fixed range of 8 to 16 Hz in 0.5 Hz bins, by default, when outputing this):

destrat out.db +SPINDLES -r F 

Plotting the output of this command, we potentially see evidence of a bimodality in spindle frequencies for this individual:


Relationship to individualized spindle detectors

Rather than using fixed target frequencies (e.g. 13.5 Hz, 11 Hz, or 15 Hz, etc) to define "spindles", some approaches to spindle detection seek to find, empirically, the "peak" or "peaks" (which may or may not be obviously present in the power spectra) for each particular individual, and then target spindles of those particular frequencies for that individual. In contrast, this approach is effectively agnostic to an individual's "true" characteristic spindle frequency (or frequencies), or to having to define spindle frequencies up front, by virtue of simply considering all possible frequencies in the broad sigma range, and then combining results. (Arguably this approach might be preferred a priori in any case, on the grounds that even within a single individual we can see marked shifts in spindle frequency over the course of the night.)

Other metrics

The following options report some additional spindle morphology characteristics (these are unlikely to be desired as part of any primary spindle analysis, and so can be ignored by most users). The if option estimates spindle instantaneous frequency, averaged over spindles, as a function of the relative location within the spindle. The tlock option produces an averaged EEG signature of detected spindles, synced to the spindle peak (point of max peak-to-peak amplitude).

Parameter Example Description
if if Estimate instantaneous frequency of spindles
if-frq if-frq=1 Window around target frequency (default 2 hz)
tlock tlock Flag to request (verbose) average, peak-locked waveforms

Instantaneous frequency is estimated via the filter-Hilbert method, where the filter is FC +/- H Hz (where H is set by if-frq, default 2).


Instantaneous frequency (IF) per spindle (option: if, strata: CH x F x SPINDLE)

Variable Description
IF Mean frequency per spindle over duration

Mean IF stratified by relative location in spindle (option: if, strata CH x F x RELLOC)

Variable Description
IF Mean frequency of all spindles, per relative position within the spindle (five bins)

With the tutorial data:

luna s.lst 2 -o out.db -s "MASK ifnot=NREM2 & RE & \
                           SPINDLES sig=EEG fc=11,15 if tlock q=0.3 cycles=12"
To get the estimate of IF per spindle (which will should typically be very highly correlated with FRQ and FFT, as they are all estimating essentially the same thing):

destrat out.db +SPINDLES -r F CH SPINDLE -v FFT FRQ IF 

To get the mean IF stratified by relative position with the spindle (i.e. dividing each spindle into five, equally-sized time bins):

destrat out.db +SPINDLES -r F CH RELLOC -p 3
nsrr02  EEG 11  1   11.321
nsrr02  EEG 11  2   11.218
nsrr02  EEG 11  3   11.041
nsrr02  EEG 11  4   10.912
nsrr02  EEG 11  5   10.818
nsrr02  EEG 15  1   14.287
nsrr02  EEG 15  2   14.119
nsrr02  EEG 15  3   13.959
nsrr02  EEG 15  4   14.015
nsrr02  EEG 15  5   14.181

To plot the mean signal time-locked to each spindle peak (in lunaR):

k <- ldb("out.db")

d <- lx( k , "SPINDLES" , "CH" , "F" , "MSEC" ) 

par( mfcol=c(1,2) )
for (f in c(11,15)) { 
 plot( d$MSEC[ d$F==f ] , d$TLOCK[ d$F==f ] , ylim=c(-20,20) , type="l"
       xlab="Time (sec)",ylab="Mean EEG",main=paste(f,"Hz spindles"), lwd=2)



Luna can generate FTR files representing the spindles detected in a given run. These can subsequently be attached to an EDF via lunaR, for example, in order to visualize spindles, as in the examples above.

Parameter Example Description
ftr ftr=tag Produce FTR files for all spindles, with the tag in the filename
ftr-dir ftr-dir=/path/to/folder Folder for FTR files
show-coef show-coef Flag to request (very verbose) coefficient output (to stdout)

The ftr option generates one file per channel/target frequency combination, with a filename in the form:


where ID is the EDF ID, CH is the channel label, FC is the target frequency, and tag is the user-supplied tag (from the ftr option).

Spindle/SO coupling

The sw option will use a heuristic method to detect slow oscillations (SO) in the sleep EEG, along with spindles, and will quantify the coupling between spindles and SO, via event-related phase locking measures. There are a number of additional parameters that modify the behavior of sw.


The documentation for spindle/SO coupling, and for SO detection in general is still skeletal. We're working on it...


Bandpass filtering for SO detection:

Parameter Example Description
f-lwr f-lwr=0.2 Bandpass filter lower transition frequency
f-upr f-upr=4.5 Bandpass filter upper transition frequency

Time criteria for SO detection:

Parameter Example Description
t-lwr t-lwr= Minimum duration of a SO
t-upr t-upr= Maximum duration of a SO
t-neg-lwr t-neg-lwr=0.5 Lower threshold for negative peak (seconds)
t-neg-upr t-neg-lwr=2 Upper threshold for negative peak (seconds)

Amplitude criteria for SO detection:

Parameter Example Description
uV-neg uV-neg=-75 Negative peak amplitude threshold ( default 0, i.e. no threshold )
uV-p2p uV-p2p=140 Peak-to-peak amplitude (default 0, i.e. no threshold)
mag mag=1 Use a relative amplitude threshold instead (see below)
sw-mean sw-mean Use mean (instead of median)

Additional SO detection options (currently not fully supported, do not use):

Parameter Example Description
neg2pos neg2pos Define full SO based on consecutive negative-to-positive zero-crossings
half-wave half-wave Detect all half waves instead of full waves
negative-half-wave negative-half-wave Detect negative half-waves instead of full waves
positive-half-wave positive-half-wave Detect positive half-waves instead of full waves

See the output tables below for the SO command for a description of the additional output generated by the sw option of SPINDLES. In addition, there are a number of additional tables detailing spindle/SO coupling generated by the sw option.

Primary individual-level spindle/SO coupling output (option: sw, strata: CH x F)

Variable Description
ITPC_ANGLE_PEAK Circular mean of SO phase at spindle peak (SW_SPHASE_PEAK below)
ITPC_PEAK Intertrial (phase) coherence (ITPC) statistic
PROP_ON_SW Proportion of spindles temporally overlapping a detected SO
PROP_PEAK_ON_SW Proportion of spindles whose peak temporally overlaps a detected SO

Spindle-level spindle/SO coupling output (option: sw, strata: CH x F x SPINDLE)

Variable Description
SW_SPHASE_PEAK SO phase at spindle peak
NEAREST_SW Time to nearest SO
NEAREST_SW_NUM Number of nearest SO

SO morphology (option: sw, strata: CH x PHASE)

Variable Description
SWPL_EEG Phase-locked mean EEG value (as a function of SO phase, in 20-deg bins)

SO morphology (option: sw, strata: CH x SP)

Variable Description
SWTL_EEG Time-locked mean EEG value (as a function of sample-points before/after SO negative peak)

Spindle/SO phase coupling (option: sw, strata: CH x F x PHASE

Variable Description
SWPL_CWT SO-phased locked mean spindle wavelet coefficient (18 20-deg bins)

Spindle/SO phase coupling (option: sw, strata: CH x F x SP

Variable Description
SWTL_CWT Time-locked mean spindle wavelet coefficient

In the literature, different investigators have used different definitions for SOs in different contexts. To detect SO based on relative amplitude of waveforms (relative to the whole night/stage for that individual), you could use something like (note, these are options added to a SPINDLES command, although we are not showing the full command here):

sw f-lwr=0.3 f-upr=4 t-lwr=0.5 t-upr=2 mag=1

Here, the EEG channel is first band-pass filtered with transition frequencies 0.3 and 4 Hz. After counting all subsequent positive-to-negative zero-crossings, only intervals between 0.5 and 2 seconds in duration are retained as putative SO. Further, only intervals where both the negative peak amplitude and peak-to-peak amplitude are at least 1.0 (mag) times the median for that individual are designated as the final set of SO.

Instead of using a relative amplitude threshold (mag), you can set absolute thresholds. For example, delta waves with a negative peak of at least -40 uV, and at least 75 uV peak-to-peak:

sw f-lwr=0.3 f-upr=4 uV-neg=-40 uV-p2p=75
For larger slow waves, one might use something like:
sw f-lwr=0.3 f-upr=4 uV-neg=-75 uV-p2p=140

You can also specify separate thresholds for the duration of the negative and positive half waves:

t-neg-lwr=0.3 t-neg-upr=1.5  t-pos-lwr=0.0 t-pos-upr=1.0


This section is still under development. We'll add a specific vignette on using Luna for spindle/SO coupling analyses in the near future.



The show-coef is primarily for development purposes and should not be part of a standard spindle-analysis workflow. It produces a lot of output: five values for every time-point for every channel for every spindle frequency: this will be hard work for the output database, which is not really designed for this type of data: you are better off using this without the -o database option; rather, extract the desired information directly from standard output:

luna s.lst 1 < cmd.txt \
  | awk ' $2 == "SPINDLES" && $5 == "CWT" { print $3,$4,$6 } ' OFS="\t"
Variable Description
RAWCWT Absolute, raw wavelet coefficient value
CWT Normalized wavelet coefficient value
CWT_TH Wavelet threshold for core spindle
CWT_TH2 Wavelet threshold for flanking spindle
CWT_THMAX Wavelet max threshold for spindle


Detects slow oscillations via a simple heuristic approach

The SO command uses a a heuristic approach to detect slow oscillations in the EEG, combined with a filter-Hilbert approach to estimate SO phase. Currently can only be applied to a single channel at a time.


see spindle/SO coupling section above

Additional options;

Parameter Example Description
tl tl=sig1,sig2 Other signals to be summarized w.r.t. time-locked SO
onset onset
pos pos
window window=1 Window around SO peak (default +/- 3 seconds)


Note: the output uses SW (slow wave) as equivalent to SO (slow oscillation). Also, additional variables starting SW_MD reflect the median rather than the mean across events.

SO information per channel (option: sw, strata: CH)

Variable Description
SW Number of SO detected
SW_RATE SO per minute
SW_TH_NEG Threshold used for negative peak
SW_TH_P2P Threshold used for peak-to-peak
SW_DUR Median SO duration
SW_NEG_DUR Median negative peak duration
SW_POS_DUR Median positive peak duration
SW_AMP Median amplitude (of negative peak)
SW_P2P Median peak-to-peak amplitude
SW_SLOPE_NEG1 Median slope (positive-to-negative zero-crossing to negative peak)
SW_SLOPE_NEG2 Median slope (negative peak to negative-to-positive zero-crossing)
SW_SLOPE_POS1 Median slope (negative-to-positive zero crossing to positive peak)
SW_SLOPE_POS2 Median slope (positive peak to positive-to-negative zero-crossing)

Per SO output (strata: CH x N)

Variable Description
START_IDX Start of SO (in sample-point units)
STOP_IDX Stop of SO (in sample-point units)
START Start of SO (in seconds elapsed from start of EDF)
STOP Stop of SO (in seconds elapsed from start of EDF)
DUR SO duration
UP_AMP SO positive peak amplitude
DOWN_AMP SO negative peak amplitude
P2P_AMP SO peak-to-peak amplitude
UP_IDX Sample index of positive peak
DOWN_IDX Sample index of negative peak
SLOPE_POS1 Slope (as above) (unless half-waves)
SLOPE_POS2 Slope (as above) (unless half-waves)
SLOPE_NEG1 Slope (as above) (unless half-waves)
SLOPE_NEG2 Slope (as above) (unless half-waves)

Verbose SO-level out (strata: CH x N x SP)

Variable Description
FLT Filtered wave-form for each SO
PH Instantaneous phase for each SO

Epoch-level counts of SO (strata: CH x E)

Variable Description
N Number of SOs detected in that epoch

From tl ( strata: CH x CH2 x SP)

Variable Description


to be added