Применение теоремы Парсевала с методом Уэлча

Question

Для моего проекта мне пришлось вручную кодировать метод Уэлча, используя код ниже. В значительной степени это включает в себя нахождение спектральной плотности через FFT и включает в себя окна и сегментирование с перекрытием.

Мне нужно знать, правильный ли мой PSD, поэтому я проверяю двумя способами.

1. корень области под PSD должен быть среднеквадратическим (E[x^2]), и, поскольку я убрал свое среднее значение из данных, это должно равняться стандартному отклонению данных
2. Проверка выполнения теоремы Парсеваля.

У меня две основные проблемы. Во-первых, я получаю пик 0 Гц даже после удаления среднего значения в начале. Единственный способ удалить пик 0 Гц - убрать среднее значение после того, как я сегментировал и умножил окно. (#### линия тренда ####)

С пиком точка 1 хорошо держится, но без пика стандартное отклонение немного отличается от среднеквадратичного значения, но это не представляет большой проблемы.

Вторая проблема заключается в том, что я не могу получить теорему Парсеваля. Я что-то не так делаю? Нужно ли придерживаться теории Парсеваля? Пожалуйста, дайте мне знать, если моя теория неверна.

Пара заметок о коде:

• извиняюсь, если это не эффективное кодирование, я все еще новичок
• деление fft на (fs*(win*win).sum()) - это то, что я нашел в ходе исследований, а также часть встроенного кода Уэлча
• наконец, я использую правило трапеции для интеграции
• моя энергия БПФ на пару величин меньше (т.е. E-09)

PSD с удаленным пиком 0 Гц

PSD с пиком 0 Гц

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)