Thursday, 31 December 2020

Difference between similar FFT code on complete signal or averaged on chunks

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?



from Difference between similar FFT code on complete signal or averaged on chunks

No comments:

Post a Comment