0
$\begingroup$

The power spectral density St of a signal u may be computed as the product of the FFT of the signal, u_fft with its complex conjugate u_fft_c. In Python, this would be written as:

import numpy as np

u = # Some numpy array containing signal
u_fft = np.fft.rfft(u-np.nanmean(u))

St = np.multiply(u_fft, np.conj(u_fft))

However, the FFT definition in Numpy requires the multiplication of the result with a factor of 1/N, where N=u.size in order to have an energetically consistent transformation between u and its FFT. This leads to the corrected definition of the PSD using numpy's fft:

St = np.multiply(u_fft, np.conj(u_fft))
St = np.divide(St, u.size)

On the other hand, Scipy's function signal.welch computes the PSD directly from input u:

from spicy.signal import welch

freqs_st, St_welch = welch(u-np.nanmean(u), 
     return_onesided=True, nperseg=seg_size, axis=0)

The resulting PSD, St_welch, is obtained by performing several FFTs in segments of the array u with size seg_size. Thus, my question is:

Should St_welch be multiplied by a factor of 1/seg_size to give an energetically consistent PSD? Should it be multiplied by 1/N? Should it not be multiplied at all?

PD: Comparison by performing both operations on a signal is not straightforward, since the Welch method also introduces smoothing of the signal and changes the display in the frequency domain.

Information on the necessity of the prefactor when using numpy.fft :

Journal article on the matter

$\endgroup$
0
$\begingroup$

I would like to post here the link to the accepted answer that I obtained to the question in StackOverflow, as I think it is relevant to this forum.

The original answer by francis is here: post in StackOverflow

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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