4

I am using Python numpy's ftt.ftt() method to generate the fourier transform of a signal. However I want to calculate the bandpower over a range of frequencies. MATLAB has the method bandpower(x,fs,freqrange), I am trying to simulate specifically this syntax of the function. Source: https://www.mathworks.com/help/signal/ref/bandpower.html

It doesn't look like numpy has an equivalent function, but does anyone know a code snippet I can use to mimic bandpower(x,fs,freqrange)? It's not clear to me what exactly is going on behind the scenes in the function.

Note: If you know some non-Python pseudocode that would achieve the Matlab function, that would also be helpful.

HexTree
  • 148
  • 1
  • 7

2 Answers2

6

The following snippet for computing the power in the band [fmin, fmax] worked for me: 

import scipy 

def bandpower(x, fs, fmin, fmax):
    f, Pxx = scipy.signal.periodogram(x, fs=fs)
    ind_min = scipy.argmax(f > fmin) - 1
    ind_max = scipy.argmax(f > fmax) - 1
    return scipy.trapz(Pxx[ind_min: ind_max], f[ind_min: ind_max])
JNM
  • 108
  • 4
  • Why do you use default scaling = 'density' instead of scaling = 'spectrum' – John Sep 28 '21 at 09:07
  • That's a good point! This may be off by a constant factor depending on matlab's implementation. I just checked the matlab docs and they're not as explicit about units as scipy's. I no longer have a matlab license to verify. Are you able to compare matlab's output vs this scipy implementation to verify whether scaling = 'spectrum' or scaling = 'density' match what matlab does? – JNM Sep 29 '21 at 15:59
0
def bandpower(x, fs, fmin, fmax, time):
    # f, Pxx = scipy.signal.periodogram(x, fs=fs)
    f, Pxx = signal.welch(x, fs, nperseg=fs*time)
    ind_min = np.argmax(f > fmin) - 1
    ind_max = np.argmax(f > fmax) - 1
    return np.trapz(Pxx[ind_min: ind_max], f[ind_min: ind_max])
user12084390
  • 1
  • 1