Определение длительности высокочастотного импульса
Немного из текущей работы.
Строю вот такие картинки
Задача состоит в том, чтобы найти параметры, при которых генерируется самый короткий импульс. Цветом на картинке выше передана величина, равная обратной длительности, посчитанная разными способами. Проблема в том, что вместе с импульсом генерируется ещё очень много всего. Характерный сигнал имеет приблизительно такой вид :
Короткий импульс - это всплеск в центре.
Таким образом, надо придумать автоматизированный способ как-то оценить эту длительность. У меня родилось три идеи, как это сделать.
Первый: давайте просто найдём расстояние по x между максимумом и минимумом. В python-библиотеке numpy это делается просто:
def pulse_duration_1(d):
return numpy.argmin(d) - numpy.argmax(d)
Очевидная проблема возникает в таких случаях:
Второй способ заключается в том, чтобы определить максимальную (вернее, минимальную, потому что в интересущей области она отрицательна) производную заданной функции. В numpy это делает следующей функцией:
def pulse_duration_2(d):
grad = numpy.gradient(d)
return -1./grad.min()
Этот способ плох тем, что из знания производной извлечь длительность не так-то просто, да и наличие однозначной связи между ними не очевидно.
Тогда я решил работать с отфильтрованным сигналом. Если убрать низкочастотную часть, сигнал принимает следующий вид:
Надо найти длительность импульса. Для этого я, во-первых, нахожу огибающую преобразованием Гильберта, а затем ширину получившейся функции. Сначала я попробовал сделать это по ширине на полувысоте. То есть вычитаю из анализируемой функции половину её максимума, а затем нахожу нули. В numpy/scipy соответствующая функция выглядит следующим образом:
def pulse_duration_3(d):
f = numpy.fft.rfft(d)
for i in range(1000):
f[i] = 0.
filtered = numpy.fft.irfft(f)
envelope = abs(scipy.signal.hilbert(d))
zero_crossings = numpy.where(numpy.diff(numpy.sign(envelope - envelope.max()/2.)))[0]
return zero_crossings[1] - zero_crossings[0]
Проблема в том, что такой подход может завышать длительность в некоторых сложных случаях.
Альтернативно, можно определять длительность импульса по огибащей, делая фиттинг, например, гауссом. Код в numpy/scipy здесь сложнее (подсмотрено здесь):
def pulse_duration_4(d):
f = numpy.fft.rfft(d)
for i in range(1000):
f[i] = 0.
filtered = numpy.fft.irfft(f)
envelope = abs(scipy.signal.hilbert(d))
def gaussian(x, sig):
return envelope.max()*numpy.exp(-(x-envelope.argmax())**2/(sig**2))
def fit(p, x):
return numpy.sum([gaussian(x, p) for i in xrange(len(p))], axis=0)
err = lambda p, x, y: fit(p,x)-y
x = xrange(len(envelope))
results, value = scipy.optimize.leastsq(err, [2], args=(x, envelope))
return results[0]
По моим наблюдениям, этот подход завышает длительность почти всегда. Можно сравнить, например, огибающую и фиттинг на иллюстрации ниже:
Собственно, на первой картинке и приведено сравнение этих четырёх подходов (порядок соответствует изложению в посте). И видно, что разные подходы дают разные результаты (даже если не обращать внимания на абсолютные значения). Мне, однако, кажется, что самый адекватный - третий. А что посоветовали бы вы?
Строю вот такие картинки
Задача состоит в том, чтобы найти параметры, при которых генерируется самый короткий импульс. Цветом на картинке выше передана величина, равная обратной длительности, посчитанная разными способами. Проблема в том, что вместе с импульсом генерируется ещё очень много всего. Характерный сигнал имеет приблизительно такой вид :
Короткий импульс - это всплеск в центре.
Таким образом, надо придумать автоматизированный способ как-то оценить эту длительность. У меня родилось три идеи, как это сделать.
Первый: давайте просто найдём расстояние по x между максимумом и минимумом. В python-библиотеке numpy это делается просто:
def pulse_duration_1(d):
return numpy.argmin(d) - numpy.argmax(d)
Очевидная проблема возникает в таких случаях:
Второй способ заключается в том, чтобы определить максимальную (вернее, минимальную, потому что в интересущей области она отрицательна) производную заданной функции. В numpy это делает следующей функцией:
def pulse_duration_2(d):
grad = numpy.gradient(d)
return -1./grad.min()
Этот способ плох тем, что из знания производной извлечь длительность не так-то просто, да и наличие однозначной связи между ними не очевидно.
Тогда я решил работать с отфильтрованным сигналом. Если убрать низкочастотную часть, сигнал принимает следующий вид:
Надо найти длительность импульса. Для этого я, во-первых, нахожу огибающую преобразованием Гильберта, а затем ширину получившейся функции. Сначала я попробовал сделать это по ширине на полувысоте. То есть вычитаю из анализируемой функции половину её максимума, а затем нахожу нули. В numpy/scipy соответствующая функция выглядит следующим образом:
def pulse_duration_3(d):
f = numpy.fft.rfft(d)
for i in range(1000):
f[i] = 0.
filtered = numpy.fft.irfft(f)
envelope = abs(scipy.signal.hilbert(d))
zero_crossings = numpy.where(numpy.diff(numpy.sign(envelope - envelope.max()/2.)))[0]
return zero_crossings[1] - zero_crossings[0]
Проблема в том, что такой подход может завышать длительность в некоторых сложных случаях.
Альтернативно, можно определять длительность импульса по огибащей, делая фиттинг, например, гауссом. Код в numpy/scipy здесь сложнее (подсмотрено здесь):
def pulse_duration_4(d):
f = numpy.fft.rfft(d)
for i in range(1000):
f[i] = 0.
filtered = numpy.fft.irfft(f)
envelope = abs(scipy.signal.hilbert(d))
def gaussian(x, sig):
return envelope.max()*numpy.exp(-(x-envelope.argmax())**2/(sig**2))
def fit(p, x):
return numpy.sum([gaussian(x, p) for i in xrange(len(p))], axis=0)
err = lambda p, x, y: fit(p,x)-y
x = xrange(len(envelope))
results, value = scipy.optimize.leastsq(err, [2], args=(x, envelope))
return results[0]
По моим наблюдениям, этот подход завышает длительность почти всегда. Можно сравнить, например, огибающую и фиттинг на иллюстрации ниже:
Собственно, на первой картинке и приведено сравнение этих четырёх подходов (порядок соответствует изложению в посте). И видно, что разные подходы дают разные результаты (даже если не обращать внимания на абсолютные значения). Мне, однако, кажется, что самый адекватный - третий. А что посоветовали бы вы?