Multi-taper Spectral Analysis with spectrum in python

Question

I'm using the multi-taper analysis using the spectrum library on python (https://pyspectrum.readthedocs.io/en/latest/install.html), but I can't understand fully the amplitude of the output.

Here a piece of code for illustration:

from spectrum import *
N=500
dt=2*10**-3
# Creating a signal with 2 sinus waves.
x = np.linspace(0.0, N*dt, N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)

# classical FFT
yf = fft.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*dt), N//2)

# The multitapered method
NW=2.5
k=4
[tapers, eigen] = dpss(N, NW, k)
Sk_complex, weights, eigenvalues=pmtm(y, e=eigen, v=tapers, NFFT=500, show=False)

Sk = abs(Sk_complex)
Sk = np.mean(Sk * np.transpose(weights), axis=0)

# ploting both the results
plt.plot(xf,abs(yf[0:N//2])*dt*2)

plt.plot(xf,Sk[0:N//2])

Both the results are similar and find frequency peak at 50 and 80 Hz. The classical FFT finds as well the good amplitude (1 and 0.5) But the multi taper method do not find the proper amplitude. In this example it is around 5 times to important. Do anyone knows actually how to properly display the results ? thanks

Answer 1

From my understanding, there is a couple of factors that are at play here.

First, to get the multitaper estimate of the power spectrum density, you should compute like this:

Sk = abs(Sk_complex)**2
Sk = np.mean(Sk * np.transpose(weights), axis=0) * dt

I.e., you need to average the power spectrum, not the Fourier components.

Second, to get the power spectrum, you just need to divide the energy spectrum by N of your estimation using fft and multiply by dt as you did (and you need the **2 to get the power from the fourier components):

plt.plot(xf,abs(yf[0:N//2])**2 / N * dt)
plt.plot(xf,Sk[0:N//2])

Finally, what should be directly comparable, is not so much the amplitude in the power spectrum density, but the total power. You can look at:

print(np.sum(abs(yf[0:N//2])**2/N * dt), np.sum(Sk[0:N//2]))

Which match very closely.

So your whole code becomes:

from spectrum import *
N=500
dt=2*10**-3
# Creating a signal with 2 sinus waves.
x = np.linspace(0.0, N*dt, N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)

# classical FFT
yf = fft.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*dt), N//2)

# The multitapered method
NW=2.5
k=4
[tapers, eigen] = dpss(N, NW, k)
Sk_complex, weights, eigenvalues=pmtm(y, e=eigen, v=tapers, NFFT=N, show=False)

Sk = abs(Sk_complex)**2
Sk = np.mean(Sk * np.transpose(weights), axis=0) * dt

# ploting both results
plt.figure()
plt.plot(xf,abs(yf[0:N//2])**2 / N * dt)
plt.plot(xf,Sk[0:N//2])

# ploting both results in log scale
plt.semilogy(xf, abs(yf[0:N // 2]) ** 2 / N * dt)
plt.semilogy(xf, Sk[0:N // 2])

# comparing total power
print(np.sum(abs(yf[0:N//2])**2 / N * dt), np.sum(Sk[0:N//2]))