python列表转化为数字信号_python – 使用numpy/scipy识别数字信号的斜率变化？

我试图在

Python中提出一种通用的方法来识别在一组计划的航天器机动过程中发生的俯仰旋转.您可以将其视为

shift detection问题的特例.

让我们考虑一组测量中的solar_elevation_angle变量,确定从航天器仪器测量的太阳仰角.对于那些可能想要使用数据的人,我保存了solar_elevation_angle.txt文件here.

import numpy as np

import matplotlib.pyplot as plt

from matplotlib import gridspec

from scipy.signal import argrelmax

from scipy.ndimage.filters import gaussian_filter1d

solar_elevation_angle = np.loadtxt("solar_elevation_angle.txt", dtype=np.float32)

fig, ax = plt.subplots()

ax.set_title('Solar elevation angle')

ax.set_xlabel('Scanline')

ax.set_ylabel('Solar elevation angle [deg]')

ax.plot(solar_elevation_angle)

plt.show()

扫描线是我的时间维度.斜率变化的四个点确定航天器俯仰旋转.

正如您所看到的,航天器机动区域外的太阳高度角演变与时间的关系几乎是线性的,对于这个特定的航天器而言应该始终如此(主要故障除外).

请注意,在每次航天器操作过程中,坡度变化显然是连续的,尽管在我的角度值组中是离散的.这意味着：对于每次机动,尝试定位单个扫描线并进行机动动作并没有多大意义.我的目标是为每个操纵识别扫描线范围内的“代表性”扫描线,该扫描线限定机动发生的时间间隔(例如,中间值或左边界).

一旦我获得了一组“代表性”扫描线索引,其中所有操作都已发生,我可以使用这些索引粗略估计操纵持续时间,或自动在标绘上放置标签.

到目前为止我的解决方案是：

>使用计算太阳高度角的二阶导数

np.gradient.

>计算绝对值并剪切结果

曲线.由于我的假设,裁剪是必要的

线性段中的离散噪声,这将严重影响点4中“实际”局部最大值的识别.

>对生成的曲线应用平滑,以消除多个峰.我正在使用scipy的1d高斯滤波器,其具有试错法sigma值.

>识别局部最大值.

这是我的代码：

fig = plt.figure(figsize=(8,12))

gs = gridspec.GridSpec(5, 1)

ax0 = plt.subplot(gs[0])

ax0.set_title('Solar elevation angle')

ax0.plot(solar_elevation_angle)

solar_elevation_angle_1stdev = np.gradient(solar_elevation_angle)

ax1 = plt.subplot(gs[1])

ax1.set_title('1st derivative')

ax1.plot(solar_elevation_angle_1stdev)

solar_elevation_angle_2nddev = np.gradient(solar_elevation_angle_1stdev)

ax2 = plt.subplot(gs[2])

ax2.set_title('2nd derivative')

ax2.plot(solar_elevation_angle_2nddev)

solar_elevation_angle_2nddev_clipped = np.clip(np.abs(np.gradient(solar_elevation_angle_2nddev)), 0.0001, 2)

ax3 = plt.subplot(gs[3])

ax3.set_title('absolute value + clipping')

ax3.plot(solar_elevation_angle_2nddev_clipped)

smoothed_signal = gaussian_filter1d(solar_elevation_angle_2nddev_clipped, 20)

ax4 = plt.subplot(gs[4])

ax4.set_title('Smoothing applied')

ax4.plot(smoothed_signal)

plt.tight_layout()

plt.show()

然后我可以使用scipy的argrelmax函数轻松识别局部最大值：

max_idx = argrelmax(smoothed_signal)[0]

print(max_idx)

# [ 689 1019 2356 2685]

这正确地标识了我正在寻找的扫描线索引：

fig, ax = plt.subplots()

ax.set_title('Solar elevation angle')

ax.set_xlabel('Scanline')

ax.set_ylabel('Solar elevation angle [deg]')

ax.plot(solar_elevation_angle)

ax.scatter(max_idx, solar_elevation_angle[max_idx], marker='x', color='red')

plt.show()

我的问题是：有没有更好的方法来解决这个问题？

我发现必须手动指定削波阈值以消除噪声和高斯滤波器中的西格玛,这大大削弱了这种方法,使其无法应用于其他类似的情况.