0
$\begingroup$

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()
$\endgroup$
8
  • $\begingroup$ What do you mean by "phase shift"? That terms is usually associated with systems, not signals. The phase of a signal is not particularly easy to define, so you need to be specific about what exactly you want to calculate. $\endgroup$
    – Hilmar
    Nov 10, 2022 at 14:51
  • $\begingroup$ I assumed he just meant the phase response, but I could be wrong... $\endgroup$
    – Jdip
    Nov 10, 2022 at 15:10
  • 1
    $\begingroup$ "Response" implies system to me (like in impulse response). Relative phase between two signals is very meaningful. Absolute phase of a signal not so much, especially if you segment the signal. $\endgroup$
    – Hilmar
    Nov 10, 2022 at 16:23
  • $\begingroup$ Right, I mean absolute phase, not phase response. I assumed because of "I am performing FFT on a signal". And agreed, as mentioned at the top of my answer, no sense averaging... $\endgroup$
    – Jdip
    Nov 10, 2022 at 16:55
  • $\begingroup$ Are you looking for a short-term Fourier transform? docs.scipy.org/doc/scipy/reference/generated/… $\endgroup$ Nov 10, 2022 at 17:01

2 Answers 2

0
$\begingroup$

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.

$\endgroup$
7
  • $\begingroup$ How do I perform overlapping and averaging on f1z? $\endgroup$ Nov 10, 2022 at 13:36
  • $\begingroup$ "This makes sense for magnitude response, but not for phase response". You can use welch's method plt.plot(f, Pxx) for the magnitude, and angle(f1z) for the phase. $\endgroup$
    – Jdip
    Nov 10, 2022 at 13:38
  • $\begingroup$ You're asking now for custom code. I'm not going to write custom code for you that does overlapping and averaging, since that's the point of having functions such as welch. (Also, you can find that easily on google ;) ). If you want to write custom code from scratch, first study what welch's method is, try to write it yourself, and compare the output of your custom implementation vs the output of the built-in welch function. If you get stuck somewhere, come back and ask for help! $\endgroup$
    – Jdip
    Nov 10, 2022 at 13:44
  • $\begingroup$ But as far as your original question(s), those have been answered. Good luck! $\endgroup$
    – Jdip
    Nov 10, 2022 at 13:48
  • $\begingroup$ The welch method gives the output in form of power spectral density. I want the output in the form of magnitude like the one obtained by numpy.fft, but the problem I am facing is to perform overlapping and averaging. $\endgroup$ Nov 11, 2022 at 7:41
0
$\begingroup$

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)
$\endgroup$
2
  • $\begingroup$ Function STFT does not perform overlapping and averaging which I require for the output. $\endgroup$ Nov 11, 2022 at 7:42
  • $\begingroup$ it performs overlapping. You'd only have to do the averaging manually, which is a basic coding problem, but it's not the purpose of this website to write custom code for you. Zxx_mag and Zxx_phase have everything you need for what you're trying to achieve. $\endgroup$
    – Jdip
    Nov 11, 2022 at 9:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.