This page was created by the IDL library routine
mk_html_help. For more information on
this routine, refer to the IDL Online Help Navigator
or type:
? mk_html_help
at the IDL command line prompt.
Last modified: Fri Sep 18 10:39:01 2020.
: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
x = 0.5*cos(2.0*!pi*t/8.0) + randomn(3,512,1) ; time series
data = [[t],[x]]
spd_mtm_mtm, data=data, NW=3, Ktpr=5, padding=1, dpss=dpss,
flim=[0,1], smoothing=’all’, psmooth=2,
model=’wht’, procpeak=’all’, conf=[0.90,0.95d],
/makeplot,
x_label=’Time’, y_units=’##’,
x_units=’s’, x_conv=1,
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 /Users/sdimatte/spd_mtm_v1/spd_mtm.pro)
: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 /Users/sdimatte/spd_mtm_v1/spd_mtm_adapt.pro)
:Description:
First order autoregressive (AR(1)) function values
at the given frequencies with the ar1_par parameters.
:Params:
INPUTS:
f - frequency vector
pl_par - power law parameters:
ar1_par[0] --> c, constant factor
ar1_par[1] --> rho, lag one autoregressive coefficient
fny - Nyquist frequency
OUTPUTS:
ar1fun - first order autoregressive function
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/models/ar1/spd_mtm_ar1_fun.pro)
:Description:
Evaluate the log-likelihood of
the first order autoregressive (AR(1)) PSD model.
:Params:
INPUTS:
ar1_par - first order autoregressive PSD parameters
OUTPUTS:
M - log-likelihood
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
:Uses:
spd_mtm_ar1_norm.pro
spd_mtm_ar1_fun.pro
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/models/ar1/spd_mtm_ar1_lglk.pro)
:Description:
Constant normalization factor for the AR(1) model estimated with
the maximum likelihood theoretical solution.
:Params:
INPUTS:
rho - lag-one autoregressive coefficient
ff - Fourier frequencies vector
psd - power spectral density
alpha - half degree of freedom of Smtm for each frequency
fny - Nyquist frequency
OUTPUTS:
c - constant factor
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/models/ar1/spd_mtm_ar1_norm.pro)
:Description:
Evaluate the binned PSD.
:Params:
INPUTS:
PSD - Power Spectral Density
ff - frequency vector
win - running window number of points
OUTPUTS:
binPSD - bin-smoothed PSD
fbin - binned frequencies
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/smoothings/bin/spd_mtm_bin.pro)
:Description:
Kolmogorov-Smirnov test for the smoothing procedure:
binned spectrum.
:Params:
INPUTS:
smo_par - percentages of the frequency range
defining the smoothing window
OUTPUTS:
KS - distance between the estimated and theoretical
cumulative distribution function of gamma_j
COMMON VARIABLES:
max_lklh
psd_ml - PSD for the maximum likelihood model fit
ff_ml - frequency vector for the maximum likelihood model fit
:Uses:
spd_mtm_bin.pro
spd_mtm_ks.pro
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/smoothings/bin/spd_mtm_bin_ks.pro)
:Description:
Bending Power Law (BPL) function values at the given frequencies
with the bpl_par parameters.
:Params:
INPUTS:
f - frequency vector
bpl_par - bending power law parameters:
bpl_par[0] --> c, constant factor
bpl_par[1] --> beta, spectral index
bpl_par[2] --> gamma, spectral index
bpl_par[3] --> fb, frequency break
OUTPUTS:
bplfun - bending power law function
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/models/bpl/spd_mtm_bpl_fun.pro)
:Description:
Evaluate the log-likelihood of
the bending power law (BPL) PSD model.
:Params:
INPUTS:
bpl_par - bending power law PSD parameters
OUTPUTS:
M - log-likelihood
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
:Uses:
spd_mtm_bpl_norm.pro
spd_mtm_bpl_fun.pro
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/models/bpl/spd_mtm_bpl_lglk.pro)
:Description:
Constant normalization factor for the BPL model estimated with
the maximum likelihood theoretical solution.
:Params:
INPUTS:
ff - Fourier frequencies vector
bpl_par - bending power law parameters:
bpl_par[0] --> beta, spectral index
bpl_par[1] --> gamma, spectral index
bpl_par[2] --> fb, frequency break
bet - spectral index
ff - Fourier frequencies vector
psd - power spectral density
alpha - half degree of freedom of Smtm for each frequency
OUTPUTS:
c - constant factor
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/models/bpl/spd_mtm_bpl_norm.pro)
:Description:
Low pass filter of the PSD with the
butterworth filter.
:Params:
INPUTS:
PSD - Power Spectral Density
ncutoff - number of the frequency cutoff [1,nfreq/2]
norder - order of the filter
OUTPUTS:
butPSD - but-smoothed PSD
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/smoothings/but/spd_mtm_but.pro)
:Description:
Kolmogorov-Smirnov test for the smoothing procedure:
butterworth filtered PSD.
:Params:
INPUTS:
smo_par - percentages of the frequency range
defining the smoothing window
OUTPUTS:
KS - distance between the estimated and theoretical
cumulative distribution function of gamma_j
COMMON VARIABLES:
max_lklh
psd_ml - PSD for the maximum likelihood model fit
ff_ml - frequency vector for the maximum likelihood model fit
:Uses:
spd_mtm_but.pro
spd_mtm_ks.pro
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/smoothings/but/spd_mtm_but_ks.pro)
:Description:
Kolmogorov-Smirnov significance.
:Params:
INPUTS:
psd - power spectral density
OUTPUTS:
CKS - confidence level
:Uses:
spd_mtm_ks.pro
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/pdf_functions/spd_mtm_cks.pro)
: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 /Users/sdimatte/spd_mtm_v1/spd_mtm_confthr.pro)
: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 /Users/sdimatte/spd_mtm_v1/spd_mtm_dbl2flt.pro)
: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 /Users/sdimatte/spd_mtm_v1/spd_mtm_dispar.pro)
: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 /Users/sdimatte/spd_mtm_v1/spd_mtm_dpss.pro)
: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 /Users/sdimatte/spd_mtm_v1/spd_mtm_findpeaks.pro)
: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 /Users/sdimatte/spd_mtm_v1/spd_mtm_fitmodel.pro)
:Description:
Return the cumulative distribution funtion (CDF) value at a given point
for the distribution of the random variable gamma.
:Params:
INPUTS:
x - gamma values
OUTPUTS:
gamj_cdf - CDF at the x values
COMMON VARIABLES:
gammaj
alpha_gj - half degree of freedom of the PSD at each frequency
Ktpr - number of tapers (max value is 2*NW - 1)
gamj_cdf_xx - CDF at fixed gamma values
xx - fixed gamma values (cover up to Ktpr=20)
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/pdf_functions/spd_mtm_gammaj_cdf.pro)
:Description:
Return the cutting value at a given confidence threshold pp
for the distribution of the random variable gamma.
:Params:
INPUTS:
x - gamma values
COMMON VARIABLES:
gammaj
alpha_gj - half degree of freedom of the PSD at each frequency
Ktpr - number of tapers (max value is 2*NW - 1)
pp - confidence threshold percentage
gamj_cdf_xx - CDF at fixed gamma values
xx - fixed gamma values (cover up to Ktpr=20)
:Uses:
spd_mtm_gammaj_cdf.pro
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/pdf_functions/spd_mtm_gammaj_cvf.pro)
:Description:
Kolmogorov-Smirnov test.
:Params:
INPUTS:
Sj - Power Spectral Density
OUTPUTS:
KS - distance between the estimated and theoretical
cumulative distribution function of gamma_j
COMMON VARIABLES:
max_lklh
psd_ml - PSD for the maximum likelihood model fit
:Uses:
spd_mtm_gammaj_cdf.pro
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/pdf_functions/spd_mtm_ks.pro)
: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 /Users/sdimatte/spd_mtm_v1/spd_mtm_makeplot.pro)
:Description:
Evaluate the running median of the PSD.
:Params:
INPUTS:
PSD - Power Spectral Density
win - running window number of points
OUTPUTS:
medPSD - med-smoothed PSD
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/smoothings/med/spd_mtm_med.pro)
:Description:
Kolmogorov-Smirnov test for the smoothing procedure:
running median.
:Params:
INPUTS:
smo_par - percentages of the frequency range
defining the smoothing window
OUTPUTS:
KS - distance between the estimated and theoretical
cumulative distribution function of gamma_j
COMMON VARIABLES:
max_lklh
psd_ml - PSD for the maximum likelihood model fit
ff_ml - frequency vector for the maximum likelihood model fit
:Uses:
spd_mtm_med.pro
spd_mtm_ks.pro
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/smoothings/med/spd_mtm_med_ks.pro)
:Description:
Evaluate the running median of the PSD
with a uniform moving window in the log-frequency space.
:Params:
INPUTS:
PSD - Power Spectral Density
ff - frequency vector
psmooth - percentage of the frequency range defining
the smoothing window
OUTPUTS:
mlogPSD - mlog-smoothed PSD
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/smoothings/mlog/spd_mtm_mlog.pro)
:Description:
Kolmogorov-Smirnov test for the smoothing procedure:
running median with uniform window in log-frequency space.
:Params:
INPUTS:
smo_par - percentages of the frequency range
defining the smoothing window
OUTPUTS:
KS - distance between the estimated and theoretical
cumulative distribution function of gamma_j
COMMON VARIABLES:
max_lklh
psd_ml - PSD for the maximum likelihood model fit
ff_ml - frequency vector for the maximum likelihood model fit
:Uses:
spd_mtm_mlog.pro
spd_mtm_ks.pro
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/smoothings/mlog/spd_mtm_mlog_ks.pro)
: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 /Users/sdimatte/spd_mtm_v1/spd_mtm_modelgof.pro)
:Description:
Power Law (PL) function values at the given frequencies
with the pl_par parameters.
:Params:
INPUTS:
f - frequency vector
pl_par - power law parameters:
pl_par[0] --> c, constant factor
pl_par[1] --> beta, spectral index
OUTPUTS:
plfun - power law function
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/models/pl/spd_mtm_pl_fun.pro)
:Description:
Evaluate the log-likelihood of
the power law (PL) PSD model.
:Params:
INPUTS:
pl_par - power law PSD parameters:
OUTPUTS:
M - log-likelihood
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
:Uses:
spd_mtm_pl_norm.pro
spd_mtm_pl_fun.pro
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/models/pl/spd_mtm_pl_lglk.pro)
:Description:
Constant normalization factor for the PL model estimated with
the maximum likelihood theoretical solution.
:Params:
INPUTS:
bet - spectral index
ff - Fourier frequencies vector
psd - power spectral density
alpha - half degree of freedom of Smtm for each frequency
OUTPUTS:
c - constant factor
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/models/pl/spd_mtm_pl_norm.pro)
: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 /Users/sdimatte/spd_mtm_v1/spd_mtm_regre.pro)
: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 /Users/sdimatte/spd_mtm_v1/spd_mtm_reshape.pro)
: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 /Users/sdimatte/spd_mtm_v1/spd_mtm_smoothing.pro)
: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 /Users/sdimatte/spd_mtm_v1/spd_mtm_spec.pro)
:Description:
Evaluate the log-likelihood of
the white noise (WHT) PSD model.
:Params:
INPUTS:
wht_par - power law PSD parameters:
OUTPUTS:
M - log-likelihood
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
:Author:
Simone Di Matteo, Ph.D.
8800 Greenbelt Rd
Greenbelt, MD 20771 USA
E-mail: simone.dimatteo@nasa.gov
(See /Users/sdimatte/spd_mtm_v1/models/pl/spd_mtm_wht_lglk.pro)