4

I have the the noisy curve defined by numpy 2D array: mEPSC

As you can see, it has the first flat segment, then rise, peak and decay phases. I need to find the starting point of the rise phase, marked here by the red dot. How do I do that in python?

Axon
  • 547
  • 2
  • 6
  • 14
  • Can you take an average of the y over a certain distance, get the difference from the previous value and the use the first negative value? – Eyrofire Apr 15 '14 at 15:12
  • @Eyrofire, I have conidered that, but this method is very noise-sensitive. The curve I have shown here is relatively smooth, but I may need to handle the data with a couple of orders of magnutude smaller signal to noise ratio. – Axon Apr 15 '14 at 15:21
  • 1
    It seems you need a good algorithm first - the best way to smooth/filter the data and still preserve the *inflection* point, You may want to ask over in http://dsp.stackexchange.com/. Once you have that, you can return here with your Python implementation if you need to. They are ```NumPy``` and ```SciPy``` *aware* over there. – wwii Apr 15 '14 at 15:31
  • Are you in need of finding the actual inflection point (i.e., where the concavity changes, up or down)? Or are you trying to identify the starting point the spectrum where the signal is strong? The red point is not an inflection point, strictly speaking. – Taro Sato Apr 15 '14 at 17:51
  • @Taro Sato, I need the starting point. Sorry, I am not thery good at math terminology even in my native language. ☺ – Axon Apr 15 '14 at 18:24

2 Answers2

8

If the data look like the one in the example figure, you could estimate the background and its noise level and apply some threshold to extract the portion of data that are above the background. The example follows:

#!/usr/bin/env python2.7
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import gaussian_filter


def generate_fake_data():
    """Generate data that looks like an example given."""
    xs = np.arange(0, 25, 0.05)
    ys = - 20 * 1./(1 + np.exp(-(xs - 5.)/0.3))
    m = xs > 7.
    ys[m] = -20.*np.exp(-(xs - 7.)[m] / 5.)

    # add noise
    ys += np.random.normal(0, 0.2, xs.size)
    return xs, ys


def main():
    xs, ys = generate_fake_data()

    # smooth out noise
    smoothed = gaussian_filter(ys, 3.)

    # find the point where the signal goes above the background noise
    # level (assumed to be zero here).
    base = 0.
    std = (ys[xs < 3] - base).std()
    m = smoothed < (base - 3. * std)
    x0 = xs[m][0]
    y0 = ys[m][0]

    plt.plot(xs, ys, '.')
    plt.plot(xs, smoothed, '-')
    plt.plot(x0, y0, 'o')
    plt.show()


if __name__ == '__main__':
    main()

enter image description here

Taro Sato
  • 1,444
  • 1
  • 15
  • 19
  • Well, it is not exact, but close, thank you. Only, gaussian filter affects rise kinetics too much, I'd rather not use any filter at all. – Axon Apr 15 '14 at 19:08
  • Smoothing is not really necessary here, but I just used it to estimate where the signal starts contributing more to the data. If you deal with noisier data than used in this example, however, you probably need some way of modeling that transition; smoothing is just a quickie. – Taro Sato Apr 15 '14 at 19:12
8

Well, I calculated the local differentials along the curve for the small dt and the extremum of the derivative curve pointed out the "inflection point" quite well. I think, I'll settle with that. enter image description here

Axon
  • 547
  • 2
  • 6
  • 14
  • And what is the code that you used? This is nice but not useful to others. Please update your answer. – My Work Apr 19 '23 at 13:16