Аппроксимация экспериментально полученных кривых

[Sanyok]

Привет!

Задача формулируется следующим образом:
Есть экспериментально полученная кривая F0(t)
Найти (или рассчитать) гладкую кривую F1(t) такую, что бы
sum((F1-F0)^2) -> min, то есть F1 должна быть близка к F0 в смысле наименьших квадратов.

Про кривую F0 известно следующее:
F0' < 0 для всех t > 0,
F0'' > 0 для всех t > 0.

[Dan_Te]

Цитата:

sum((F1-F0)^2) -> min

что такое сумма функции? Может, интеграл?

Приближаем на конечном отрезке или на всей положительной полупрямой? Известно что-нибудь про функцию, которой аппроксимируем (например, что это экспонента, или дробно-рациональная функция)?

[Sanyok]

Dan_Te писал(а):

Цитата:

sum((F1-F0)^2) -> min

что такое сумма функции? Может, интеграл?

Для непрерывной функции - это, конечно, интеграл. Просто экспериментально полученные кривульки - это всегда набор дискретных значений, вот я и написал sum (просто индексы не хотелось писать).

Цитата:

Приближаем на конечном отрезке или на всей положительной полупрямой?

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

Цитата:

Известно что-нибудь про функцию, которой аппроксимируем (например, что это экспонента, или дробно-рациональная функция)?

Кроме того, что я написал - очень немного. Эти кривые хорошо аппроксиммируются функциями y(t) = A*t^(f(t)), где f(t) - слабо меняется на начальном участке отрезка (где-то на 60 % его длины).

[незванный гость (1)]

Увы, сказать можно немного.

1) Очевидно, через любое конечное множество точек проходит соответствующий интерполяционный полином. Они [полиномы сии] крайне нестабильны (для высоких порядков) и сейчас почти не применяются (насколько я знаю).

2) Если речь идет о практических подсчетах, часто применяют сплайны. В двух ипостасях - интерполирующие и аппроксимирующие. Первые дают точное решение в заказанных точках, вторые более гладкие на всей области определеня. Мне - нравятся, использовал часто и много.

3) Если речь идет о природе функции, то тут начинается высокое искуство, за которое платят немалые деньги (страховые, финансовые и маркетинговые компании в США). Нужно попытаться угадать, какие факторы могут быть. Факторы образуют "базис" решения, и затем идет подбор параметров.

Например, мы можем искать решение в виде:

f(x) = a + bx + c x^2 + d exp(g x)

Далее, применяем метод наименьших квадратов, т.е. минимизируем сумму S(y -f(x)) по всем x. Теперича эта S(a,b,c,d,g) суть функция гладкая, и минимум ее можно найти дифферецируя по параметрам и решая уравнения (численно).

~~

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

~~~~~~~~~~~~~~~~~~~~~~~

Простите великодушно, коли вопрос не так понял, да банальные истины понаписал... Если нужны более подробные ссылки по какому пункту - пишите, подберу, что могу.

[Dan_Te]

Абстрактно построить гладкую кривую - так задача обычно не ставится.
Всегда можно интерполировать полиномами Лагранжа (если у нас имеется N точек, то степень многочлена будет не более N-1), который будет попросту совпадать с теоретической функцией во всех точках х_1...х_n, в которых эта функция измерена. Естественно, при этом ошибка (сумма квадратов разностей значений теоретический функции и многочлена в точках х_1...х_n) будет равна нулю, но вряд ли такой результат можно назвать приемлемым, потому что между точками х_1...х_n многочлен может очень сильно колбасить.

Обычно все-таки делают предположение о виде теоретической функции (например, как у вас, y(t) = A*t^(f(t)), только надо еще как-нибудь расписать f(t) ), получается зависимость от нескольких параметров. Теперь надо найти такие значения параметров, чтобы ошибка была минимальной. В качестве ошибки обычно берут сумму квадратов разностей, но можно и сумму модулей брать, и еще что-нибудь.

Вам это требуется? Тогда вам в помощь программа Statistica, пакет Nonlinear estimation. Задаете вид теоретической функции, выбираете функцию ошибок, и программа вам рассчитывает значения параметров.

[незванный гость (1)]

Цитата:

например, как у вас, y(t) = A*t^(f(t)

Позвольте добавить, что в этом случае часто логарифмируют, и затем ищут аппроксимацию вида:

Log(y(t)) ~~ Log(A) + Log(t) * f(t)

В этом виде применять метод наименьших квадратов проще. (Ну а Log(A), понятное дело, при поиске на b заменяют.)

[Sanyok]

В-общем, пора карты раскрывать, дабы все стало яснее...

Речь идет о банальной очистке сигнала от помех. Сигнал убывает во времени , и с увеличением времени -> 0.
Полиномы и сплайны использовать нет смысла - частота выборки намного больше частоты помехи, и поскольку результирующая кривая будет проходить через все точки исходной, то она будет не лучше исходной.

Пока что из всех способов, что я испробовал, наилучшим оказалась фильтрация КИХ-фильтром с окном переменной ширины (в начале, где сигнал быстро меняется - окно узкое, и расширяется по мере уменьшения скорости изменения сигнала).

Мой вопрос состоял в том - можно ли, руководствуясь только требованиями к производным функции(условия - F0' < 0 и F0'' > 0) и условием "наименьших квадратов", разработать сходящийся алгоритм (т.е., на каждом шаге которого ошибка (сумма квадратов разностей) будет уменьшаться), позволяющий численно рассчитать функцию, удовлетворяющую этим условиям. При этом не делая никаких предположений о типе функции.

[Sanyok]

незванный гость (1) писал(а):

:evil:

Цитата:

например, как у вас, y(t) = A*t^(f(t)

Позвольте добавить, что в этом случае часто логарифмируют, и затем ищут аппроксимацию вида:

Log(y(t)) ~~ Log(A) + Log(t) * f(t)

В этом виде применять метод наименьших квадратов проще. (Ну а Log(A), понятное дело, при поиске на b заменяют.)

С логарифмированием там туго... Дело в том, что значения искомого сигнала (они все положительные, тут проблем нет) меньше в несколько раз, чем амплитуда помехи (а помеха - знакопеременная), короче - почти половина значений исходного (грязного) сигнала - отрицательная.

[Sanyok]

Цитата:

2) Если речь идет о практических подсчетах, часто применяют сплайны. В двух ипостасях - интерполирующие и аппроксимирующие. Первые дают точное решение в заказанных точках, вторые более гладкие на всей области определеня. Мне - нравятся, использовал часто и много.

Не понял - разве аппроксиммирующие сплайны не проходят через "заказанные" точки?

[незванный гость (1)]

Sanyok писал(а):

Не понял - разве аппроксиммирующие сплайны не проходят через "заказанные" точки?

Сглаживающие сплайны ( - по-видимому, я изобрел термин "аппроксимирующие":oops:) не обязаны проходить через исходные точки, в отличие от интерполяционных. Если я правильно помню, то интерполяционный слайн минимизирует Integral(f''(x),{x}] при f(x[k]) == y[k]. Сглаживающий-же минимизирует Integral(f''(x),{x}] + M Summa[(y[k]-f(x[k]))^2, {k}]. От выбора M зависит, естественно плотность прилегания. При больших M страдает гладкость, при малых - точность аппроксимации. Понятное дело, количество сегментов может быть заметно меньше, чем количество данных.

Если речь идет об обработке экспериментальных данных, то, может, Вам стоит посмотреть книжки по цифровой обработке сигналов. Разного рода фильтры имеют разные характеристики, и дают часто результаты зело удивительные. Тут уж простор необычайный - и Butterworth, и Чебышев, и Calman - всего не перечислишь. Извлекают сигнал из-под помехи поразительно.

[незванный гость (1)]

Во-первых, простите, что не всю историю прочитал. Посему, кое-что в моем предыдущем post'е не в тему - отвечено до вопроса, заданного Вами.

1)

незванный гость писал(а):

Сглаживающий-же минимизирует Integral(f''(x),{x}] + M Summa[(y[k]-f(x[k]))^2, {k}].

Позвольте уточнить - Integral(f''(x),{x}] + Summa[m[k] (y[k]-f(x[k]))^2, {k}].[/u]. Т.е. в общем случае вес зависит от координаты. Таким образом, в разных зонах можно требовать разное приближение.

2) Вы не могли бы уточнить - требуется ли обработка в реальном времени или нет. От этого кое-что может зависеть. В частности, стоит ли применять КИХ-фильтры, или те, что я упоминал.

3) Если это возможно, попробуйте обрабатывать функцию с конца, т.е. x -> XX - x. В Вашем случае, это даст фильтру устаканиться около нуля. При этом параметры фильтра (например, размер окна, или временную постоянную) можно менять наблюдая рост и скорость роста результата. Еще один интересный параметр наблюдения, который стоит учитывать - это уровень шума. Его тоже проще ловить, когда сигнал близок к нулю...

4) Если известна природа помехи, то можно применять специализированные фильтры. Например, резонансный подавитель 50 и 100 герцовых частот. У меня был весьма успешный опыт с ним...

5) Осмелюсь посоветовать - будьте осторожнее с изменением размеров окна. Изменение размера есть скачкообразное изменение функции, и может привести к появлению артефактов - разрывов, потери гладкости. Лучше фильтры с плавным изменением параметов...

[Гость]

незванный гость (1) писал(а):

:evil:
Во-первых, простите, что не всю историю прочитал. Посему, кое-что в моем предыдущем post'е не в тему - отвечено до вопроса, заданного Вами.

1)

незванный гость писал(а):

Сглаживающий-же минимизирует Integral(f''(x),{x}] + M Summa[(y[k]-f(x[k]))^2, {k}].

Позвольте уточнить - Integral(f''(x),{x}] + Summa[m[k] (y[k]-f(x[k]))^2, {k}].[/u]. Т.е. в общем случае вес зависит от координаты. Таким образом, в разных зонах можно требовать разное приближение.

Я вообще не знал, что есть сплайны, не проходящие через точки исходной кривой... К своему стыду... Может, это как раз то, что мне нужно... Где можно об этом почитать?

незванный гость писал(а):

2) Вы не могли бы уточнить - требуется ли обработка в реальном времени или нет. От этого кое-что может зависеть. В частности, стоит ли применять КИХ-фильтры, или те, что я упоминал.

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

незванный гость писал(а):

4) Если известна природа помехи, то можно применять специализированные фильтры. Например, резонансный подавитель 50 и 100 герцовых частот. У меня был весьма успешный опыт с ним...

Помеха - широкополосная 50-герцовка (куча гармоник, и четные, и нечетные). Иногда бывает, что на одной кривой характер помехи меняется (спектр гармоник, их фазы) - тогда вообще караул - даже широкое окно плохо справляется (но это бывает нечасто). Интересно, как подавить все ЧЕТНЫЕ гармоники... С нечетными даже простым усреднением по периоду можно справиться....
Еще один отрицательный момент - частота этой промышленной помехи редко бывает стабильной (в пределах +/- 1 Гц уж точно гуляет).

незванный гость писал(а):

5) Осмелюсь посоветовать - будьте осторожнее с изменением размеров окна. Изменение размера есть скачкообразное изменение функции, и может привести к появлению артефактов - разрывов, потери гладкости. Лучше фильтры с плавным изменением параметов...

Ширина окна у меня счас - это гладкая функция от времени. Подбирал ее методом "проб и ошибок" довольно долго.. Я в обычной фильтрации потихоньку разочаровываюсь - частенько спектр помех и спектр полезного сигнала перекрываются (на кривой это самое проблемное место - там, где измеряемый сигнал имеет наименьший радиус кривизны - КИХ-фильтр это место искажает больше всего). [/b]

[Sanyok]

Sorry - предыдущий пост - мой.. Зайти забыл....

[незванный гость (1)]

Sanyok писал(а):

Где можно об этом почитать?

В книгах - увы, порекомендовать могу немного.

[1] Carl De Boor, A Practical Guide to Splines, Springer-Verlag, 1978.
Сразу видно - старая. И программы на Фортране .

[2] Математика и САПР, Мир, 1988 (том 1)
Здесь и вовсе - больше упоминание, чем описание. Может книгу найти проще, однако.

Можно поискать на интернете - "smoothing spline" на Google выдает хренову тучу страниц. Может, что и полезное найдется. И еще - где-то я недавно видел ссылку на страницу NIST'а, посвященную алгоритмам. Чуть ли не на этом сайте. NIST обычно очень добротен.

В целом - народ, похоже, ушел в фильтры и DSP. "Там - чужие слова, там - дурная молва, там ненужные встречи случаются..."

Sanyok писал(а):

Помеха - широкополосная 50-герцовка

Ну с ней-то и впрямь резонансный подавитель работает. И хорошо работает! Попробуйте, очень советую... Особенно с конца - он дает артефакт в момент инициализации.

Кстати, еще вопрос - а много ли у Вас данных. То есть, надо ли Вам делать полностью автоматическую обработку, или Вас устроит полуавтоматический вариант - ну, например, прогнали раз, в точке излома узлов подкинули, прогнали двас .

[Sanyok]

Цитата:

Ну с ней-то и впрямь резонансный подавитель работает. И хорошо работает! Попробуйте, очень советую... Особенно с конца - он дает артефакт в момент инициализации.

Я бы рад, но такой термин ("резонансный подавитель") мне встречается впервые... Может, поясните слегка, что это за алгоритм (или ссылочку, если можно)? Это не обычный режекторный фильтр, насколько я понимаю?

Цитата:

Кстати, еще вопрос - а много ли у Вас данных. То есть, надо ли Вам делать полностью автоматическую обработку, или Вас устроит полуавтоматический вариант - ну, например, прогнали раз, в точке излома узлов подкинули, прогнали двас.

Точек на одной кривой - много. Несколько десятков тысяч. Да и самих таких кривых не одна сотня... Лучше, конечно, все сделать в автомате, дабы пользователи не парились сильно (насколько это возможно).

Кстати, на работу резонансного подавителя не повлияет тот факт, что на каждой кривой период 50-герцовки может отличаться? Т.е быть не 20 мс, а 20.5 или 19.2 мс?