Vector and RMS averaging in FFT

Question

I have a data array on which I have performed an FFT. This is the code that I have applied.

import numpy as np

# "data" is a column vector on which FFT needs to be performed
# N = No. of points in "data" 
# dt = time interval between two corresponding data points

FFT_data         = np.fft.fft(data)       # Complex values
FFT_data_real    = 2/N*abs(FFT_data)      # Absolute values

However, I went through following link: https://www.dsprelated.com/showarticle/1159.php

Here it says, to enhance the SNR we can apply "RMS-averaged FFT" and "Vector Averaged FFT".

Can somebody please let me know how to we go about doing these two methodologies in Python or is there any documentation/links to which we can refer ?

Answer 1

As your reference indicates:

If you take the square root of the average of the squares of your sample spectra, you are doing RMS Averaging. Another alternative is Vector Averaging in which you average the real and complex components separately.

Obviously to be able to perform either averaging you'd need to have more than a single data set to average. In your example code, you have a single column vector data. Let's assume you have multiple such column vectors arranged as a 2D NxM matrix, where N is the number of points per dataset and M is the number of datasets. Since the datasets are stored in columns, when computing the FFT you will need to specify the parameter axis=0 to compute the FFT along columns.

RMS-averaged FFT

As the name suggests, for this method you need to take the square-root of the mean of the squared amplitudes. Since the different sets are stored in columns, you'd need to do the average along the axis 1 (the other axis than the one used for the FFT).

FFT_data        = np.fft.fft(data, axis=0)  # Complex values
FFT_data_real   = 2/N*abs(FFT_data)         # Absolute values
rms_averaged    = np.sqrt(np.mean(FFT_data_real**2, axis=1))

Vector Averaged FFT

In this case you need to obtain the real and imaginary components of the FFT data, then compute the average on each separately:

FFT_data        = np.fft.fft(data, axis=0)  # Complex values
real_part_avg   = 2/N*np.mean(np.real(FFT_data),axis=1)
imag_part_avg   = 2/N*np.mean(np.imag(FFT_data),axis=1)
vector_averaged = np.abs(real_part_avg+1j*imag_part_avg)

Note that I've kept the 2/N scaling you had for the absolute values.

But what can I do if I really only have one dataset?

If that dataset happens to be stationary and sufficiently large then you could break down your dataset into smaller blocks. This can be done by reshaping your vector into an NxM matrix with the following:

data = data.reshape(N,M)
...

Then you could perform the averaging with either method.