1
$\begingroup$

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.

  1. 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.
  2. 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 > 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)

$\endgroup$
1
$\begingroup$

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)
$\endgroup$
  • $\begingroup$ Hi, thanks for this, it was really helpful. However I still cant get the energies to match. Note I did multiply the FFT_data by the conjugate which gives me a 0j component. I did try your method @SleuthEye with a couple of different changes (ie with and without the scaling factor) and still cant get the energies to equal $\endgroup$ – P. Dis Aug 29 '18 at 23:40

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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