Understanding the frequency axis of multitaper method in Python

Question

I am using the spectrum package in Python, using the multitaper method to calculate the power spectral density (PSD) (see http://pyspectrum.readthedocs.io/en/latest/_modules/spectrum/mtm.html) of a sample of 60 magnetometer local magnetic north component readings, sampled every minute, measured in nanotesla.

My aim is to look at the PSD in the 1.5-6mHz band only, since this is the only band important for my research.

Running the function with the adaptive method (default NFFT=256, NW=1.4, default k) I get the following plot

My understanding is that the y axis has units (nT)^2/Hz which is fine, but how do I get the specific mHz frequency band of interest in the x axis? I initially thought that the x axis was in Hz, but for whatever NFFT value i choose the graph looks the same but just spans from 0 to the NFFT value in the x-axis.

Does anybody know how I go about finding the power in the 1.5 - 6 mHz range (or if I'm going about the wrong way entirely?

Thanks!

If it helps here is the code for what the function is doing (for other functions that it refers to please see the link I posted above, I have just posted this particular function since it gives more insight into what's going on (using np.FFT.FFT for example):

def pmtm(x, NW=None, k=None, NFFT=None, e=None, v=None, method='adapt', show=False):
"""Multitapering spectral estimation

:param array x: the data
:param float NW: The time half bandwidth parameter (typical values are
    2.5,3,3.5,4). Must be provided otherwise the tapering windows and
    eigen values (outputs of dpss) must be provided
:param int k: uses the first k Slepian sequences. If *k* is not provided,
    *k* is set to *NW*2*.
:param NW:
:param e: the window concentrations (eigenvalues)
:param v: the matrix containing the tapering windows
:param str method: set how the eigenvalues are used. Must be
    in ['unity', 'adapt', 'eigen']
:param bool show: plot results
:return: Sk (complex), weights, eigenvalues

Usually in spectral estimation the mean to reduce bias is to use tapering window.
In order to reduce variance we need to average different spectrum. The problem
is that we have only one set of data. Thus we need to decompose a set into several
segments. Such method are well-known: simple daniell's periodogram, Welch's method
and so on. The drawback of such methods is a loss of resolution since the segments
used to compute the spectrum are smaller than the data set.
The interest of multitapering method is to keep a good resolution while reducing
bias and variance.

How does it work? First we compute different simple periodogram with the whole data
set (to keep good resolution) but each periodgram is computed with a different
tapering windows. Then, we average all these spectrum. To avoid redundancy and bias
due to the tapers mtm use special tapers.

.. plot::
    :width: 80%
    :include-source:

    from spectrum import data_cosine, dpss, pmtm

    data = data_cosine(N=2048, A=0.1, sampling=1024, freq=200)
    # If you already have the DPSS windows
    [tapers, eigen] = dpss(2048, 2.5, 4)
    res = pmtm(data, e=eigen, v=tapers, show=False)
    # You do not need to compute the DPSS before end
    res = pmtm(data, NW=2.5, show=False)
    res = pmtm(data, NW=2.5, k=4, show=True)


.. versionchanged:: 0.6.2

    APN modified method to return each Sk as complex values, the eigenvalues
    and the weights

"""
assert method in ['adapt','eigen','unity']

N = len(x)

# if dpss not provided, compute them
if e is None and v is None:
    if NW is not None:
        [tapers, eigenvalues] = dpss(N, NW, k=k)
    else:
        raise ValueError("NW must be provided (e.g. 2.5, 3, 3.5, 4")
elif e is not None and v is not None:
    eigenvalues = e[:]
    tapers = v[:]
else:
    raise ValueError("if e provided, v must be provided as well and viceversa.")
nwin = len(eigenvalues) # length of the eigen values vector to be used later

# set the NFFT
if NFFT==None:
    NFFT = max(256, 2**nextpow2(N))

Sk_complex = np.fft.fft(np.multiply(tapers.transpose(), x), NFFT)
Sk = abs(Sk_complex)**2

# si nfft smaller thqn N, cut otherwise add wero.
# compute
if method in ['eigen', 'unity']:
    if method == 'unity':
        weights = np.ones((nwin, 1))
    elif method == 'eigen':
        # The S_k spectrum can be weighted by the eigenvalues, as in Park et al.
        weights = np.array([_x/float(i+1) for i,_x in enumerate(eigenvalues)])
        weights = weights.reshape(nwin,1)

elif method == 'adapt':
    # This version uses the equations from [2] (P&W pp 368-370).

    # Wrap the data modulo nfft if N > nfft
    sig2 = np.dot(x, x) / float(N)
    Sk = abs(np.fft.fft(np.multiply(tapers.transpose(), x), NFFT))**2
    Sk = Sk.transpose()
    S = (Sk[:,0] + Sk[:,1]) / 2    # Initial spectrum estimate
    S = S.reshape(NFFT, 1)
    Stemp = np.zeros((NFFT,1))
    S1 = np.zeros((NFFT,1))
    # Set tolerance for acceptance of spectral estimate:
    tol = 0.0005 * sig2 / float(NFFT)
    i = 0
    a = sig2 * (1 - eigenvalues)

    # converges very quickly but for safety; set i<100
    while sum(np.abs(S-S1))/NFFT > tol and i<100:
        i = i + 1
        # calculate weights
        b1 = np.multiply(S, np.ones((1,nwin)))
        b2 = np.multiply(S,eigenvalues.transpose()) + np.ones((NFFT,1))*a.transpose()
        b = b1/b2

        # calculate new spectral estimate
        wk=(b**2)*(np.ones((NFFT,1))*eigenvalues.transpose())
        S1 = sum(wk.transpose()*Sk.transpose())/ sum(wk.transpose())
        S1 = S1.reshape(NFFT, 1)
        Stemp = S1
        S1 = S
        S = Stemp   # swap S and S1
    weights=wk

if show is True:
    from pylab import semilogy
    if method == "adapt":
        Sk = np.mean(Sk * weights, axis=1)
    else:
        Sk = np.mean(Sk * weights, axis=0)
    semilogy(Sk)

return Sk_complex, weights, eigenvalues

Answer 1

If it's properly sampled, a discrete-time signal sampled at 1/60 Hz represents frequencies from -1/120 Hz through 1/120 Hz: https://en.wikipedia.org/wiki/Nyquist_frequency

If the samples are all real-valued, then the negative and positive frequency components are the same.

Frequency in discrete-time signals is circular/modular, i.e., in your signal, frequencies f and f + 1/60 Hz are the same.

In your PSD result, the x-axis runs over this range, starting at 0 Hz, up to 1/120 Hz, which is the same as -1/120Hz, and continuing up to just before 0 Hz again.

The frequency represented by bin x is f = x / (NFFT*60 Hz), and you can ignore all x >= NFFT/2.

If you're only interested in a particular band, then you can just pick the values you want and discard the rest, choosing NFFT to determine the resolution.