冲击可以通过峭度指标来检测，是如何来检测的，python示例代码

问题

冲击可以通过峭度指标来检测，是如何来检测的，可以给1个示例代码吗

思路

• 带冲击的信号其峭度指标>3
• 不带冲击的信号其峭度指标在3左右
• 可以通过滑动窗来检测在哪一段

示例代码

带冲击的信号峭度指标值

import numpy as np
import matplotlib.pyplot as plt

def generate_signal_with_impact(impact_time=0.5, impact_magnitude=10, duration=1, sampling_rate=1000):
    """
    生成含有单个冲击的信号
    :param impact_time: 冲击发生的时间，以秒为单位
    :param impact_magnitude: 冲击的幅度
    :param duration: 信号的总持续时间，以秒为单位
    :param sampling_rate: 采样率
    :return: 时间序列和信号值
    """
    t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
    signal = np.random.normal(0, 1, size=t.shape)  # 生成随机噪声信号
    impact_index = int(impact_time * sampling_rate)
    signal[impact_index] += impact_magnitude  # 在指定时间添加冲击
    return t, signal

def calculate_kurtosis(signal):
    """
    计算信号的峭度
    :param signal: 输入信号
    :return: 峭度值
    """
    x = signal
    X_rms2 = np.sum(x**2)/len(x)  # 3.均方值
    X_rms = np.sqrt(X_rms2)   # 4.均方根值(有效值)
    X_beta = np.mean( np.power(x, 4) )   # 12.峭度
    X_kf = X_beta/X_rms ** 4     # 18.峭度指标
    return X_kf

# 生成含有冲击的信号
t, signal = generate_signal_with_impact()

# 计算信号的峭度
kurtosis_value = calculate_kurtosis(signal)
print(f"Kurtosis of the signal: {kurtosis_value}")

# 绘制信号
plt.figure(figsize=(10, 4))
plt.plot(t, signal, label='Signal with Impact')
plt.title(f"Signal and its Kurtosis: {kurtosis_value:.2f}")
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.show()

Kurtosis of the signal: 12.830772158435776

不带冲击的信号峭度指标值

import numpy as np
import matplotlib.pyplot as plt

def generate_signal_with_impact(impact_time=0.5, impact_magnitude=10, duration=1, sampling_rate=1000):
    """
    生成含有单个冲击的信号
    :param impact_time: 冲击发生的时间，以秒为单位
    :param impact_magnitude: 冲击的幅度
    :param duration: 信号的总持续时间，以秒为单位
    :param sampling_rate: 采样率
    :return: 时间序列和信号值
    """
    t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
    signal = np.random.normal(0, 1, size=t.shape)  # 生成随机噪声信号
    impact_index = int(impact_time * sampling_rate)
#     signal[impact_index] += impact_magnitude  # 在指定时间添加冲击
    return t, signal

def calculate_kurtosis(signal):
    """
    计算信号的峭度
    :param signal: 输入信号
    :return: 峭度值
    """
    n = len(signal)
    mean = np.mean(signal)
    std_dev = np.std(signal)
    sum_diff = np.sum((signal - mean) ** 4)
    kurtosis = (n * sum_diff) / ((n - 1) * (n - 2) * (std_dev ** 4))
    return kurtosis

# 生成含有冲击的信号
t, signal = generate_signal_with_impact()

# 计算信号的峭度
kurtosis_value = calculate_kurtosis(signal)
print(f"Kurtosis of the signal: {kurtosis_value}")

# 绘制信号
plt.figure(figsize=(10, 4))
plt.plot(t, signal, label='Signal with Impact')
plt.title(f"Signal and its Kurtosis: {kurtosis_value:.2f}")
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.show()

Kurtosis of the signal: 3.0646772222342067

基于滑动窗的冲击信号检测

import numpy as np
import matplotlib.pyplot as plt

def generate_signal_with_impacts(impacts=[(0.2, 5), (0.5, 7)], sampling_rate=1000, duration=1):
    """
    生成含有多个冲击的信号
    :param impacts: 冲击列表，每个元素为(冲击位置, 冲击幅度)
    :param sampling_rate: 采样率
    :param duration: 信号持续时间
    :return: 时间数组和信号数组
    """
    t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
    signal = np.random.normal(0, 1, size=t.shape)  # 生成随机噪声信号
    for impact_position, impact_magnitude in impacts:
        impact_index = int(len(t) * impact_position)
        signal[impact_index] += impact_magnitude  # 添加冲击
    return t, signal

def calculate_sliding_kurtosis(signal, window_size=50):
    """
    计算信号的滑动窗口峭度
    :param signal: 输入信号
    :param window_size: 滑动窗口大小
    :return: 峭度值数组
    """
    kurtosis_values = []
    for i in range(len(signal) - window_size + 1):
        window = signal[i:i + window_size]
        n = len(window)
        mean = np.mean(window)
        std_dev = np.std(window)
        sum_diff = np.sum((window - mean) ** 4)
        kurtosis = (n * sum_diff) / ((n - 1) * (n - 2) * (std_dev ** 4))
        kurtosis_values.append(kurtosis)
    return np.array(kurtosis_values)

# 生成信号
t, signal = generate_signal_with_impacts()

# 计算滑动窗口峭度
kurtosis_values = calculate_sliding_kurtosis(signal)

# 绘制信号和峭度
plt.figure(figsize=(14, 6))

plt.subplot(2, 1, 1)
plt.plot(t, signal, label='Signal')
plt.title('Signal with Impacts')
plt.xlabel('Time')
plt.ylabel('Amplitude')

plt.subplot(2, 1, 2)
plt.plot(t[:len(kurtosis_values)], kurtosis_values, label='Sliding Kurtosis', color='orange')
plt.axhline(y=np.mean(kurtosis_values) + 2*np.std(kurtosis_values), color='r', linestyle='--', label='Threshold')
plt.title('Sliding Window Kurtosis')
plt.xlabel('Time')
plt.ylabel('Kurtosis')

plt.tight_layout()
plt.legend()
plt.show()