Application of Parsevals theorem with Welch method For my project I had to manually code the welch method using the code below. Pretty much it involves finding the spectral density through fft and incorporating windowing and segmenting with overlap.

I need to know if my PSD is correct so I am checking through 2 methods.


root of area under the PSD should be the RMS (E[x^2]) and since I have taken away my mean from the data, this should equal the standard deviation of the data
Checking if Parseval's theorem holds.
There are two main issues, I am having trouble with. Firstly I am getting a 0Hz peak even after removing the mean at the start. The only way I can remove the 0Hz peak is to take the mean away after I have segmented and multiplied the window. (#### Detrend line ####)


With the peak, point 1 holds up well but without the peak the standard deviation is slightly different from the RMS but thats not a huge concern.

My second issue is that I cant get the Parseval's theorem to hold. Is there something I'm doing wrong? Is it necessary for the Parseval's theory to hold? Please let me know if my theory is wrong.

Couple of notes about the code:


apologies if it isnt efficient coding, I'm still a novice
dividing the fft by (fs*(win*win).sum()) is something I found through research and also a part of the inbuilt Welch code
lastly I am using trapezoidal rule for integration
my FFT energy is a couple of magnitudes smaller (ie E-09)


With 0Hz peak removed

With 0Hz peak



import numpy as np
import matplotlib.pyplot as plt
import os
import math
from scipy.fftpack import fft, ifft, fftfreq, rfft
from scipy.signal import detrend, get_window

plt.clf()
plt.close()

###############################################################################

config = '1perc_damping'

dir0 = os.getcwd()
dir1 = dir0+'/'+config
dir2 = dir1+'/Summary files'

data = np.loadtxt(dir1 + '/000deg_500rpm_An_Mxx.dat')

###############################################################################

fs = 625
nperseg = 1024
noverlap = 512
window_meth = 'hanning'


###############################################################################

#detrend
data = data - np.mean(data)


#get windowing method
win = get_window(window_meth, nperseg)

# number of segments
nseg = np.ceil((len(data)-nperseg)/(nperseg-noverlap)) + 1 

total_n = nperseg + (nseg-1)*(nperseg-noverlap)

# padding with zeros
if total_n &gt; len(data):
    n = total_n - len(data)
    padding = np.zeros(np.int(n))
    data = np.concatenate((data,padding))


# fft 

fft_data = 0

for i in range(0,np.int(nseg)):

    start = i*(nperseg-noverlap)

    f_data = np.multiply(data[start:start+nperseg],win) 
    f_data = f_data- np.mean(f_data)   #### Detrend line ### 

    transf = fft(f_data)
    transf_c = (transf*np.conj(transf))/(fs*(win*win).sum())

    fft_data = fft_data + transf_c


Pxx = np.real(fft_data)/(nseg)
fx = fftfreq(nperseg, 1/fs)

plt.figure(1)
plt.plot(fx[0:512],2*Pxx[0:512])
#plt.plot(fx,Pxx)

plt.figure(2)
x = np.linspace(0,len(data),len(data))
plt.plot(x,data)

#Finding RMS
l = np.int(len(fx)/2 -1)
c = 2
A = 0

for i in range(1,l+1):
    A = A + 0.5*(c*Pxx[i]+c*Pxx[i-1])*(fx[i]-fx[i-1])

RMS = np.sqrt(A)

# Checking Parseval's theorem
data_energy = 0
FFT_energy = 0

for i in range(1,l+1):
    FFT_energy = FFT_energy + 0.5*(np.power(c*Pxx[i],2)+np.power(c*Pxx[i-1],2))*(fx[i]-fx[i-1])

for i in range(1,len(data)):
    data_energy = data_energy + 0.5*(np.power(abs(data[][1][i]),2)+np.power(abs(data[i-1]),2))*(i-(i-1))*(1/625)
RMS value


  Firstly I am getting a 0Hz peak even after removing the mean at the start. The only way I can remove the 0Hz peak is to take the mean away after I have segmented and multiplied the window.


The mean doesn't have to remain at 0 when you apply a window. Changing the weight of a particular sample would shift the mean accordingly (e.g. a weight greater than 1 on a positive sample would increase the mean). If you must remove the 0Hz peak, the only guaranteed way to do so is right before the FFT (or right after the FFT, but that second option wouldn't remove the spectral leakage to adjacent bins).

Parseval's relationship


  My second issue is that I cant get the Parseval's theorem to hold. Is there something I'm doing wrong? Is it necessary for the Parseval's theory to hold? Please let me know if my theory is wrong


Since you are dealing with discrete time signals, you should be using the specific formulation of Parseval's theorem for discrete time signals. This formulation is derived directly from the definition of the Discrete Fourier Transform. It then has to hold between the input and output data to the fft call. Note however that different implementations may use different normalization which can affect the exact formulation by a constant scaling factor. With the normalization used in numpy.fft.fft, you should have the following relationship:

$$
  \sum_{n=0}^{N-1} |x_n|^2 = \frac{1}{N} \sum_{k=0}^{N-1} |X_k|^2
$$

Note however that your code currently incorrectly throws away the imaginary parts of your computed spectrum on the following line:

Pxx = np.real(fft_data)/(nseg)


You final computation should instead make use of squared-magnitudes. Since you are later taking the square of each frequency bin value, you may take the absolute value:

Pxx = np.abs(fft_data)/(nseg)


Your will also need to remove the frequency step part in your final integration of FFT_energy and stick with the straight summation given in the discrete time formulation of Parseval's theorem given above.

for i in range(0,len(Pxx)):
    FFT_energy = FFT_energy + np.power(Pxx[i],2)
FFT_energy = FFT_energy / len(Pxx)

for i in range(0,len(data)):
    data_energy = data_energy + np.power(abs(data[i]),2)