Unable to find the right energy using scipy.signal.welch

Question

For a given discrete time signal x(t) with an interval dt (which is equal to 1/fs , fs is the sampling frequency), the energy is equal to:

 E[x(t)] = sum(abs(x)**2.0)/fs 

Then I do the DFT x(t) :

 x_tf = np.fft.fftshift( np.fft.fft( x ) ) / ( fs * ( 2.0 * np.pi ) ** 0.5 ) 

and again calculate the energy:

 E[x_tf] = sum( abs( x_tf ) ** 2.0 ) * fs * 2 * np.pi / N 

(here the factor fs*2*np.pi/N = ripple distance dk , the fftfreq documentation gives more detailed information about the distance in the frequency domain), I have the same energy:

 E[x(t)] = E[x_tf] 

BUT ... when I calculate the spectral power density x(t) using scipy.signal.welch , I cannot find the energy I need. scipy.signal.welch returns the vector of frequencies f and energy Pxx (or energy per frequency, depending on which scaling we enter the arguments scipy.signal.welch ).

How can I find the same energy as E[x(t)] or E[x_tf] using Pxx ? I tried to calculate:

 E_psd = sum(Pxx_den) / nperseg 

where nperseg is the length of each segment of the Welch algorithm, factors such as fs and np.sqrt(2*np.pi) will be canceled and drag E [x (t)] with nperseg , but without any success (orders less E[x(t)] )

I used the following code to generate my signal:

 #Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by 0.001 V**2/Hz of white noise sampled at 10 kHz. fs = 10e3 #sampling rate, dt = 1/fs N = 1e5 amp = 2*np.sqrt(2) freq = 1234.0 noise_power = 0.001 * fs / 2 time = np.arange(N) / fs x = amp*np.sin(2*np.pi*freq*time) x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape) 

and I did the following to get the power spectral density:

 f, Pxx_den = signal.welch(x, fs ) 

Answer 1

The resolution of this apparent discrepancy lies in a thorough understanding and application.

• continuous with respect to discrete Fourier transforms and
• energy, power and spectral power density of a given signal

I also struggled with this exact question, so I will try to be as explicit as possible in the discussion below.

Discrete Fourier Transform (DFT)

A continuous signal x (t) satisfying some integrability conditions has the Fourier transform X (f). However, when working with a discrete signal x [n], it is often customary to work with a discrete time Fourier transform (DTFT). I will denote the DTFT as X_ {dt} (f), where dt is equal to the time interval between adjacent samples. The key to answering your question requires that you acknowledge that the DTFT does not match the corresponding Fourier transform! In fact, the two are related as

X_ {dt} (f) = (1 / dt) * X (f)

Furthermore, the discrete Fourier transform (DFT) is simply a discrete DTFT sample. Of course, DFT is what Python returns when using np.fft.fft(...) . So your calculated DFT is not equal to the Fourier transform!

Power spectral density

scipy.signal.welch(..., scaling='density', ...) returns the power spectral density (PSD) estimate of the discrete signal x [n]. A full discussion of PSD is slightly beyond the scope of this publication, but for a simple periodic signal (for example, in your example) PSD S_ {xx} (f) is given as

S_ {xx} = | X (f) | ^ 2 / T

where | X (f) | represents the Fourier transform of the signal, and T is the total duration (in time) of the signal (if your signal x (t) instead was a random process, we would have to take the average of many implementations of the system ...). The total power in the signal is simply the integral of S_ {xx} over the system bandwidth. Using your code above, we can write

 import scipy.signal # Estimate PSD `S_xx_welch` at discrete frequencies `f_welch` f_welch, S_xx_welch = scipy.signal.welch(x, fs=fs) # Integrate PSD over spectral bandwidth # to obtain signal power `P_welch` df_welch = f_welch[1] - f_welch[0] P_welch = np.sum(S_xx_welch) * df_welch 

To contact your calculations np.fft.fft(...) (which return the DFT), we must use the information from the previous section, namely that

X [k] = X_ {dt} (f_k) = (1 / dt) * X (f_k)

Thus, in order to calculate the power spectral density (or apparent power) from the FFT calculations, we need to recognize that

S_ {xx} = | X [k] | ^ 2 * (dt ^ 2) / T

 # Compute DFT Xk = np.fft.fft(x) # Compute corresponding frequencies dt = time[1] - time[0] f_fft = np.fft.fftfreq(len(x), d=dt) # Estimate PSD `S_xx_fft` at discrete frequencies `f_fft` T = time[-1] - time[0] S_xx_fft = ((np.abs(Xk) * dt) ** 2) / T # Integrate PSD over spectral bandwidth to obtain signal power `P_fft` df_fft = f_fft[1] - f_fft[0] P_fft = np.sum(S_xx_fft) * df_fft 

Your values ​​for P_welch and P_fft should be very close to each other, and also close to the expected signal strength, which can be calculated as

 # Power in sinusoidal signal is simply squared RMS, and # the RMS of a sinusoid is the amplitude divided by sqrt(2). # Thus, the sinusoidal contribution to expected power is P_exp = (amp / np.sqrt(2)) ** 2 # For white noise, as is considered in this example, # the noise is simply the noise PSD (a constant) # times the system bandwidth. This was already # computed in the problem statement and is given # as `noise_power`. Simply add to `P_exp` to get # total expected signal power. P_exp += noise_power 

Note: P_welch and P_fft will not be exactly equal and probably will not even be equal to the exact digits. This is due to the fact that there are random errors associated with the estimation of the power spectral density. To reduce such errors, the Welch method splits your signal into several segments (the size of which is controlled by the nperseg ), calculates the PSD of each segment and averages the PSD to get a better estimate of the PSD signal (the more segments are averaged, the resulting random error). The FFT method, in fact, is equivalent only to calculating and averaging over one large segment. Thus, we expect some differences between P_welch and P_fft , but we should expect P_welch be more accurate.

Signal energy

As you said, the signal energy can be obtained from a discrete version of the Parseval theorem

 # Energy obtained via "integrating" over time E = np.sum(x ** 2) # Energy obtained via "integrating" DFT components over frequency. # The fact that `E` = `E_fft` is the statement of # the discrete version of Parseval theorem. N = len(x) E_fft = np.sum(np.abs(Xk) ** 2) / N 

We would now like to understand how S_xx_welch calculated above via scipy.signal.welch(...) refers to the total energy E in the signal. Above, S_xx_fft = ((np.abs(Xk) * dt) ** 2) / T Reordering the terms in this expression, we see that np.abs(Xk) ** 2 = (T / (dt ** 2)) * S_xx_fft . Moreover,

We know from above that np.sum(S_xx_fft) = P_fft / df_fft and that P_fft and P_welch approximately equal. Next, P_welch = np.sum(S_xx_welch) / df_welch , so we get

np.sum(S_xx_fft) = (df_welch / df_fft) * np.sum(S_xx_welch)

Next, S_xx_fft = ((np.abs(Xk) * dt) ** 2) / T Substituting S_xx_fft into the above equation and rearranging the terms, we arrive at

np.sum(np.abs(Xk) ** 2) = (T / (dt ** 2)) * (df_welch / df_fft) * np.sum(S_xx_welch)

The left-hand side (LHS) in the above equation should now be very close to the expression for the total energy in the signal calculated from the DFT components. Now notice that T / dt = N , where N is the number of sample points in your signal. Separating N , we now have LHS, which by definition is equal to E_fft calculated above. Thus, we can get the full energy in the signal from the Welch PSD through

 # Signal energy from Welch PSD E_welch = (1. / dt) * (df_welch / df_fft) * np.sum(S_xx_welch) 

E , E_fft and E_welch should be very close in value :) As discussed at the end of the previous section, we expect slight differences between E_welch compared to E and E_fft , but this is due to the fact that the values ​​obtained from the Welch method reduced the random error (i.e. more precisely).