I have 2 snippets of code that should produce the same result. In the first one, I consider a 2500 sample signal on 64 channels sampled at 500 Hz; and in the second one, I cut down this signal into 5 chunks of 500 samples each.
# Snippet 1:
data = get_data() # data shape is 64x2500
# Snippet 2:
data = get_data().reshape(64, 5, 500) # data shape is 64x5x500
Where get_data() is just a placeholder for this example, feel free to generate a random sample array to use with both pieces of code. In practice, to test my code I duplicated a sample signal with copy.deepcopy().
My goal is to extract the frequency power on 2 bands by averaging on the band and across all 64 channels.
- alpha: (8, 13) Hz
- delta: (1, 4) Hz
To do so, I apply a fast Fourier transform through numpy: (The code below should work on both snippets with both array dimensions.)
import numpy as np
alpha = (8, 13)
delta = (1, 4)
fs = 500. # sampling frequency
frequencies = np.fft.rfftfreq(data.shape[-1], 1.0/fs)
alpha_band = np.where(np.logical_and(frequencies>=alpha[0],
frequencies<=alpha[1]))[0]
delta_band = np.where(np.logical_and(frequencies>=delta[0],
frequencies<=delta[1]))[0]
fftval = np.abs(np.fft.rfft(data, axis=-1) / fs)
The absolute values are then squared. Moreover, each of the 64 channels also has a weight applied to it via multiplication. To test the code, I set all weights to 1.
weights = np.ones(shape=(64,))
alpha = np.average(np.multiply(np.square(fftval[..., alpha_band]).T, weights))
delta = np.average(np.multiply(np.square(fftval[..., delta_band]).T, weights))
My understanding was that averaging across chunks should not change the result. Yet, here are the 2 outputs I got on the same data:
# Snippet 1, 64x2500
alpha
Out: 5.294149596657551e-13
delta
Out: 9.372696436349552e-13
alpha/delta
Out: 0.564848081084284
# Snippet 2: 64x5x500
alpha
Out: 6.326672955916193e-12
delta
Out: 7.584706602278469e-13
alpha/delta
Out: 8.34135489699185
Anyone knows what I am doing wrong here and why both results are completely different? What did I miss?