I noticed that the final spectrum estimate returned by a call to mtspec() is scaled appropriately, i.e., the Riemann sum of the spectrum (np.sum(spec)*df) is equal to the variance of the data. To get this result there is a factor of two (non-negative frequencies returned only) and also scaling by dt or df (lines 1819-1832 of mtspec.f90).
However, the eigenspectra are neither scaled by 2, although only the non-negative frequencies are returned, nor scaled by dt (or df). The average of the eigenspectra should be close to the returned adaptively weighted multitaper estimate, especially if using Gaussian white noise, but they will be off as a result (a factor of approximately 2*dt).
The eigencoefficients also do not appear to be scaled by dt or df (no 2 needed as the full frequency band is returned).
This results in the bjk quantity (jackspec.f90, Line 123) being off by a factor of (K-1)ln(2) I think?
If this is the desired / expected behaviour, then disregard this please.
If not, one potential way to address this is to scale the DPSS's by sqrt(dt) as soon as they are calculated - this is the method used in the R package "multitaper", which is also based on Dave Thomson's original Fortran code (https://github.com/wesleyburr/multitaper/blob/master/R/multitaper.R#L250).
I noticed that the final spectrum estimate returned by a call to mtspec() is scaled appropriately, i.e., the Riemann sum of the spectrum (np.sum(spec)*df) is equal to the variance of the data. To get this result there is a factor of two (non-negative frequencies returned only) and also scaling by dt or df (lines 1819-1832 of mtspec.f90).
However, the eigenspectra are neither scaled by 2, although only the non-negative frequencies are returned, nor scaled by dt (or df). The average of the eigenspectra should be close to the returned adaptively weighted multitaper estimate, especially if using Gaussian white noise, but they will be off as a result (a factor of approximately 2*dt).
The eigencoefficients also do not appear to be scaled by dt or df (no 2 needed as the full frequency band is returned).
This results in the bjk quantity (jackspec.f90, Line 123) being off by a factor of (K-1)ln(2) I think?
If this is the desired / expected behaviour, then disregard this please. If not, one potential way to address this is to scale the DPSS's by sqrt(dt) as soon as they are calculated - this is the method used in the R package "multitaper", which is also based on Dave Thomson's original Fortran code (https://github.com/wesleyburr/multitaper/blob/master/R/multitaper.R#L250).