import numpy as np
import matplotlib.pyplot as plt
import pywt
from scipy import signal as sig
from scipy.signal import find_peaks
from scipy.stats import skew, kurtosis
import pandas as pd
from prettytable import PrettyTable
plt.rcParams['figure.figsize'] = (12,6)
def process(df):
df.drop(['TIME AXIS - Plot 0',
'TIME AXIS - Plot 1',
'TIME AXIS - Plot 2',
'TIME AXIS - Plot 3',
'TIME AXIS - Plot 4',
'TIME AXIS - Plot 5'],axis=1,inplace=True)
df.dropna(inplace=True)
def autoscales(N, dt, dj, wf='morlet', p=6):
"""Compute scales as fractional power of two."""
if wf == 'dog':
s0 = (dt * np.sqrt(p + 0.5)) / np.pi
elif wf == 'paul':
s0 = (dt * ((2 * p) + 1)) / (2 * np.pi)
elif wf == 'morlet':
s0 = (dt * (p + np.sqrt(2 + p ** 2))) / (2 * np.pi)
else:
raise ValueError('wavelet function not available')
J = int(np.floor(dj ** -1 * np.log2((N * dt) / s0)))
s = np.empty(J + 1)
for i in range(s.shape[0]):
s[i] = s0 * 2 ** (i * dj)
return s
def smooth_and_find_peaks(signal):
smoothed_signal = sig.savgol_filter(signal, 51, 2)
peaks, _ = find_peaks(smoothed_signal)
return smoothed_signal, peaks
def plot_signal_with_peaks(ax, t, ns_signal, sm_signal, peaks, title):
ax.plot(t, ns_signal, 'navy', label='noisy')
ax.plot(t, sm_signal, 'firebrick', label='smooth')
ax.plot(t[peaks], sm_signal[peaks], 'ro')
ax.set_title(title)
ax.set_xlabel('Time')
ax.set_ylabel('Amplitude')
ax.legend()
def plot_spectrogram(ax, signal, title):
# Spectrogram
freq, time, sxx = sig.spectrogram(signal)
plot = ax.contourf(time, freq, np.log10(sxx), cmap='jet', aspect='auto')
ax.set_title(f'Spectrogram of {title}')
ax.set_ylabel('Frequency')
ax.set_xlabel('Time')
plt.colorbar(plot, ax=ax)
def plot_cwt(ax,signal,title):
# Continuous Wavelet Transform
wavelet = 'morl' # Complex Morlet wavelet
widths = autoscales(len(signal), 1, 0.6)
cwmatr, freq_cwt = pywt.cwt(signal, widths, wavelet)
cwtmatr = np.abs(cwmatr)**2
plot = ax.imshow(cwtmatr, cmap='jet', aspect='auto', extent=[0, 1, min(widths), max(widths)])
ax.set_title(f'CWT of {title}')
ax.set_ylabel('Scales')
ax.set_xlabel('Time')
plt.colorbar(plot, ax=ax)
def calculate_power_density_spectrum(signal, fs=1000):
freq, psd = sig.welch(signal, fs, nperseg=1024)
return freq, psd
def calculate_temporal_features(signal, fs):
mean_amp = np.mean(signal)
std_dev = np.std(signal)
skewness = skew(signal)
kurt = kurtosis(signal)
return mean_amp, std_dev, skewness, kurt
def calculate_spectral_features(freq, power_spectral_density):
max_freq = freq[np.argmax(power_spectral_density)]
total_power = np.sum(power_spectral_density)
normalized_power = power_spectral_density / total_power
spectral_centroid = np.sum(normalized_power * freq)
spectral_spread = np.sqrt(np.sum(normalized_power * (freq - spectral_centroid)**2))
return max_freq, spectral_centroid, spectral_spread
# def calculate_vector_strength(peaks, t):
# # Calculate angular phases from peak times
# phases = 2 * np.pi * t[peaks]
# # Calculate mean resultant vector
# mean_resultant_vector = np.mean(np.exp(1j * phases))
# # Calculate vector strength
# vector_strength = np.abs(mean_resultant_vector)
# return vector_strength
def analyze_and_plot_signals(ns_signals,sig_names, t, fs=1000):
num_signals = len(ns_signals)
# Create subplots
fig, axs = plt.subplots(4, num_signals, figsize=(20, 15))
axs = axs.flatten()
# Create tables for numerical values
columns_temporal = ['Signal', 'Mean Amplitude', 'Std Deviation', 'Skewness', 'Kurtosis']
columns_spectral = ['Signal', 'Max Frequency', 'Spectral Centroid', 'Spectral Spread']
# columns_vector_strength = ['Signal', 'Vector Strength']
table_temporal = PrettyTable(columns_temporal)
table_spectral = PrettyTable(columns_spectral)
# table_vector_strength = PrettyTable(columns_vector_strength)
for i, ns_signal in enumerate(ns_signals, start=1):
sm_signal, peaks = smooth_and_find_peaks(ns_signal)
# Plot the signals with detected peaks
plot_signal_with_peaks(axs[i - 1], t, ns_signal, sm_signal, peaks, f'Smoothed {sig_names[i-1]}')
# Plot CWT and Spectrogram
plot_spectrogram(axs[i + num_signals - 1], sm_signal, f' {sig_names[i-1]}')
plot_cwt(axs[i + 2 * num_signals - 1], sm_signal, f' {sig_names[i-1]}')
# Calculate and plot Power Density Spectrum
freq_psd, psd = calculate_power_density_spectrum(sm_signal, fs)
plt.ylim(0.5e-9,0)
axs[i + 3 * num_signals - 1].semilogy(freq_psd, psd)
axs[i + 3 * num_signals - 1].set_title(f'Power Spectral Density ({sig_names[i-1]})')
axs[i + 3 * num_signals - 1].set_xlabel('Frequency (Hz)')
axs[i + 3 * num_signals - 1].set_ylabel('Power/Frequency (dB/Hz)')
# Calculate and print Temporal Features
temporal_features = calculate_temporal_features(sm_signal, fs)
temporal_features = [round(value, 6) for value in temporal_features]
# Calculate and print Spectral Features
spectral_features = calculate_spectral_features(freq_psd, psd)
spectral_features = [round(value, 6) for value in spectral_features]
# Calculate and print Vector Strength
# vector_strength = calculate_vector_strength(peaks, t)
# vector_strength = round(vector_strength, 6)
# Add data to PrettyTable
# Add data to PrettyTables
table_temporal.add_row([f' {sig_names[i-1]}'] + temporal_features)
table_spectral.add_row([f'{sig_names[i-1]}'] + spectral_features)
# table_vector_strength.add_row([f' {sig_names[i]}', vector_strength])
plt.tight_layout()
plt.show()
# Print PrettyTables
print("Temporal Features:")
print(table_temporal)
print("\nSpectral Features:")
print(table_spectral)High Aerated Slug
HAS_5_600 = pd.read_excel('run1.xlsx')HAS_5_600
process(HAS_5_600)
HAS_5_600_sen4 = HAS_5_600['PRESSURE AXIS RANGE 0- 0.1 BAR - Plot 4'][:60000]
time = HAS_5_600_sen4.index[:60000]
plt.plot(time,HAS_5_600_sen4);
HAS_5_600_sen4_sm = sig.savgol_filter(HAS_5_600_sen4,51,2)
plt.plot(time,HAS_5_600_sen4_sm)
plt.plot(time,HAS_5_600_sen4,'navy')
plt.plot(time,HAS_5_600_sen4_sm,'firebrick')
plt.xlabel('Time (ms)')
plt.ylabel('Pressure (bar)')
plt.title('HAS_6_500 Smoothed Data')
plt.show()
def autoscales(N, dt, dj, wf='morlet', p=6):
"""Compute scales as fractional power of two.
:Parameters:
N : integer
number of data samples
dt : float
time step
dj : float
scale resolution (smaller values of dj give finer resolution)
wf : string
wavelet function ('morlet', 'paul', 'dog')
p : float
omega0 ('morlet') or order ('paul', 'dog')
:Returns:
scales : 1d numpy array
scales
"""
if wf == 'dog':
s0 = (dt * sqrt(p + 0.5)) / np.pi
elif wf == 'paul':
s0 = (dt * ((2 * p) + 1)) / (2 * np.pi)
elif wf == 'morlet':
s0 = (dt * (p + np.sqrt(2 + p ** 2))) / (2 * np.pi)
else:
raise ValueError('wavelet function not available')
# See (9) and (10) at page 67.
J = int(np.floor(dj ** -1 * np.log2((N * dt) / s0)))
s = np.empty(J + 1)
for i in range(s.shape[0]):
s[i] = s0 * 2 ** (i * dj)
return s
widths = autoscales(len(HAS_5_600_sen4_sm),1,0.6)
cwtmatr1, freqs1 = pywt.cwt(HAS_5_600_sen4_sm,widths,'morl')
plt.imshow(np.abs(cwtmatr1)**2, cmap='jet',aspect='auto')
plt.colorbar()
plt.title('Scalogram (HAS_5.5_500)')
plt.xlabel('Time')
plt.ylabel('Scales')
plt.show()
np.arange(cwtmatr1.shape[0])
np.arange(cwtmatr1.shape[1])
# Plotting Scalogram in 3D
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')
X, Y = np.meshgrid(np.arange(cwtmatr1.shape[1]),np.arange(cwtmatr1.shape[0]))
sur=ax.plot_surface(X, Y, np.abs(cwtmatr1)**2, cmap='jet')
ax.set_xlabel('Time')
ax.set_ylabel('frequency')
ax.set_zlabel('Amplitude')
ax.set_title('Scalogram in 3D')
plt.colorbar(sur,shrink=0.5)
plt.show()
Low Aerated Slug
LAS_3_600 = pd.read_excel('run1.xlsx')
LAS_3_600
process(LAS_3_600)
time = LAS_3_600.index[:60000]
LAS_3_600_sen4 = LAS_3_600['PRESSURE AXIS RANGE 0- 0.1 BAR - Plot 4'][:60000]
plt.plot(time, LAS_3_600_sen4);
LAS_3_600_sen4_sm = sig.savgol_filter(LAS_3_600_sen4,51,2)
plt.plot(time,LAS_3_600_sen4_sm);
plt.plot(time,LAS_3_600_sen4,'navy')
plt.plot(time,LAS_3_600_sen4_sm,'firebrick')
plt.xlabel('Time (ms)')
plt.ylabel('Pressure (bar)')
plt.title('LAS_3_600 Smoothed Data')
plt.show()
widths = autoscales(len(LAS_3_600_sen4_sm),1,0.6)
print(widths)
[1.93602662e+00 2.93446762e+00 4.44782118e+00 6.74163625e+00
1.02184098e+01 1.54882129e+01 2.34757409e+01 3.55825695e+01
5.39330900e+01 8.17472780e+01 1.23905704e+02 1.87805928e+02
2.84660556e+02 4.31464720e+02 6.53978224e+02 9.91245629e+02
1.50244742e+03 2.27728445e+03 3.45171776e+03 5.23182579e+03
7.92996503e+03 1.20195794e+04 1.82182756e+04 2.76137421e+04
4.18546064e+04]
cwtmatr2, freqs2 = pywt.cwt(LAS_3_600_sen4_sm,widths,'morl')
plt.imshow(np.abs(cwtmatr2)**2, cmap='jet',aspect='auto')
plt.colorbar()
plt.title('Scalogram (LAS_3_600)')
plt.xlabel('Time')
plt.ylabel('Scales')
plt.show()
# Plotting Scalogram in 3D
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')
X, Y = np.meshgrid(np.arange(cwtmatr2.shape[1]),np.arange(cwtmatr2.shape[0]))
sur= ax.plot_surface(X, Y, np.abs(cwtmatr2)**2, cmap='jet')
ax.set_xlabel('Time')
ax.set_ylabel('frequency')
ax.set_zlabel('Amplitude')
ax.set_title('Scalogram in 3D')
plt.colorbar(sur,shrink=0.5)
plt.show()
担任《Mechanical System and Signal Processing》审稿专家，担任《中国电机工程学报》，《控制与决策》等EI期刊审稿专家，擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。
知乎学术咨询：https://www.zhihu.com/consult/people/792359672131756032?isMe=1重庆成都太热了，贵阳安顺避暑几天。
1082
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
