Python numpy equivalent of bandpower() from MATLAB

Question

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.

Are you trying to simulate all of the features or just with one input argument syntax? — Divakar
– Divakar, Commented Jun 14, 2017 at 14:33
No idea mate, I've completely forgotten what this all about, haha. — HexTree
– HexTree, Commented Mar 12, 2020 at 13:05

JNM · Accepted Answer · 2017-06-26 23:37:41Z

7

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])

answered Jun 26, 2017 at 23:37

JNM

1184 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

John Over a year ago

Why do you use default scaling = 'density' instead of scaling = 'spectrum'

JNM Over a year ago

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?

Al Grant Over a year ago

@JNM would you mind adding some more comments to your code? I would like to understand what x, fs, Pxx, ind_min and ind_max are and any other explanations you care to add.

user12084390 · Accepted Answer · 2021-12-06 12:21:01Z

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])

answered Dec 6, 2021 at 12:21

user12084390

11 bronze badge

2 Comments

user12084390 Over a year ago

by according to the last answer from John Maidens

Al Grant Over a year ago

Can you add some comments to your code? what does time represent? If I am feeding in a np.sample of IQ then what are the freq values?

Collectives™ on Stack Overflow

Python numpy equivalent of bandpower() from MATLAB

2 Answers 2

3 Comments

2 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

3 Comments

2 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related