title: Using numpy/scipy to identify slope changes in digital signals?
body: <p>I am trying to come up with a generalised way in Python to identify pitch rotations occurring during a set of planned spacecraft manoeuvres. You could think of it as a particular case of a <a href="https://en.wikipedia.org/wiki/Step_detection" rel="noreferrer">shift detection</a> problem.</p>
  
  <p>Let's consider the <code>solar_elevation_angle</code> variable in my set of measurements, identifying the elevation angle of the sun measured from the spacecraft's instrument. For those who might want to play with the data, I saved the <code>solar_elevation_angle.txt</code> file <a href="https://pastebin.com/xyY11isP" rel="noreferrer">here</a>.</p>
  
  <pre><code>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()
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/ffyBT.png" rel="noreferrer"><img src="https://i.stack.imgur.com/ffyBT.png" alt="Solar elevation angle plot"></a></p>
  
  <p>The scanline is my time dimension. The four points where the slope changes identify the spacecraft pitch rotations. </p>
  
  <p>As you can see, the solar elevation angle evolution outside the spacecraft manoeuvres regions is pretty much linear as a function of time, and that should always be the case for this particular spacecraft (except for major failures). </p>
  
  <p>Note that during each spacecraft manoeuvre, the slope change is obviously continuous, although discretised in my set of angle values. That means: for each manoeuvre, it does not really make sense to try to locate a single scanline where a manoeuvre has taken place. My goal is rather to identify, for each manoeuvre, a "representative" scanline in the range of scanlines defining the interval of time where the manoeuvre occurred (e.g. middle value, or left boundary).  </p>
  
  <p>Once I get a set of "representative" scanline indexes where all manoeuvres have taken place, I could then use those indexes for rough estimations of manoeuvres durations, or to automatically place labels on the plot.  </p>
  
  <p>My solution so far has been to:</p>
  
  <ol>
  <li>Compute the 2nd derivative of the solar elevation angle using
  <code>np.gradient</code>.   </li>
  <li>Compute absolute value and clipping of resulting
  curve. The clipping is necessary because of what I assume to be
  discretisation noise in the linear segments, which would then severely affect the identification of the "real" local maxima in point 4. </li>
  <li>Apply smoothing to the resulting curve, to get rid of multiple peaks. I'm using scipy's 1d gaussian filter with a trial-and-error sigma value for that.</li>
  <li>Identify local maxima.</li>
  </ol>
  
  <p>Here's my code:</p>
  
  <pre><code>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()
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/b0Roj.png" rel="noreferrer"><img src="https://i.stack.imgur.com/b0Roj.png" alt="enter image description here"></a></p>
  
  <p>I can then easily identify the local maxima by using scipy's <code>argrelmax</code> function:</p>
  
  <pre><code>max_idx = argrelmax(smoothed_signal)[0]
  print(max_idx)
  # [ 689 1019 2356 2685]
  </code></pre>
  
  <p>Which correctly identifies the scanline indexes I was looking for:</p>
  
  <pre><code>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()
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/ZBMQn.png" rel="noreferrer"><img src="https://i.stack.imgur.com/ZBMQn.png" alt="enter image description here"></a></p>
  
  <p>My question is: Is there a better way to approach this problem?<br>
  I find that having to manually specify the clipping threshold values to get rid of the noise and the sigma in the gaussian filter weakens this approach considerably, preventing it to be applied to other similar cases.</p>