Свернуть экпоненту гауссом

[SVD-d]

День добрый!

По ходу исследований полупроводников столкнулся с такой задачей. Есть, значится, релаксация, описывается она (если правильно помню термин) затухающей экспонентой. То есть,
$f(t) = A \exp \left( -\frac{t}{\tau} \right)$
Если мы говорим о физике дела, то эта релаксация описывает эмиссию носителей с какого-то центра. Такие релаксации довольно легко рассчитываются, анализируются, для них есть свой вполне работающий математический аппарат (вроде анализа производной по логарифму времени) и так далее.
Но если бы простая формула всегда работала, жизнь была бы неинтересной. Если у нас в одном образце есть два уровня, с которых идёт эмиссия, то и экспонент становится две, а формула усложняется:
$f(t) = A_{1} \exp \left( -\frac{t}{\tau_{1}} \right) + A_{2} \exp \left( -\frac{t}{\tau_{2}} \right)$
Если мы смотрим на релаксацию, то в общем виде её уже не так-то просто проанализировать. Впрочем, задачу можно упростить, если рассматривать не саму релаксацию, которая монотонная и непонятно когда выйдет на насыщение, а производную по логарифму времени: $df/dln(t) = - \frac{At}{\tau} \exp \left( -\frac{t}{\tau} \right)$, которая имеет явно выраженный пик в точке $t = \tau$ с ожидаемой формой и известной полушириной.

Это было вступление. При работе с образцами мы обнаружили, что уровень, который у нас есть, не простой, а размытый, то есть, на самом деле состоящий из большого количества находящихся рядышком уровней, составляющий целую зону. Самое простое здесь - предположить, что уровень размыт по гауссу. То есть, мы имеем некоторую функцию амплитуды
$A(\tau) = \frac{A_{0}}{2\pi \sigma} \cdot \exp \left( -\frac{(\ln^2(\tau)-\ln^2(\tau_{0}))^2}{2\sigma^2} \right)$
То что размытие по логарифму тау - это так и надо, оно, по-видимому, так и происходит.

И лично мне кажется, что итоговая релаксация будет иметь вид свёртки функций f и A. То есть,
$ f(t) = \frac{A_{0}}{2\pi \sigma} \cdot \int_{0}^{\infty}{\exp \left( -\frac{(\ln^2(\tau)-\ln^2(\tau_{0}))^2}{2\sigma^2} \right) \cdot \exp \left( -\frac{t}{\tau} \right) d \tau} $

И чтобы этого не показалось мало, практически всегда у нас не один размытый уровень, а два или три.

Конечо, понятно, что производная свёртки будет свёрткой этого самого пика. Но даже зная это, когда мы имеем дело с реальными измерениями (пятнадцать тысяч точек, шум порядка пяти процентов, а иногда и выше), численными методами подбирать три пика, размытых по гауссу, с неизвестными амплитудами, временами и сигмами - удовольствие то ещё. Задача, как правило, решается неоднозначно, и хорошо ещё если решается вообще. Самая главная проблема - это, конечно, численное решение. Если скормить Левенбергу набор точек, он довольно быстро застывает далеко от правильного ответа. Оно и понятно: когда в том же ориджине мы пишем программу для фиттинга, мы не можем брать настоящий интеграл и ограничиваемся простым суммированием сотни экспонент в пределах трёх сигма. И эта сотня экспонент считается долго и ненадёжно.

Я в своё время подумал насчёт фурье-преобразования. Это бы позволило уйти от интегрирования, что увеличило бы точность. Но как должен выглядеть фурье-образ такой свёртки - рассчитать не могу.
Может, кто-нибудь здесь знает, что в таких случаях полагается делать? Или уже занимался подобным? А то я покопался, но ничего не нашёл.

[SVD-d]

Вообще, если построить кривую $f(t) = A \exp \left( -\frac{t}{\tau} \right)$ с какими-нибудь заданными параметрами, и получить из неё FFT, то можно получить вполне ожидаемую кривую, которую даже можно продифференцировать по логарифму времени и получить пик в точке $\lambda = \frac{1}{2 \pi \tau}$. В той же точке есть пик у мнимой части FFT, даже без производной. Но я никак не могу это использовать, потому что не могу откопать аналитический вид функции фурье-образа затухающей экспоненты.
Раньше я думал, что раз $f(t) = A \exp \left( -\frac{t}{\tau} \right)$, то $\hat{f}(\lambda) = \frac{1/\tau}{1/\tau^2 + t^2}$. Но это почему-то не так. ЧЯДНТ?

[Munin]

SVD-d в сообщении #1169678 писал(а):

Но я никак не могу это использовать, потому что не могу откопать аналитический вид функции фурье-образа затухающей экспоненты.

Поищите фурье-образ экспоненты, помноженной на функцию Хевисайда.

[SVD-d]

Хммм... Вы, наверное, правы.
То есть, не то чтобы мне понравился внешний вид функции, которая получается при перемножении. Но с другой стороны, затухающая экспонента, умноженная на Хевисайда визуально равна функции $Ce^{-a|x|}$, умноженной на него же. Соответственно, когда я взял свой набор точек и симметрично достроил его в минус, прибавив правильный ноль, FFT для неё стал легко подгоняться функцией $\hat{f}(\lambda) = Ca/(a^2+\labmda^2)$.

Правда, теперь мне стало совсем непонятно, как свёртка такой функции гауссом должна выглядеть в фурье-пространстве. То есть, понятно как должна, но непонятно, почему она выглядит не так, как я рассчитывал, когда мне было понятно. Но это ладно, надо будет всё аккуратно пересчитать. Тогда, может быть, всё сойдётся.

[SVD-d]

В итоге, как я понимаю, фурье-образ свёртки затухающей экспоненты гауссом
$f(t) = A \int^{\infty}_{0} {\exp \left( -\frac{(\ln(\tau)-\ln(\tau_{0}))^2}{2\sigma^2} \right) \cdot \exp \left( -\frac{|t|}{\tau} \right) d\ln(\tau)}$
выглядит примерно вот так:
$\hat{f}(\lambda) = A_{0} \left( 1 - \frac{1}{\sqrt{2\pi}} \exp \left( -\frac{(\ln(1/\lambda)-\ln(\tau_{0}))^2}{\sigma^2/2} \right) \right) \cdot \left( \frac{1/\tau_{0}}{1/\tau_{0}^{2} + \lambda^{2}} \right)$
Кажется, я что-то забыл, но с точностью до константы сойдёт.
Всем спасибо, все расходимся.

[SVD-d]

Называется, рано радовался.

Я наивно надеялся, что раз фурье преобразует гаусса в гаусса, то FFT свёртки затухающей экспоненты функцией гаусса будет равно поточечному перемножению фурье-образа экспоненты и гаусса.
Но я как-то упустил из виду тот факт, что гаусс у меня не простой, а по логарифму, то есть,
$g(\tau) = \exp \left( -\frac{(\ln(\tau)-\ln(\tau_{0}))^2}{2\sigma^2} \right)$
Я искренне верил в то, что FFT такой функции будет в логарифмическом масштабе тоже гауссом. На всякий случай даже попробовал смоделировать пару кривых, сделал им Фурье, разделил на FFT классической экспоненты и убедился в том, что результат в логарифмическом масштабе похож на колокол. Даже проверил, чтобы при $\sigma \rightarrow 0$ кривая схлопывалась в одну классическую экспоненту.
А теперь я проверяю, сходится ли моделированная кривая с ожидаемой при достаточно больших $\sigma$ (в районе 1.5-2.5, например) - и нифига она не сходится.

Как думаете, проблема в том, что гаусс - не гаусс? Или гаусс на самом деле гаусс, а значит я что-то неправильно моделирую?

[10110111]

$\exp(-\ln(x)^2)=\exp(\ln(x))^{-\ln(x)}=x^{-\ln(x)}$ — не очень похоже на гауссиану.

[mserg]

Данные можете выложить?

[SVD-d]

10110111 в сообщении #1171718 писал(а):

$\exp(-\ln(x)^2)=\exp(\ln(x))^{-\ln(x)}=x^{-\ln(x)}$ — не очень похоже на гауссиану.

В логарифмическом масштабе вполне похоже : )
Между прочим, я довольно сильно удивился, когда увидел, что FFT набора точек, представляющих гаусса в логарифмическом масштабе, в логарифмическом масштабе выглядел как гаусс. Но потом подумал и решил, что вроде всё логично. А потом ещё подумал, и решил, что что-то здесь не так. А потом подумал ещё, и решил спросить.

С данными несколько сложнее. Данные - это пятнадцать тысяч точек, и в них ещё довольно много шума. Чтобы проверить, работает ли система вообще, я беру модель: точки распределены равномерно от 2e-4 до 3 секунд, а выходной сигнал - $f(t) = A \exp \left( -\frac{t}{\tau} \right)$, где $\tau = 0.3, A = 1$, например.
А потом, чтобы получить интеграл, беру сумму из полутора сотен таких релаксаций, где амплитуду A рассчитываю по гауссу, меняя $\ln(\tau)$ в пределах трёх сигма.