8

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

Eric O. Lebigot
  • 91,433
  • 48
  • 218
  • 260
markov zain
  • 11,987
  • 13
  • 35
  • 39

2 Answers2

4

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])
Wajdan Ali
  • 179
  • 13
0

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
Community
  • 1
  • 1
maxymoo
  • 35,286
  • 11
  • 92
  • 119
  • I fail to see how this answers the question: the energy E_F in the formula sums the squared Fourier coefficients *around frequency F*, over a *frequency width* Delta. What this answer does is completely different: it restricts the sum to the exact (one-point) local maxima (there is no Delta, even, in this solution, so it cannot even work!). – Eric O. Lebigot Jun 08 '15 at 23:29
  • oy ok so i added F and delta into my answer, I think that it was implicit in my first answer though... it's not required to answer the question, since it just has a single frequency. – maxymoo Jun 09 '15 at 02:21
  • 1
    This is good, but why all the lines that start at `maximum_ix`: you answer the question at line `A_f_delta = …` (which is actually called E_F in the original question). Side note: with your argument about the single frequency justifying looking only at the peaks, you could go all the way and even sum all the squared amplitudes, since your `argrelextrema()` picks up the only non-negligible amplitudes. That would obviously be quite useless, in the general case, because of Parseval's theorem (since the result would be the same as doing the same on the original signal). – Eric O. Lebigot Jun 09 '15 at 02:43