This page was created by the IDL library routine `mk_html_help2`.

Last modified: Thu Sep 11 18:17:40 2025.

# Directory Listing of Routines

# general/science/spd_mtm

- SPD_MTM
- SPD_MTM_ADAPT
- SPD_MTM_CONFTHR
- SPD_MTM_DBL2FLT
- SPD_MTM_DISPAR
- SPD_MTM_DPSS
- SPD_MTM_FINDPEAKS
- SPD_MTM_FITMODEL
- SPD_MTM_MAKEPLOT
- SPD_MTM_MODELGOF
- SPD_MTM_REGRE
- SPD_MTM_RESHAPE
- SPD_MTM_SMOOTHING
- SPD_MTM_SPEC
- TPLOT_SPD_MTM

# Routine Descriptions

## SPD_MTM

[Next Routine] [List of Routines]

     :Description:
       Perform a spectral analysis following
       Di Matteo, S., Viall, N., Kepko, L. (2020),
       "Power Spectral Density Background Estimate and
       Signals Detection via the Multitaper Method."
       Journal of Geophysical Research Space Physics.

       PSD stands for Power Spectral Density

     :Params:
       INPUTS:
          data - 2 column data: [time vector, time series]
                 N.B. The average value of the data is removed by default.
                 If the data are not evenly sampled, the average time step
                 is used when it is lower than its standard deviation.
            NW - time-halfbandwidth product, DEFAULT: NW = 3
          Ktpr - Number of tapers to apply, DEFAULT: Ktpr = 2*NW - 1
       padding - Amount of padding to apply (proportional to the original data)
                 DEFAULT: padding=1 so that the frequency step (df)
                 equals the Rayleigh frequency [1/(N*dt)]
          dpss - structure with tapers and eigenvalues
                 evaluated by spd_mtm_dpss.pro or previous call of spd_mtm.pro
          flim - limits in percentages of the frequency range to be analyzed,
                 DEFAULT: [df,fny-df], where fny is the Nyquist frequency [1/(2*dt)]
     smoothing - possible smoothing:
                 'raw'(0) no smoothing (DEFAULT)
                 'med'(1) running median
                 'mlog'(2) running median (constant window on log(f) )
                 'bin'(3) binned PSD
                 'but'(4) low pass filtered PSD (Butterworth filter)
                 'all'(9) use all the smoothing procedures
       psmooth - value between (0,0.5] defining the "smoothing window"
                 psmooth=2 -> search the optimum window according
                              to the Kolmogorov-Smirnov test (DEFAULT)
        model  - implemented models:
                 'wht'(0) white noise, [N]
                 'pl'(1) power law, [N, beta]
                 'ar1'(2) first order autoregressive process, [N, rho]
                 'bpl'(3) bending power law, [N, beta, gamma, fb] (DEFAULT)
                 'all'(9) use all the models
      procpeak - possible choices:
                 '' do not search for peaks
                 'gt'(0) gamma test only
                 'ft'(1) F test only
                 'gft'(2) gamma and F test
                 'gftm'(3) only max F value for each PSD enhancement
                 'all'(9) use all the selection procedures (DEFAULT)
          conf - confidence levels in increasing order
                 (DEFAULT is conf=[0.90,0.95,0.99])
          resh - value between (0,1.0], reshape the PSD from peaks identified
                 according to the 'gft' procedure at the "resh" confidence level
                 (DEFAULT: resh=0, do not perform the reshaping)
           gof - PSD background chosing criterium:
                 'MERIT' - MERIT function (DEFAULT)
                 'CKS' - Kolmogorov-Smirnov test
                 'AIC' - Akaike Information Criteria

       x_label - label for x-axis or defined set of choice
                 (default choice 'Time', see spd_mtm_makeplot.pro)
       x_units - x variable unit of measures
       y_units - y variable unit of measures
       f_units - "frequency" variable unit of measures
        x_conv - conversion factor for the x variable
                 (N.B. IS UP TO THE USER TO BE CONSISTENT WITH X_UNITS)
        f_conv - conversion factor for the "frequency" variable
                 (N.B. IS UP TO THE USER TO BE CONSISTENT WITH F_UNITS)

       OUTPUTS:
       spec.
           .ff       Fourier frequencies
           .raw      adaptive multitaper PSD
           .back     best PSD background among the probed ones
           .ftest    values of the F test
           .resh     reshaped PSD (if resh in spd_mtm is imposed)
           .dof      degree of freedom at each Fourier frequency
           .fbin     binned frequencies for the bin-smoothed PSD
           .smth     smoothed PSD, spec.smth[#smth, *]:
                     #smth = 0->'raw'; 1->'med'; 2->'mlog', 3->'bin'; 4->'but'
           .modl     PSD background models, spec.modl[#smth, #modl, *]:
                     #modl = 0->'wht'; 1->'pl'; 2->'ar1'; 3->'bpl'
           .frpr     model parameters resulting from the fitting procedures
                     spec.frpr[#smth, #modl, #par]:
                     #par = 0->'c'; 1->'rho' or 'beta'; 2->'gamma'; 3->'fb'
           .conf     confidence threshold values for the PSD (gamma test)
           .fconf    confidence threshold values for the F statistic (F test)
           .CKS      C_KS value for each probed model spec.CKS[#smth, #modl]
           .AIC      AIC value for each probed model spec.AIC[#smth, #modl]
           .MERIT    MERIT value for each probed model spec.MERIT[#smth, #modl]
           .Ryk      real part of the eigenspectra spec.Ryk[#Ktpr,*]
           .Iyk      imaginary part of the eigenspectra spec.Iyk[#Ktpr,*]
           .psmooth  percentage of the frequency range defining
                     the smoothing window spec.psmooth[#smth]
           .muf      complex amplitude at each Fourier frequency
           .indback  indices of the selected PSD background;
                     smoothing/model = spec.indback[0]/spec.indback[1]
           .poor_MTM flag for failed convergence of the adaptive MTM PSD

       peak.
           .ff       Fourier frequencies
           .pkdf     for each peak selection method and confidence level
                     a value greater than zero at a specific frequency
                     indicate the occurence of a signal at that frequency:
                     peak.pkdf[#peakproc, #conf, #freq]
                     'gamma test' -> peak.pkdf[0,*,*] contains the badwidth
                     of the PSD enhancements at the identified frequencies
                     'F test', 'gft', and 'gftm' -> peak.pkdf[1:3,*,*]
                     is equal to par.df at the identified frequencies

        par.
           .npts     time series number of points
           .dt       average sampling time
           .dtsig    standard deviation of the sampling time
           .datavar  variance of the time series
           .fray     Rayleigh frequency: fray = 1/(npts*dt)
           .fny      Nyquist frequency: fny = 1/(2*dt)
           .npad     time series number of points after padding
           .nfreq    number of frequencies
           .df       frequency resolution (after padding), it corresponds to
                     the Rayleigh frequency for no padding (padding = 1)
           .fbeg     beginning frequency of the interval under analysis
           .fend     ending frequency of the interval under analysis
           .psmooth  value imposed in spd_mtm
           .conf     confidence thresholds percentages
           .NW       time-halfbandwidth product
           .Ktpr     number of tapers (max value is 2*NW - 1)
           .tprs     discrete prolate spheroidal sequences (dpss)
           .v        dpss eigenvalues

       ipar.
           .hires    keyword /hires is selected (1) or not (0)
           .power    keyword /power is selected (1) or not (0)
           .specproc array[#smth,#modl] smoothing and model combinations
                     1 (0) the smoothing + model combination is (not) probed
           .peakproc array[4] referring to 'gt', 'ft', 'gft', and 'gftm'
                     1 (0) peaks according to this procedure are (not) saved
           .allpkwd  keyword /allpeakwidth is selected (1) or not (0)
           .pltsm    keyword /pltsmooth is selected (1) or not (0)
           .resh     confidence threshold percentage chosen to select the
                     PSD enhancements to be removed in the PSD reshaping
           .gof      criterium used to select the PSD background

       rshspec  provide the spec structures based on the reshaped PSD
                (If it is not specified, the results are overwritten on spec)

       rshpeak  provide the peak structures based on the reshaped PSD
                (If it is not specified, the results are overwritten on peak)

     :Keywords:
         /hires - spec.raw is the high-resolution PSD estimate
         /power - spec.raw is the integrated PSD: df*PSD [##^2]
         /quiet - do not display parameters and results
         /allpeakwidth - no constrains on the PSD enhancements width
         /pltsmooth - plot only the smoothed PSDs and stop
         /makeplot - plot the results
         /double - output in double precision

     :Uses:
         spd_mtm_dpss.pro
         spd_mtm_spec.pro
         spd_mtm_regre.pro
         spd_mtm_smoothing.pro
         spd_mtm_fitmodel.pro
         spd_mtm_modelgof.pro
         spd_mtm_conflvl.pro
         spd_mtm_findpeaks.pro
         spd_mtm_reshape.pro
         spd_mtm_dispar.pro
         spd_mtm_dbl2flt.pro
         spd_mtm_makeplot.pro

     :Example:

         t = [0:511] ; time vector, suppose dt = 1s
         x = 0.5*cos(2.0*!pi*t/8.0) + randomn(3,512,1) ; time series
         data = [[t],[x]]

         spd_mtm, data=data, NW=3, Ktpr=5, padding=1, dpss=dpss, $
                  flim=[0,1], smoothing=’all’, psmooth=2, $
                  model=’wht’, procpeak=[’gt’,'gft'], $
                  conf=[0.90,0.95d], $
                  /makeplot, $
                  x_label=’Time’, y_units=’##’, $
                  x_units=’min’, x_conv=1.0/60.0, $
                  f_units=’mHz’, f_conv=1d3, $
                  spec=spec, peak=peak, par=par, ipar=ipar

         ; plot the adaptive MTM PSD
         plot(spec.ff, spec.raw, /ylog)

         ; plot the selected PSD background
         plot(spec.ff, spec.back, 'r', /overplot)

         ; plot confidence threshold for the PSD
         plot(spec.ff, spec.back*spec.conf[0], 'r--', /overplot)

         ; plot a smoothed PSD
         ; (N.B. The smoothing approach has to be present the inputs)
         plot(spec.ff, spec.smth[1,*]) ; med
         plot(spec.ff, spec.smth[2,*]) ; mlog
         plot(spec.fbin, spec.smth[3,*]) ; bin
         plot(spec.ff, spec.smth[4,*]) ; but

         ; plot a fitted PSD model on a smoothed PSD
         ; (N.B. The model has to be present the inputs)
         plot(spec.ff, spec.modl[0,0,*]) ; raw/WHT
         plot(spec.ff, spec.modl[3,1,*]) ; bin/PL
         plot(spec.ff, spec.modl[2,2,*]) ; mlog/AR(1)
         plot(spec.ff, spec.modl[3,3,*]) ; bin/BPL

         ; to recover the identified peaks:
         indices_peaks = where(peak.pkdf[0,0,*] gt 0)

         ; N.B.
         ; peak.pkdf[0,0,*] -> gamma test, lowest confidence level (from conf)
         ; peak.pkdf[2,0,*] -> gamma and F test, lowest confidence level (from conf)
         ; peak.pkdf[0,-1,*] -> gamma test, highest confidence level (from conf)

         if (indices_peaks[0] ge 0) then begin
           signals_frequency = peak.ff[indices_peaks]
         endif

     :Version:
         Version 1.0

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/spd_mtm.pro)

## SPD_MTM_ADAPT

[Previous Routine] [Next Routine] [List of Routines]

     :Description:
       Evaluate the adaptive multitaper spectrum from
       Thomson, D. J. (1982). "Spectrum estimation and harmonic analysis."
       Proceedings of the IEEE,70(9), 1055-1096. doi:10.1109/PROC.1982.12433

     :Params:
       INPUTS:
            spec.Ryk - real part of the multitaper eigenspectra
            spec.Iyk - imaginary part of the multitaper eigenspectra
               par.V - dpss eigenvalues
            par.Ktpr - number of tapers to be used (max value is 2*NW - 1)
           par.nfreq - number of frequencies
         par.datavar - variance of the data
          ipar.power - keyword /power in spd_mtm is selected (1) or not (0)

       OUTPUTS:
            spec.raw - adaptive MTM spectrum
            spec.dof - degree of freedom at each Fourier frequency
       spec.poor_MTM - flag for failed convergence

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/spd_mtm_adapt.pro)

## SPD_MTM_CONFTHR

[Previous Routine] [Next Routine] [List of Routines]

     :Description:
       Determine the confidence thresholds for the gamma and F test.

     :Params:
       INPUTS:
           par.conf - confidence thresholds percentages
           par.Ktpr - number of tapers (max value is 2*NW - 1)

       OUTPUTS:
          spec.conf - confidence threshold values for the PSD (gamma test)
         spec.fconf - confidence threshold values for the F statistic (F test)

       COMMON VARIABLE:
            gammaj
                 pp - confidence threshold percentage
     :Uses:
         spd_mtm_gammaj_cvf.pro

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/spd_mtm_confthr.pro)

## SPD_MTM_DBL2FLT

[Previous Routine] [Next Routine] [List of Routines]

     :Description:
       Turn the double precision outputs in single precision (float).

     :Params:
       INPUTS:
            spec - structure with results of the spectral analysis
                   (as defined in spd_mtm.pro)
             par - properties of the time series and parameters defining
                   the spectral analysis (as defined in spd_mtm.pro)
            peak - identified signals (as defined in spd_mtm.pro)
         rshspec - spectral analysis results for the reshaped PSD
                   (same structure as spec)
         rshpeak - signals identified with the reshaped PSD
                   (same structure as peak)

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/spd_mtm_dbl2flt.pro)

## SPD_MTM_DISPAR

[Previous Routine] [Next Routine] [List of Routines]

     :Description:
       Display parameters and results of spd_mtm.pro on IDL console.

     :Params:
       INPUTS:
           smth_label - labels of the smoothing procedures
           modl_label - labels of the PSD models
               par.dt - average sampling time
            par.dtsig - standard deviation of the sampling time
             par.fray - Rayleigh frequency: fray = 1/(npts*dt)
              par.fny - Nyquist frequency: fny = 1/(2*dt)
             par.npts - time series number of points
            ipar.resh - confidence threshold percentage chosen to select the
                        PSD enhancements to be removed in the PSD reshaping
           ipar.pltsm - keyword /pltsmooth is selected (1) or not (0)
        spec.poor_MTM - flag for failed convergence of the adaptive MTM PSD
         spec.indback - indices of the selected PSD background;
                        smoothing/model = spec.indback[0]/spec.indback[1]
            spec.frpr - model parameters resulting from the fitting procedures
                        spec.frpr[#smth, #modl, #par]:
                        #par = 0->'c'; 1->'rho' or 'beta'; 2->'gamma'; 3->'fb'

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/spd_mtm_dispar.pro)

## SPD_MTM_DPSS

[Previous Routine] [Next Routine] [List of Routines]

     :Description:
       Generate the Slepian sequences
       (or discrete prolate spheroidal sequences, dpss)
       following the tridiagonal formulation
       by Percival and Walden (1993) p.386-387.

     :Params:
       INPUTS:
                N - number of data points
               NW - time-halfbandwidth product
             Ktpr - number of tapers to be used (max value is 2*NW - 1)

       OUTPUTS:
         dpss.
             .E     discrete prolate spheroidal sequences (dpss), eigenvectors
             .V     dpss eigenvalues
             .N     number of data points
             .NW    time-halfbandwidth product
             .Ktpr  number of tapers to be used (max value is 2*NW - 1)

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/spd_mtm_dpss.pro)

## SPD_MTM_FINDPEAKS

[Previous Routine] [Next Routine] [List of Routines]

     :Description:
       Find significant enhancements in the gamma and F test.

     :Params:
       INPUTS:
            spec.ff - Fourier frequencies
           spec.raw - adaptive multitaper PSD
          spec.back - best PSD background among the probed ones
          spec.conf - confidence threshold values for the PSD (gamma test)
         spec.ftest - values of the F test
         spec.fconf - confidence threshold values for the F statistic (F test)
           par.conf - confidence thresholds percentages
          par.nfreq - number of frequencies
             par.df - frequency resolution (after padding), it corresponds to
                      the Rayleigh frequency for no padding (padding = 1)
             par.NW - time-halfbandwidth product
           par.fray - Rayleigh frequency: fray = 1/(npts*dt)
      ipar.peakproc - array[4] referring to 'gt', 'ft', 'gft', and 'gftm'
                      1 (0) peaks according to this procedure are (not) saved
       ipar.allpkwd - keyword /allpeakwidth is selected (1) or not (0)

       OUTPUTS:
           peak.
               .ff    Fourier frequencies
               .pkdf  for each peak selection method and confidence level
                      a value greater than zero at a specific frequency
                      indicate the occurence of a signal at that frequency:
                      peak.pkdf[#peakproc, #conf, #freq]
                      'gamma test' -> peak.pkdf[0,*,*] contains the badwidth
                      of the PSD enhancements at the identified frequencies
                      'F test', 'gft', and 'gftm' -> peak.pkdf[1:3,*,*]
                      is equal to par.df at the identified frequencies

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/spd_mtm_findpeaks.pro)

## SPD_MTM_FITMODEL

[Previous Routine] [Next Routine] [List of Routines]

     :Description:
       Evaluate the PSD background according to four different models:
       'wht'(0) white noise, [c]
       'pl'(1) power law, [c, beta]
       'ar1'(2) first order autoregressive process, [c, rho]
       'bpl'(3) bending power law, [c, beta, gamma, fb]

     :Params:
       INPUTS:
             spec.ff - Fourier frequencies
            spec.raw - adaptive multitaper PSD
            spec.dof - degree of freedom at each Fourier frequency
           spec.smth - smoothed PSD, spec.smth[#smth, *]:
                       #smth = 0->'raw'; 1->'med'; 2->'mlog', 3->'bin'; 4->'but'
           spec.fbin - binned frequencies for the bin-smoothed PSD
        spec.psmooth - percentage of the frequency range defining
                       the smoothing window spec.psmooth[#smth]
           par.nfreq - number of frequencies
            par.npts - time series number of points
              par.df - frequency resolution (after padding), it corresponds to
                       the Rayleigh frequency for no padding (padding = 1)
         par.datavar - variance of the time series
             par.fny - Nyquist frequency: fny = 1/(2*dt)
       ipar.specproc - array[#smth,#modl] smoothing and model combinations
                       1 (0) the smoothing + model combination is (not) probed

       OUTPUTS:
           spec.modl - PSD background models, spec.modl[#smth, #modl, *]:
                       #modl = 0->'wht'; 1->'pl'; 2->'ar1'; 3->'bpl'
           spec.frpr - model parameters resulting from the fitting procedures
                       spec.frpr[#smth, #modl, #par]:
                       #par = 0->'c'; 1->'rho' or 'beta'; 2->'gamma'; 3->'fb'

       COMMON VARIABLES:
         max_lklh
              psd_ml - PSD for the maximum likelihood model fit
               ff_ml - frequency vector for the maximum likelihood model fit
            alpha_ml - half degree of freedom of the PSD at each frequency
                 fny - Nyquist frequency
                itmp - indices of the frequency interval of interest
              n_itmp - number of frequencies in the interval of interest

     :Uses:
         spd_mtm_pl_fun.pro
         spd_mtm_pl_lglk.pro
         spd_mtm_pl_norm.pro
         spd_mtm_ar1_fun.pro
         spd_mtm_ar1_lglk.pro
         spd_mtm_ar1_norm.pro
         spd_mtm_bpl_fun.pro
         spd_mtm_bpl_lglk.pro
         spd_mtm_bpl_norm.pro

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/spd_mtm_fitmodel.pro)

## SPD_MTM_MAKEPLOT

[Previous Routine] [Next Routine] [List of Routines]

     :Description:
       Plot the results from spd_mtm.pro.

     :Params:
       INPUTS:
            data - 2 column data: [time vector, time series] = [x, y]
                   N.B. The average value of the data is removed by default.
                   If the data are not evenly sampled, the average time step
                   is used when it is lower than its standard deviation.
            spec - structure with results of the spectral analysis
                   (as defined in spd_mtm.pro)
            peak - identified signals (as defined in spd_mtm.pro)
             par - properties of the time series and parameters defining
                   the spectral analysis (as defined in spd_mtm.pro)
            ipar - indeces define the procedures for the identification
                   of PSD background and signals (as defined in spd_mtm.pro)
         x_label - label for x-axis or defined set of choice
                   (default choice 'Time')
         x_units - x variable unit of measures
         y_units - y variable unit of measures
         f_units - "frequency" variable unit of measures
          x_conv - conversion factor for the x variable
                   (N.B. IS UP TO THE USER TO BE CONSISTENT WITH X_UNITS)
          f_conv - conversion factor for the "frequency" variable
                   (N.B. IS UP TO THE USER TO BE CONSISTENT WITH F_UNITS)

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/spd_mtm_makeplot.pro)

## SPD_MTM_MODELGOF

[Previous Routine] [Next Routine] [List of Routines]

     :Description:
       Determine the best PSD background model according to a
       statistical criterium among the MERIT, AIC, and C_KS
       when comparing multiple smoothing+model pairs.

     :Params:
       INPUTS:
             spec.ff - Fourier frequencies
            spec.raw - adaptive multitaper PSD
            spec.dof - degree of freedom at each Fourier frequency
           spec.modl - PSD background models, spec.modl[#smth, #modl, *]:
                       #modl = 0->'wht'; 1->'pl'; 2->'ar1'; 3->'bpl'
           spec.frpr - model parameters resulting from the fitting procedures
                       spec.frpr[#smth, #modl, #par]:
                       #par = 0->'c'; 1->'rho' or 'beta'; 2->'gamma'; 3->'fb'
       ipar.specproc - array[#smth,#modl] smoothing and model combinations
                       1 (0) the smoothing + model combination is (not) probed
            ipar.gof - criterium used to select the PSD background

       OUTPUTS:
            spec.CKS - C_KS value for each probed model
            spec.AIC - AIC value for each probed model
          spec.MERIT - MERIT value for each probed model
        spec.indback - indices of the selected PSD background;
                       smoothing/model = spec.indback[0]/[1]
           spec.back - best PSD background among the probed model

       COMMON VARIABLES:
         max_lklh
              psd_ml - PSD for the maximum likelihood model fit
               ff_ml - frequency vector for the maximum likelihood model fit
            alpha_ml - half degree of freedom of the PSD at each frequency
                itmp - indices of the frequency interval of interest
              n_itmp - number of frequencies in the interval of interest

     :Uses:
         spd_mtm_cks.pro
         spd_mtm_wht_lglk.pro
         spd_mtm_pl_lglk.pro
         spd_mtm_ar1_lglk.pro
         spd_mtm_bpl_lglk.pro

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/spd_mtm_modelgof.pro)

## SPD_MTM_REGRE

[Previous Routine] [Next Routine] [List of Routines]

     :Description:
       Perform the harmonic analysis (F test) according to
       Thomson, D. J. (1982). "Spectrum estimation and harmonic analysis."
       Proceedings of the IEEE,70(9), 1055-1096. doi:10.1109/PROC.1982.12433

     :Params:
       INPUTS:
         spec.Ryk - real part of the eigenspectra
         spec.Iyk - imaginary part of the eigenspectra
         par.tprs - discrete prolate spheroidal sequences (dpss)
         par.Ktpr - number of tapers to be used (max value is 2*NW - 1)
        par.nfreq - number of frequencies

       OUTPUTS:
       spec.ftest - values of the F test
         spec.muf - complex amplitude at each Fourier frequency

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/spd_mtm_regre.pro)

## SPD_MTM_RESHAPE

[Previous Routine] [Next Routine] [List of Routines]

     :Description:
       Reshaping of the power spectral density as described in:
       Thomson, D. J. (1982). "Spectrum estimation and harmonic analysis."
       Proceedings of the IEEE,70(9), 1055-1096. doi:10.1109/PROC.1982.12433

     :Params:
       INPUTS:
            spec - structure with results of the spectral analysis
                   (as defined in spd_mtm.pro)
             par - properties of the time series and parameters defining
                   the spectral analysis (as defined in spd_mtm.pro)
            ipar - indeces define the procedures for the identification
                   of PSD background and signals (as defined in spd_mtm.pro)
            peak - identified signals (as defined in spd_mtm.pro)
          demean - time series with average value removed

       OUTPUTS:
         rshspec - provide the spec structures based on the reshaped PSD
         rshpeak - provide the peak structures based on the reshaped PSD

     :Uses:
         spd_mtm_findpeaks.pro
         spd_mtm_adapt.pro
         spd_mtm_smoothing.pro
         spd_mtm_fitmodel.pro
         spd_mtm_modelgof.pro
         spd_mtm_confthr.pro

     :File_comments:
       We perform the PSD reshaping removing the contribute of the identified
       monochromatic signals from the PSD along the entire frequency range.
       The code for the local reshaping of the PSD, that is in a frequency
       interval with the same width of the spectral window main lobe,
       is provided as comment lines.

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/spd_mtm_reshape.pro)

## SPD_MTM_SMOOTHING

[Previous Routine] [Next Routine] [List of Routines]

     :Description:
       Evaluate the smoothed PSD according to four different approach:
       'med'(1) running median
       'mlog'(2) running median (constant window on log(f) )
       'bin'(3) binned PSD
       'but'(4) low pass filtered PSD (Butterworth filter)

     :Params:
       INPUTS:
             spec.ff - Fourier frequencies
            spec.raw - adaptive multitaper PSD
            spec.dof - degree of freedom at each Fourier frequency
            par.fbeg - beginning frequency of the interval under analysis
            par.fend - ending frequency of the interval under analysis
            par.Ktpr - number of tapers (max value is 2*NW - 1)
         par.psmooth - value imposed in spd_mtm
           par.nfreq - number of frequencies
       ipar.specproc - array[#smth,#modl] smoothing and model combinations
                       1 (0) the smoothing + model combination is (not) probed

       OUTPUTS:
           spec.smth - smoothed PSD, spec.smth[#smth, *]:
                       #smth = 0->'raw'; 1->'med'; 2->'mlog', 3->'bin'; 4->'but'
           spec.fbin - binned frequencies for the bin-smoothed PSD
        spec.psmooth - percentage of the frequency range defining
                       the smoothing window spec.psmooth[#smth]

       COMMON VARIABLES:
         max_lklh
              psd_ml - PSD for the maximum likelihood model fit
               ff_ml - frequency vector for the maximum likelihood model fit
            alpha_ml - half degree of freedom of the PSD at each frequency
                 fny - Nyquist frequency
                itmp - indices of the frequency interval of interest
              n_itmp - number of frequencies in the interval of interest
         gammaj
            alpha_gj - half degree of freedom of the PSD at each frequency
                Ktpr - number of tapers (max value is 2*NW - 1)

     :Uses:
         spd_mtm_med_ks.pro
         spd_mtm_med.pro
         spd_mtm_mlog_ks.pro
         spd_mtm_mlog.pro
         spd_mtm_bin_ks.pro
         spd_mtm_bin.pro
         spd_mtm_but_ks.pro
         spd_mtm_but.pro

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/spd_mtm_smoothing.pro)

## SPD_MTM_SPEC

[Previous Routine] [Next Routine] [List of Routines]

     :Description:
       This procedure evaluates the Power Spectral Density.

     :Params:
       INPUTS:
         data          time series

          par.
             .npts     time series number of points
             .dt       average sampling time
             .npad     time series number of points after padding
             .nfreq    number of frequencies
             .df       frequency resolution (after padding), it corresponds to
                       the Rayleigh frequency for no padding (padding = 1)
             .psmooth  value imposed in spd_mtm
             .conf     confidence thresholds percentages
             .Ktpr     number of tapers (max value is 2*NW - 1)
             .tprs     discrete prolate spheroidal sequences (dpss)
             .v        dpss eigenvalues

         ipar.
             .hires    keyword /hires is selected (1) or not (0)
             .power    keyword /power is selected (1) or not (0)
             .specproc array[#smth,#modl] smoothing and model combinations
                       1 (0) the smoothing + model combination is (not) probed

       OUTPUTS:
         spec.
             .ff       Fourier frequencies
             .raw      adaptive multitaper PSD
             .dof      degree of freedom at each Fourier frequency
             .Ryk      real part of the eigenspectra spec.Ryk[#Ktpr,*]
             .Iyk      imaginary part of the eigenspectra spec.Iyk[#Ktpr,*]
             .poor_MTM flag for failed convergence of the adaptive MTM PSD

     :Uses:
         spd_mtm_adapt.pro

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/spd_mtm_spec.pro)

## TPLOT_SPD_MTM

[Previous Routine] [List of Routines]

     :Description:
       Perform a spectral analysis via the spd_mtm.pro procedure
       on a selected interval of a tplot variable.

     :Params:
         vname - tplot variable name
        trange - time range
      dec_step - decimation step: number of points between values to consider
            NW - time-halfbandwidth product, DEFAULT: NW = 3
          Ktpr - Number of tapers to apply, DEFAULT: Ktpr = 2*NW - 1
       padding - Amount of padding to apply (proportional to the original data)
                 DEFAULT: padding=1 so that the frequency step (df)
                 equals the Rayleigh frequency [1/(N*dt)]
          dpss - structure with tapers and eigenvalues
                 evaluated by spd_mtm_dpss.pro or previous call of spd_mtm.pro
          flim - limits in percentages of the frequency range to be analyzed,
                 DEFAULT: [df,fny-df], where fny is the Nyquist frequency [1/(2*dt)]
     smoothing - possible smoothing:
                 'raw'(0) no smoothing (DEFAULT)
                 'med'(1) running median
                 'mlog'(2) running median (constant window on log(f) )
                 'bin'(3) binned PSD
                 'but'(4) low pass filtered PSD (Butterworth filter)
                 'all'(9) use all the smoothing procedures
       psmooth - value between (0,0.5] defining the "smoothing window"
                 psmooth=2 -> search the optimum window according
                              to the Kolmogorov-Smirnov test (DEFAULT)
        model  - implemented models:
                 'wht'(0) white noise, [N]
                 'pl'(1) power law, [N, beta]
                 'ar1'(2) first order autoregressive process, [N, rho]
                 'bpl'(3) bending power law, [N, beta, gamma, fb] (DEFAULT)
                 'all'(9) use all the models
      procpeak - possible choices:
                 '' do not search for peaks
                 'gt'(0) gamma test only
                 'ft'(1) F test only
                 'gft'(2) gamma and F test (DEFAULT)
                 'gftm'(3) only max F value for each PSD enhancement
                 'all'(9) use all the selection procedures
          conf - confidence levels in increasing order
                 (DEFAULT is conf=[0.90,0.95,0.99])
          resh - value between (0,1.0], reshape the PSD from peaks identified
                 according to the 'gft' procedure at the "resh" confidence level
                 (DEFAULT: resh=0, do not perform the reshaping)
           gof - PSD background chosing criterium:
                 'MERIT' - MERIT function (DEFAULT)
                 'CKS' - Kolmogorov-Smirnov test
                 'AIC' - Akaike Information Criteria

       x_label - label for x-axis or defined set of choice
                 (default choice 'Time', see spd_mtm_makeplot.pro)
       x_units - x variable unit of measures
       y_units - y variable unit of measures
       f_units - "frequency" variable unit of measures
        x_conv - conversion factor for the x variable
                 (N.B. IS UP TO THE USER TO BE CONSISTENT WITH X_UNITS)
        f_conv - conversion factor for the "frequency" variable
                 (N.B. IS UP TO THE USER TO BE CONSISTENT WITH F_UNITS)

     OUTPUTS:
       spec.
          spec - structure with results of the spectral analysis
                 (as defined in spd_mtm.pro)
          peak - identified signals (as defined in spd_mtm.pro)
           par - properties of the time series and parameters defining
                 the spectral analysis (as defined in spd_mtm.pro)
          ipar - indeces define the procedures for the identification
                 of PSD background and signals (as defined in spd_mtm.pro)
       rshspec - provide the spec structures based on the reshaped PSD
                 (If it is not specified, the results are overwritten on spec)
       rshpeak - provide the peak structures based on the reshaped PSD
                 (If it is not specified, the results are overwritten on peak)

     :Keywords:
         /hires - spec.raw is the high-resolution PSD estimate
         /power - spec.raw is the integrated PSD: df*PSD [##^2]
         /quiet - do not display parameters and results
         /allpeakwidth - no constrains on the PSD enhancements width
         /pltsmooth - plot only the smoothed PSDs and stop
         /no_plot - do not plot the results
         /double - output in double precision

     :Uses:
         spd_mtm.pro

     :Author:
         Simone Di Matteo, Ph.D.
         8800 Greenbelt Rd
         Greenbelt, MD 20771 USA
         E-mail: simone.dimatteo@nasa.gov

(See general/science/spd_mtm/tplot_spd_mtm.pro)