Для моего проекта мне пришлось вручную кодировать метод Уэлча, используя код ниже. В значительной степени это включает в себя нахождение спектральной плотности через FFT и включает в себя окна и сегментирование с перекрытием.
Мне нужно знать, правильный ли мой PSD, поэтому я проверяю двумя способами.
- корень области под PSD должен быть среднеквадратическим (E [x ^ 2]), и поскольку я убрал свое среднее значение из данных, это должно равняться стандартному отклонению данных
- Проверка выполнения теоремы Парсеваля.
Есть две основные проблемы, с которыми у меня проблемы. Во-первых, я получаю пик 0 Гц даже после удаления среднего значения в начале. Единственный способ удалить пик 0 Гц - убрать среднее значение после того, как я сегментировал и умножил окно. (#### линия тренда ####)
При пике точка 1 хорошо держится, но без пика стандартное отклонение немного отличается от среднеквадратичного значения, но это не представляет большой проблемы.
Моя вторая проблема заключается в том, что я не могу получить теорему Парсеваля. Я что-то не так делаю? Нужно ли придерживаться теории Парсеваля? Пожалуйста, дайте мне знать, если моя теория неверна.
Пара замечаний по коду:
- извиняюсь, если это не эффективное кодирование, я все еще новичок
- деление fft на (fs * (win * win) .sum ()) - это то, что я нашел в результате исследований, а также часть встроенного кода Уэлча
- наконец, я использую трапециевидное правило для интеграции
- моя энергия БПФ на пару величин меньше (т.е. E-09)
PSD с удаленным пиком 0 Гц
PSD с пиком 0 Гц
</p>
<pre>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)