CODEWITHSUNDEEP
pythonsignal-processingphysics
markov zain
markov zain

Reputation: 12915

Calculate energy for each frequency band around frequency F of interest in Python

I am newbie in signal processing, in this question, i want to ask how to obtain energy for each frequency band around interested frequency F. I have found a formula, but I dont know how to implement it in Python. This is the formula and my Fourier transform plot: enter image description here

x = np.linspace(0,5,100)
y = np.sin(2*np.pi*x)

## fourier transform
f = np.fft.fft(y)
## sample frequencies
freq = np.fft.fftfreq(len(y), d=x[1]-x[0])
plt.plot(freq, abs(f)**2) ## will show a peak at a frequency of 1 as it should.

enter image description here

Upvotes: 8

Views: 6969

Answers (2)

Wajdan Ali
Wajdan Ali

Reputation: 177

You are almost there as Mike pointed but here is a different approach which is simpler to understand.you can set a variable that holds the filtered signal and return a 1d array of Af, then apply the above formula which is quite simple (squared sum of these amplitudes)

Filter out the signals like this

from scipy.signal import butter, lfilter
def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a    
def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y

now assuming y is your original signal and you need energy of 5Hz component in signal , 

 #let fs = 250
 #let order = 5
 oneD_array_of_amps_of_fiveHz_component = butter_bandpass_filter(y, 4, 6, 250, 5)
 #calculate energy like this
 energy_of_fiveHz_comp = sum([x*2 for x in oneD_array_of_amps_of_fiveHz_component])

Upvotes: 4

maxymoo
maxymoo

Reputation: 36545

Using the signal processing module from Finding local maxima/minima with Numpy in a 1D numpy array you would do the following:

from scipy.signal import argrelextrema
import numpy as np

delta = 0.1
F = 1
X_delta = abs(freq - F) < delta
A_f_delta = abs(f[X_delta])**2
maximum_ix = argrelextrema(A_f, np.greater)
E_F_delta = np.sum(A_f[maximum_ix])/2
E_F_delta
2453.8235776144866

Upvotes: 0

Related Questions