FFT of signal data with windowing, overlapping and averaging

Question

This is my first ever question here so the help is really appreciated.

I am performing FFT on a signal. I want to perform windowing, 50% overlapping and averaging to the signal. There is a function scipy.signal.welch to perform this automatically but the output is in power spectral density. I want the output in magnitude and phase shift both, but from power spectral density only magnitude is achievable. Is there a way to compute phase shift from power spectral density or a simple way to do this analysis in the form of FFT rather than in power spectral density?

I know how to apply windowing in python but I do not know how to do overlapping and averaging manually.

Below is my code:

import numpy as np
from numpy.fft import fft, ifft, fftshift, fftfreq
import pandas as pd
import matplotlib.pyplot as plt
from scipy import signal
import scipy.fft

data = pd.read_csv('lucid_1p34g_1024fps_5mins.csv')

ref = data.loc[:,"Input 0"]
sensor1x = data.loc[:,"Input 1"]
sensor1y = data.loc[:,"Input 2"]
sensor1z = data.loc[:,"Input 3"]

fs = 1024
blockSize = 1024

f, Pxx = signal.welch(sensor1z, 1024, window='hann', nperseg=blockSize, 
                     noverlap=512)

plt.plot(f, Pxx)                # power spectral density plot
plt.show()

"""Manual Calculation"""

N = len(sensor1z)
n = np.arange(N)
T = N/fs
freq = n/T

window = np.hanning(N)
f1z = fft(sensor1z)                       #fft transform of input 3
plt.plot(freq, np.abs(f1z))
plt.show()

Answer 1

I want to perform windowing, 50% overlapping and averaging to the signal

This makes sense for magnitude, but not for phase.

Is there a way to compute phase shift from power spectral density

No. The PSD is computed by averaging magnitude spectra together, so there is no phase information.

or a simple way to do this analysis in the form of FFT rather than in power spectral density

Yes, what you want is simply the phase of f1z:

plt.plot(freq, np.angle(f1z))

---

Depending on what information you want from the phase, you might want to unwrap / normalize the phase:

plt.plot(freq, np.unwrap(np.angle(f1z)))
# or
plt.plot(freq, np.unwrap(np.angle(f1z)*180/pi))

EDIT per the OP's request

What you are asking is a coding question: you already know how to perform fft and get magnitude and phase. What you want is doing this on overlapping segments of your signal, then averaging.

Again, not going to write the code for you but here would be an approach:

1. Define segment length $N$. To do this, figure out what frequency resolution $d_f$ you're happy with, and compute $N = \text{ceil}(f_s/d_f)$

2. Define overlap $R$: you can experiment with this, common ones are $R = 1/2$ or $R = 1/4$ for example.

3. Define a window function of length $N$. A Hann window for example. Or if you just want a rectangular window, go on to step 4 and dis-regard the np.multiply operation.

4. Now here comes the coding part:

   • do fft on np.multiply(data(0:N-1), window(0:N-1)).
   • extract magnitude and phase and store both somewhere
   • advance by R: do fft on np.multiply(data(R:N-1+R), window(0:N-1))
   • extract magnitude and phase and store both
   • advance by R: do fft on np.multiply(data(2R:N-1+2R), window(0:N-1))
   • extract magnitude and phase and store both
   • etc

   This can be done in a loop. Once you have all your magnitude and phase arrays, just average them together.

FYI, this is just a naive implementation of the Short Time Fourier Transform that @EricCanton mentions in his answer.

Answer 2

Are you looking for a short-term Fourier transform? This computes an FFT in a sliding window, giving an array of complex numbers from which you can extract magnitude and angle/phase. https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.stft.html

import numpy as np
import pandas as pd
from scipy import signal

data = pd.read_csv('lucid_1p34g_1024fps_5mins.csv')

ref = data.loc[:,"Input 0"]
sensor1x = data.loc[:,"Input 1"]
sensor1y = data.loc[:,"Input 2"]
sensor1z = data.loc[:,"Input 3"]

fs = 1024
blockSize = 1024

freq, ts, Zxx = signal.stft(
    sensor1z,
    fs=fs,
    nperseg=blockSize,
    noverlap=512,
    window='hann'
)

# Zxx is a numpy array of complex numbers
# get magnitude
Zxx_mag = np.abs(Zxx)
Zxx_phase = np.angle(Zxx)