Сглаживание кривой со многими пиками с помощью Gaussian

ProgramBox

Сглаживание кривой со многими пиками с помощью Gaussian

Post author:admin
Запись опубликована:19 февраля, 2023
Post category:Вопросы по программированию

#python #numpy #scipy #convolution

#python #numpy #scipy #свертка

Вопрос:

У меня есть данные спектроскопии с некоторыми очень резкими пиками, как видно на синей кривой. Я хотел бы сделать пики немного более плавными, как оранжевая кривая на графике.

Я подумал, что самый простой способ сделать это — свернуть мои точки данных с гауссианами. Я знаю оба numpy и scipy имею convolve функции, но я не уверен, нужна ли мне 1D или 2D свертка, чтобы получить то, что мне нужно. До сих пор я пробовал convolve1d и gaussian_filter1d от scipy , и convolve от numpy . Ни один из них не улучшил четкие линии, соединяющие точки данных. Я также не знаю, как выбрать правильную сигму или веса…

Текстовый файл, содержащий точки данных, находится здесь .

Оранжевая кривая генерируется из программы визуализации, и я хочу иметь возможность генерировать ее самостоятельно, python а не с помощью программы.

Редактировать:

Новая ссылка для файла

1. являются ли «пандусы» (т. Е. От 1500-3100) Частью данных или просто артефактами из plt.plot точек соединения?

2. Также какой код используется для создания этого изображения?

3. @DanielF нет, эта область просто соединяет две точки данных.

4. глядя на ваши данные, у меня есть серьезные сомнения в том, что желтый график отражает то, что происходит на самом деле. Что бы ни делала программа, это очень опасно для количественного анализа ваших данных.

5. @mikuszefski Я увеличил интенсивность оранжевого графика, чтобы он был более четким (чтобы оранжевый и синий не перекрывались и имели четкие линии). Это данные спектроскопии, поэтому нас волнует, где находятся пики. Программа — это то, что обычно используется для такого типа анализа данных многими людьми, поэтому я думаю, что все должно быть в порядке… Как вы думаете, почему это опасно?

Ответ №1:

Похоже, вам нужен «оценщик плотности ядра», который реализуется:

 from scipy.stats import gaussian_kde

X = np.random.rand(50) * 3500
Y = np.random.rand(50) * 50
xi = linspace(0, 3500, 1000)

kde = gaussian_kde(X, weights = Y, bw_method = .01)  #tune `bw_method` to get the bandwidth you want
plt.plot(xi, kde.pdf(xi))

Вам также может потребоваться настроить y масштабирование графика в соответствии с вашими требованиями

1. Итак, в вашем примере X это x координаты синей кривой, xi x координаты оранжевой кривой и kde.pdf(xi) координаты y оранжевой кривой. Однако у меня есть только координаты X для синей кривой… Поэтому я не понимаю, как ваш пример может мне помочь, если вы также не покажете мне, для чего я должен использовать xi .

2. Я не могу получить доступ к веб-сайту, на котором находится ваш текстовый файл, с моего рабочего компьютера, поэтому вам, вероятно, следует включить имеющиеся у вас данные (хотя бы пример) в сам вопрос.

Ответ №2:

Это вручную воспроизводит оранжевую кривую, указанную в OP. Оказывается, она свернута с лоренцевым, а не гауссовым.

 
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import find_peaks
from scipy.optimize import curve_fit

def gs( x, x0, a, s ):
    return a * np.exp( -( x - x0 )**2 / ( 2 * s**2 ) )

def cs( x, x0, a, s ):
    return a / ( ( x - x0 )**2   s**2 )

conrange = 40000 

### gasiian is no good
# ~condata = np.fromiter( ( gs(x, 0, 1, 1800 ) for x in np.arange( -5000, 5000 ) ), np.float )
### Cauchy looks much better
condata = np.fromiter( 
    ( 
        cs( x, 0, 1, 2000 ) for x in np.arange( -conrange, conrange ) 
    ), np.float
)
### shift can be zero. 
### Amplitude does not matter as it will be scaled later anyway
### width matters of course, but is adjusted manually for the moment.

data = np.loadtxt("ir_data.txt")
xdata = data[:, 0]
ydata = data[:, 1]

xdataint = np.fromiter( ( int( x* 100 ) for x in xdata ), int ) 
xmin = xdataint[0]
xmax = xdataint[-1]
xfilled = np.arange( xmin , xdataint[-1]   1 )
yfilled = np.zeros( len( xfilled ), dtype=np.float )
xfloat = np.fromiter( ( x / 100. for x in xfilled), float ) 


for x, y in zip( xdataint, ydata ):
    yfilled[ x - xmin ] = y
### just putting a manual scale here, but the real one can be calculated
### from the convolution properties
yc = 1e6 * np.convolve( condata, yfilled, mode="full" )

xfull = np.arange(
    -conrange   xmin, xmin   conrange   len( xfilled ) - 1
)
xfloat = np.fromiter( ( 0.01 * x for x in xfull ), float )

fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.plot( xdata, ydata, ls='', marker='o', ms=2 )
ax.plot( xfloat, yc, ls='-')
plt.show()

Отказ от ответственности

Это предварительные результаты, опубликованные только по запросу автора OP. Может быть некоторое уточнение.