Вывод функции кривой для плавного ограничения параметров, сигналов и не только в Wolfram Mathematica

Существует ряд задач, в которых диапазон выходных значений должен быть ограничен, в то время как входные данные этого гарантировать не могут. Помимо вынужденных ситуаций, ограничение сигнала может быть и целенаправленной задачей — например, при компрессии сигнала или реализации эффекта «overdrive».

Самая простая реализация ограничения — это принудительная установка в некоторое значение при превышении определённого уровня. Например, для синусоиды с возрастающей амплитудой это будет выглядеть так:

ktkxik-cgxbfy0meutgbvtjk3jw.png

В роли ограничителя здесь выступает функция Clip, в качестве аргумента которой передаётся входной сигнал и параметры ограничения, а результатом функции является выходной сигнал.

Посмотрим на график функции Clip отдельно:
0cr8_bhnjcfs5n5w7gxgtthrt3k.png

Из него видно, что пока мы не превышаем пределы ограничения, выходное значение равно входному и сигнал не меняется; при превышении же выходное значение от входного уже никак не зависит и остаётся на одном и том же уровне. По сути, мы имеем кусочно-непрерывную функцию, составленную из трёх других: y=-1, y=x и y=1, выбираемых в зависимости от аргумента, и эквивалентную следующей записи:

vb_fkn3xitnusafjnlllo26zdae.png

Переход между функциями происходит довольно резко; и выглядит заманчивым сделать его более плавным. Математически эта резкость обусловлена тем, что производные функций в точках стыковки не совпадают. Это легко увидеть, построив график производной функции Clip:

68wvo7ts5hrfff5you8ksjtc4cg.png

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

Синус


Самое простое — это использовать функцию sin на интервале от -pi/2 до pi/2, на границах которого значения производной равны нулю по определению:

egdxvlm7_qwzyrmxqezv3adpcig.png

Нужно только масштабировать аргументы, чтобы единица проецировалась на Pi/2. Теперь мы можем определить собственно ограничивающую функцию:

r-fzxh68on-ummingwg4nghun-w.png

И построить её график:

-wso0al-dlmqf459muqlxbdxk3y.png

Так как пределы ограничения у нас жёстко определены, то ограничение задаётся через масштабирование входного сигнала с последующим (при необходимости) обратном масштабировании.
Здесь также уже нет ситуации, при которой входной сигнал передаётся на выход без искажений — чем меньше уровень усиления, тем меньше уровень искажений вследствие ограничения —, но сигнал искажается в любом случае.
Влияние параметра усиления на искажение сигнала можно посмотреть и в динамике:

fvzd96eakididqx7xiht7be2nbi.png

-bie8fl15hlbrawxba6yhxdgh-0.gif

Больше гладкости


Посмотрим на производную нашей функции:

wvkki7mhi12fscostsrzmjiao9e.png

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

xdclgbcrxex2aomgr0u8aornfkg.png

Находим первообразную:

93lxjg9ea2wlksrro8bpc0crrws.png

Теперь осталось масштабировать её по оси ординат. Для этого найдём её значение в точке 1:

djh2hvrq3uiqm_zfhlx9odkdhww.png

И поделим на неё (да, конкретно здесь это элементарное умножение на 2, но далеко не всегда так бывает):

l2tgqv3pzu-qhqr6mi33fvtyufg.png

Таким образом, итоговая функция ограничения примет вид:

pk1hotuczx5fhsal6fexcaz-guo.png

Переходим на полиномы


Использование тригонометрический функций в некоторых случаях может оказаться несколько расточительным. Поэтому попробуем построить необходимую нам функцию, оставаясь в рамках элементарных математических операций.
Рассмотрим параболу:

ytz9dpncw4srtpxv_ybck_jkcvg.png

Так как у неё уже есть перегиб в точке ноль, мы можем использовать одну и ту же часть на интервале {0,1} для для стыковки с константами. Для отрицательных значений её нужно сместить вниз и влево:

j9lbunmhg7tb40vuwvlour2hfza.png

а для положительных — отразить по вертикали и горизонтали:

p2wus68r_fgkqce_j3gbmpkfikm.png

И наша функция с параболой примет вид:

b42ynjrrmblafqcitajzap0aajk.png

Немного усложним


Вернёмся к нашей параболе, перевернём её и сместим на единицу вверх:

zc0pak9jwjocs9upehf8fxqmrxe.png

Это будет производная нашей функции. Чтобы сделать её более гладкой в точках стыковки, возведём в квадрат, обнулив таким образом вторую производную:

3kb484zji-8wtjlz6csjwm-nh-w.png

Интегрируем и масштабируем:

z9ekzqvt2oqsltttiz6sbno3ojk.png

Получаем ещё более гладкую функцию:

vfcjr9bcxm_miwuhgesvvpzbswm.png

Больше гладкости богу гладкости


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

Начнём со 1-й производной:

ohjuqfq1s8miyxy4s5dxvnszjey.png

2-я:

zdgl1zrxrnju1dmy9untv3tfxlo.png

3-я:

h6ppvpclnin5_am8knkm72nvxci.png

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

i_5eqbkf288h6ttqnc-iezn11rq.png

Похоже на биномиальные коэффициенты. Сделаем смелое предположение, что это они и есть, и исходя из этого, запишем обобщённую формулу:

uie3te4qo2puxc-b99skogktfay.png

Проверим:

uyumgc9s-15ue6ki3gnrgvms1de.png

Похоже на правду [1]. Осталось только посчитать масштабный коэффициент, чтобы привести края к единице:

ir4x4ib21qyi2mrrnj_tkoa7f6o.png

А после масштабирования и упрощения мы обнаружим, что наши познания в математике несколько устарели [2]:

tgxjgpjxkfvkzmtku4kqgxkibyu.png

Таким образом, мы получили производящую функцию порядка n, в которой n-1 первых производных будут равны нулю:

ifcsl5x8nvomwcy2f0gdgkfekok.png

Посмотрим, что получилось:

jwng1sr__oi8ybymbij4glzz2hc.png

И поскольку наша обобщённая формула получилась непрерывной, при желании можно использовать и нецелые значения параметров:

cldy2azeu3nbwn6ynnw4eidltcq.png

Также можно построить графики производных, приведённых к одному масштабу:

keckxnnpj3ed53gokz5k9uplume.png

Добавляем жёсткости


Было бы заманчиво, иметь возможность регулировать и степень «жёсткости» ограничения.
Вернёмся к нашей перевёрнутой параболе и добавим коэффициент при степени x:

mrut-asek8r6qxdf2ghmwmh5f0u.png

Чем больше n, тем больше наша производная «квадратная», а её первообразная — соответственно, резкая:

9nfv5dyxws_yrc09vlhifrgj_pi.png

Посчитаем первообразную и скорректируем масштаб:

y_hb01xyoqzqt8favbkjxetljxq.png

Попробуем теперь задать дробный шаг для параметра:

zzr_efa81bti9lgqij8nzxffx6s.png

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

c8oeudj7c9ykakyrb6obkfsu0jg.png

А заменив n на n-1, можно сделать формулу чуть более красивой:

zd8fhm6o_kufvtfnyqotb37jdfq.png

Поскольку при n равным единице мы получаем деление на ноль, то попробуем найти предел:

8j_3tajousxrcx6lbl4qzl5_gbu.png

Предел находится, а значит, теперь можно доопределить [3] функцию для n равным 1 и рассматривать её для всех n больших нуля:

1r5fye54zeqcywmgpyrs9x_sp10.png

Если же мы изначально возведём нашу перевёрнутую параболу в квадрат, то получим ещё более гладкую функцию:

4yoewctdmaoxo1sr3ii_1oc1ioi.png

И можем сравнить их на одном графике:

guy6humvdmuklaxd4sgvkp60nso.png

Рационализируй это


Посмотрим на следующую функцию:

r_6dpot4wilqhjmoacqvxend44c.png

Появилась она не случайно.
Если убрать из неё единицу, x2 сократится и останется просто x, т.е наклонная прямая. Таким образом, чем меньше значение x, тем большее влияние оказывает единица в знаменателе, создавая необходимое нам искривление. А рассматривая эту функцию в разных масштабах, можно контролировать степень этого искривления:

32tgvq33kppw2mykapjvbv8ulim.png

Таким образом, мы можем переписать предыдущую функцию с контролем жёсткости, используя только рациональный полином 3-порядка:

_emlex4lxbl2vo5vrar7qmy18po.png

Автоматизируй это


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

Если наша функция уже обладает диагональной симметрией и выровнена по центру координат (как синусоида), то можно сделать просто

v7ynzoykxwrgavh9st26zozjlhc.png

Пример использования:

cucpsu5zkhqy0rejqchvdb49jye.png

Если же нужно собирать из кусочков, как в случае с параболой, и центр координат определяет точки стыковки, то формула слегка усложнится:

g2des_d0oozeon_eeanhwnn0ulc.png

Пример использования:

ywndf6wis3bxqsqk1yqwlxjfsek.png

Перейдём на экспоненту


Совершенно любая функция может быть донором для решения этой задачи, нужно лишь только обеспечить ей точки перегиба. Возьмём, например, сдвинутую вниз на единицу экспоненту:

6eqp8g4hecwf3bg9yvij5-7dkka.png

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

zn-41ubtnnsp4cgzipkk_3qhy-e.png

В зависимости от того, с какой стороны мы будем брать кривую, будет и зависеть итоговый вид функции. Теперь, воспользовавшись ранее определённой вспомогательной функцией и выбрав одну из сторон, получим:

9kugljxkkn4wqnxonhkx5fnxd08.png

Либо

gfl0agc1n4kziemm02nryxkap-y.png

И теперь можем сравнить их на одном графике:

xyeq2z9j2sdpppoahp9z93u_lck.png

Видно, что при k→0 они стремятся к совпадению; и так как напрямую посчитать их значения мы не можем, поскольку получим деление на ноль, то воспользуемся пределом:

xveymim695sl6dkkwmiy6plyx1k.png

И получили уже известную нам кусочную функцию из параболы.

Нарушая симметрию


Пока что мы рассматривали исключительно симметричные функции. Однако бывают случаи, когда симметрия нам не нужна — например, для имитации искажений при звучании ламповых усилителей.

Возьмём экспоненту и умножим её на перевёрнутую параболу в квадрате — чтобы получить пересечение с осью абсцисс в точках -1 и 1, а заодно и обеспечить гладкость второй производной; параметризацию же осуществим через масштабирование аргумента экспоненты:

3s0-zhgrrr_exzjqdphozhaentm.png

Найдём первообразную и масштабируем её:

hbbsaddm5vfi-r-u3cxwi_ul6x8.png

Так как при k=0 получим деление на ноль:

zdq15cnt6vkhmtbqrw-svlnagca.png

То дополнительно найдём предел,

va4d_mi7nqeqqyfb0rtwat0j7ru.png

который представляет из себя уже известный нам гладкий полином 3-го порядка. Соединив всё в одну функцию, получим

-y1nhdpvrlk2syka-f2opst6sbs.png

Вместо того, чтобы изначально проектировать асимметричную функцию, можно пойти и другим путём — использовать готовую симметричную, но «искривлять» значение этой функции с помощью дополнительной функции кривой, определённой на промежутке {-1,1}.

Рассмотрим, например, гиперболу:

-15eskytmag256r7fykt5dyv0xq.png

Рассматривая её отрезок в разных масштабах, можно регулировать степень искривления в обе стороны. Как же найти этот отрезок? Исходя из графика, можно было бы искать пересечения гиперболы с прямой. Однако, поскольку такое пересечение существует не всегда, это создаёт некоторые сложности. Поэтому мы пойдём другим путём.

Для начала добавим в гиперболу масштабирующие коэффициенты:

rmlbptoswocefui6xyeki6hqj-o.png

затем составим систему уравнений, задающих условия прохождения гиперболы через заданные точки — и её решение даст интересующие нас коэффициенты:

y-hp72vr9u_lcagpd8wjnoffyta.png

Теперь подставим решение в исходную формулу и упростим:

sommfjjypnhaq0mysrt_fg-hsm0.png

Посмотрим, что у нас получилось в зависимости от параметра k:

55ymh2cr6mhi1ugki_fuyz5jnsy.png

Примечательно, что при k=0 формула естественным образом схлопывается в x и никаких особых ситуаций не происходит — хотя применительно к исходной гиперболе это равносильно отрезку нулевой длины, причём двум сразу. Не менее примечательно, что обратной к ней функцией является она же самая, но с отрицательным параметром k:

ffhbme8oa8xdcofo93brdapjsw4.png

Теперь мы можем использовать её для модификации произвольной функции ограничения, а параметр k таким образом будет задавать точку пересечения с осью ординат:

6s0cr-iagfq5xgczfyi_fwucpgm.png

Аналогичным образом можно строить кривые и из других функций, например, степенной с переменным основанием:

6_ckox4togzkd05i3024v_cvgme.png

Или обратной к ней логарифмической:

gxlxepkauvevwv9qjihc1pheerg.png

Нужно больше точности


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

ep-4muzhdl4qxge1z-jssrvknrw.png

пустые места в которой необходимо заполнить какой-нибудь функцией. Очевидно, что для гладкой стыковки с линейным участком её первая производная должна быть равна единице;, а все последующие (по возможности) нулю. Чтобы не не выводить такую функцию заново, мы можем взять уже готовую и адаптировать под эту задачу. Также можно заметить, что крайние точки отстоят чуть дальше единицы — это необходимо, чтобы сохранить наклон линейного участка.

Возьмём выведенную ранее функцию PolySoft и сместим её так, чтобы в центре координат получить единицу:

azt5al4l8j-gfh-7offgynyx_rg.png

Из её свойств следует, что n-1 последующих производных в точках 0 и 2 будут равны нулю:

rx5lygmui7me5svsy6wyezkbdcc.png

Теперь проинтегрируем её:

b0w2ydyxd9h1nolap6zdakvy9be.png

Функция оказалась сдвинутой вниз относительно оси абсцисс. Поэтому необходимо добавить константу (равную значению функции в точке 0), чтобы совместить центры координат:

uyhjd6sebrc0thj0fstiqpryz3s.png

Здесь у нас появился ноль в степени n. Он не сократился, так как значение ноль в степени ноль не определено; мы можем его удалить вручную, а можем при упрощении явно указать, что n у нас больше нуля:

a6fmraosk7kdnypmqe9xsw2s9au.png

Проверим на всякий случай. Значение в точках 0 и 2 для всех n:

fkdbzkwa6jcgnny7bzlpdnolixm.png

Производные на краях интервала (для полинома порядка 5):

_efhmmh8nr4rvdbf3u921u_qzeq.png

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

1yjegsufbn1ib1vvonv7zzsi82s.png

И вот теперь ею можно заполнить свободное пространство:

x6kyrwgaf47bjaig9s5k_5ll_5k.png

Проверим:

ejajbo72zeabhiilyuoug8txg30.png

Уходим в бесконечность


Иногда может оказаться потребность в функциях, которые стремятся к единице, но не достигают её. Википедия подсказывает несколько известных решений:

puhj4hrsljuwbuqrbkjqrn17qng.png

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

q82jtyy0vc3thqoqz-un0bi23di.png

Эта функция, к слову, также является обратной самой к себе, т.е.
speahbxon7pelj37jliv-2c3k-4.png

И, применительно к арктангенсу в качестве примера, получим

qfgkkqarxwqn-qh5a65xhktemvk.png

что, в частности, с параметром k=1 даст нам функцию Гудермана.

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

Из степенной функции:

ssfz76sycvuu5gtilb4r3rv63rg.png

Из суммы двух v-образных функций со смещением:

e3qunyzxsvzuuwmcce24plj-bgq.png

Из обобщённой функции ошибок:

h3daxpyc8i67bcgy_k32pz5bbsm.png

Интегрированием рационального полинома:

snglw1xjo0riujdtna8m9gr004a.png

Интересно, что её частным случаем является арктангенс:

flvtzk2kvy-i9yvug9qrf75eibe.png

Заключение


Построение подобного рода функций может быть увлекательным занятием, в ходе которого будут получаться как простые, так и сложные, как красивые, так и не очень, формулы. Может показаться, что все они сильно друг на друга похожи и надобности в подобном разнообразии нет. Это не обязательно так.

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

В завершение стоит уточнить ещё несколько моментов.
Все функции здесь были определены в диапазоне от -1 до 1. В случае, если нужен другой диапазон (например, от 0 до 1), его легко можно пересчитать либо вручную:

lgnr1d4q95ieqwb-151g3manyio.png

Либо используя встроенную функцию масштабирования:

6iyoujbtrj5peypvpccflthjals.png

А для облегчения экспорта полученных формул в программный код может пригодиться функция CForm:

jcgaabfe5b5ooc-a6a9_mw_kkq8.png

Исходный документ Mathematica можно скачать здесь.

Примечания:

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

© Habrahabr.ru