Здравствуйте, мне нужно найти точку перегиба (IP) зашумленного сигнала, чтобы определить, где функция изменяет вогнутость. Сигнал следующий: 
Я намерен вычислить как первую, так и вторую производные (d1, d2), чтобы получить d1 = min значение и d2 = 0 при том же значении X. В этом случае, если мы вычисляем d1 для каждого X от 51K до 68K, будет значение X, при котором мы получим минимальное значение d1. Это минимальное значение d1 объясняет, где находится точка перегиба (IP). Но минимального значения d1 недостаточно, чтобы определить, где расположен IP, фактически при том же значении X (d1 = min значение) d2 должно быть d2 = 0 .
ПРИМЕЧАНИЕ: Если d2 <0 </em> у нас есть вогнутость вниз, если d2> 0 , у нас есть извилистость вверх и d2 = 0 имеем точку перегиба (IP).
Для моих данных я ожидал, что значения d1 и d2 будут различаться, как на изображении ниже, но я получаю очень несогласованные значения :
Мой код:
df = pd.read_csv("data.csv", sep=',', decimal=",")
df
# seconds
x = np.array(df.iloc[:, 0].astype('float'))
# pressure
y = np.array(df.iloc[:, 1].astype('float'))
print(x)
array([51641., 51743., 51845., 51966., 52081., 52196., 52317., 52439.,
52553., 52662., 52771., 52873., 52988., 53097., 53205., 53320.,
53537., 53652., 53767., 53882., 53991., 54106., 54221., 54329.,
54425., 54540., 54655., 54776., 54891., 55006., 55121., 55223.,
55338., 55453., 55568., 55689., 55804., 55926., 56041., 56137.,
56245., 56354., 56475., 56584., 56692., 56807., 56909., 57030.,
57145., 57260., 57394., 57509., 57630., 57732., 57841., 57962.,
58077., 58198., 58306., 58421., 58549., 58651., 58766., 58887.,
59008., 59123., 59238., 59353., 59455., 59564., 59686., 59807.,
59928., 60043., 60158., 60260., 60362., 60477., 60592., 60707.,
60841., 60956., 61084., 61199., 61301., 61422., 61544., 61665.,
61773., 61888., 62003., 62112., 62227., 62335., 62444., 62559.,
62667., 62782., 62897., 63011., 63114., 63229., 63344., 63458.,
63586., 63701., 63797., 63906., 64014., 64123., 64237., 64359.,
64467., 64582., 64697., 64914., 65029., 65144., 65252., 65361.,
65469., 65565., 65674., 65782., 65897., 66019., 66134., 66248.,
66363., 66478., 66587., 66702., 66810., 66938., 67059., 67161.,
67257., 67379., 67494., 67609., 67730., 67852., 67966., 68081.,
68196., 68299., 68420., 68535., 68663., 68784.])
print(y)
array([6.02, 6.02, 6.02, 6.01, 6. , 6. , 6. , 5.98, 5.98, 5.98, 5.97,
5.96, 5.95, 5.95, 5.94, 5.93, 5.92, 5.91, 5.9 , 5.89, 5.89, 5.88,
5.87, 5.86, 5.86, 5.85, 5.84, 5.83, 5.81, 5.8 , 5.79, 5.78, 5.77,
5.76, 5.75, 5.74, 5.72, 5.71, 5.7 , 5.69, 5.67, 5.66, 5.64, 5.63,
5.62, 5.6 , 5.58, 5.58, 5.55, 5.53, 5.51, 5.49, 5.48, 5.46, 5.44,
5.42, 5.4 , 5.38, 5.36, 5.34, 5.32, 5.3 , 5.27, 5.26, 5.23, 5.2 ,
5.18, 5.16, 5.14, 5.12, 5.09, 5.07, 5.04, 5.02, 5. , 4.97, 4.95,
4.93, 4.91, 4.88, 4.85, 4.83, 4.81, 4.79, 4.77, 4.74, 4.72, 4.7 ,
4.68, 4.66, 4.64, 4.63, 4.61, 4.6 , 4.58, 4.56, 4.54, 4.53, 4.51,
4.49, 4.48, 4.47, 4.45, 4.43, 4.42, 4.41, 4.4 , 4.38, 4.37, 4.35,
4.35, 4.33, 4.32, 4.31, 4.3 , 4.27, 4.27, 4.25, 4.24, 4.23, 4.22,
4.22, 4.2 , 4.19, 4.19, 4.18, 4.16, 4.16, 4.14, 4.13, 4.12, 4.11,
4.11, 4.1 , 4.09, 4.08, 4.07, 4.06, 4.05, 4.05, 4.03, 4.03, 4.02,
4.01, 4. , 4. , 3.99, 3.99, 3.97, 3.97])
Прочтите данные и попробуйте для их сглаживания
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline
# smoothing the signal
y_spl = UnivariateSpline(x, y, s=0, k=3)
x_range = np.linspace(x[0], x[-1], len(y))
# set figure
plt.figure(figsize=(16, 6))
plt.grid(True)
# plot noisy data
plt.semilogy(x,y, 'bo', label = 'data', lw=1)
# plot smoothed data
plt.semilogy(x_range, y_spl(x_range), 'r', lw=2.5)
# aesthetics
plt.legend(['Noisy Data','Smoothed Data'])
plt.title('Pressure on Seconds')
plt.xlabel('Seconds')
plt.ylabel('Pressure')
plt.show()
Вычислить d1 и d2
# compute the d1 and d2 for every single point of x
y_spl_1d = y_spl.derivative(n=1) #d1
y_spl_2d = y_spl.derivative(n=2) #d2
# d1 array
d1_arr = y_spl_1d(x_range)
d1_arr = np.round(d1_arr, 6)
# d2 array
d2_arr = y_spl_2d(x_range)
d2_arr = np.round(d2_arr, 6)
# enrich the df
df['d1'] = d1_arr
df['d2'] = d2_arr
Найти координату X, Y для d1 = минимальное значение
Примечание: Где d1 = минимальное значение, d2 ≠ 0
# getting the row with the min value of d1
min_d1 = df.loc[df['d1'].idxmin()]
print("Minimum d1 value: ", '\n', min_d1)
Minimum d1 value:
Seconds 59008
Pressure 5.23076923076923
Temperature 100.06105006105
Minutes 983
d1 -0.000265
d2 -1e-06
Name: 514, dtype: object
# determine the x and y intercept
xmin_d1 = float(min_d1[0])
ymin_d1 = float(min_d1[1])
print("X Intercept: ", xmin_d1, '\n')
print("Y Intercept: ", ymin_d1)
Постройте график тренда d1
# plot d1
plt.figure(figsize=(14, 6))
plt.grid()
plt.plot(x_range, d1_arr, 'o')
plt.axvline(x = xmin_d1, color='orange')
plt.title('First Derivative')
plt.xlabel('Seconds')
plt.ylabel('d1')
plt.show()
Как вы можете видеть, значения d1 уменьшаются с течением времени, но у них много колебаний, это должно быть вызвано тем, что данные недостаточно сглажены.
Перехват точки перегиба (IP) на исходных данных с использованием минимального значения d1
# plot d1 inteception on function
plt.figure(figsize=(16, 6))
plt.grid()
plt.semilogy(x,y, 'bo', label = 'data')
plt.semilogy(x_range, y_spl(x_range), 'r', lw=2.5)
plt.title('Pressure on Seconds')
plt.xlabel('Seconds')
plt.ylabel('Pressure')
plt.axvline(x = xmin_d1, color = 'orange', lw=1, ls='--')
plt.axhline(y = ymin_d1, color = 'orange', lw=1, ls='--')
plt.text(x = xmin_d1+2000, y = ymin_d1, s = 'IP', fontsize=15)
plt.show()
Примечание: Это хороший шаг, но его недостаточно, мне нужно d2 = 0 при том же значении X
Построение тренда d2
# plot d2
plt.figure(figsize=(14, 6))
plt.grid()
plt.plot(x_range, d2_arr)
plt.axvline(x = xmin_d1, color = 'orange', lw=1, ls='--')
plt.title('Second Derivative')
plt.xlabel('Seconds')
plt.ylabel('d2')
plt.show()
Как вы можете см. тренд d2 очень колеблется, он так часто перескакивает от <0 к> 0, что невозможно определить ни вогнутость, ни точку перегиба.
Вопросы?
- Слежу ли я за правильным процессом?
- Что не так в моем коде?
- Связана ли проблема с процессом сглаживания?
- Можете ли вы найти лучшее решение?
Большое спасибо за вашу помощь.