Skip to main content

Гладкая сшивка кусочно-постоянной функции.

Накануне мне понадобилось найти гладкую функцию, определенную на $\textbf{R}[a,c]$, равную единице на одном отрезке числовой оси, скажем, $[a, b-\epsilon]$ и нулю на другом, скажем $[b+\epsilon, c]$. Для того, чтобы гладко склеить эти два кусочка, я решил найти монотонно убывающий многочлен $f(x)$, равный $f(0) = 1$ и $f(1) = 0$, а потом "вставить" его в отрезок $[b-\epsilon, b+\epsilon]$. Для гладкости сшивки должны выполняться условия $f'(0) = f'(1) = 0$. Я решил пойти далее и положить $f^{(n)}(0) = f^{(n)}(1) = 0$ для всех $n = 1,2,\cdots,N$. То есть первые N производных должны быть равны нулю в т. $x = 0$ и $x = 1$.

Понятно, что это легко сделать при наперед известном N, но как найти общую формулу? Я некоторое время промучился над этой задачей, пока не решил её следующим способом.

Пусть искомая функция $f(x) = 1 + x^{N+1}L_{1}(x)$, где $L_{1}(x)$ — некий многочлен. Тогда у нас выполняется условие $f(0) = 1$ и $f^{(n)}(0) = 0$ (так как первая производная будет иметь нуль порядка $N$).

Пусть теперь та же функция равна $f(x) = (x-1)^{N+1}L_{2}(x)$ мы получаем выполнение условий в т. $x = 1$

Приравниваем эти два определения и получаем уравнение:
$1 + x^{N+1}L_{1}(x) = (x-1)^{N+1}L_{2}(x)$

Положим теперь для простоты N — нечетное, тогда уравнение можно переписать в виде:

$x^{N+1}L'_{1}(x) + (1-x)^{N+1}L_{2}(x) = 1$, где $L'_{1}(x) = -L_{1}(x)$

Заметим, что наибольший общий делитель у многочленов $(1-x)^{N+1}$ и $x^{N+1}$ равен единице, а значит уравнение выше есть ни что иное, как уравнение Безу. Есть алгоритм его решения — расширенный алгоритм Евклида, однако именно для этого случая известен ответ в явном виде (из теории вейвлетов). Перепишем уравнение следующим образом:

$x^{N+1}P(1-x) + (1-x)^{N+1}P(x) = 1$
Одним из решений (с наименьшей степенью) этого уравнения будет многочлен $P(x) = \sum_{n=0}^{N}C_{N+n}^{n}x^n$

Теперь найдем $f(x)$ пользуясь первым определением. Получаем $L_{1}(x) = -L'_{1}(x) = -P(1-x)$. Тогда $f(x) = 1 - x^{N+1}\sum_{n=0}^{N}C_{N+n}^{n}(1-x)^n$
Это и есть ответ. Замечу, что ничего не изменится, если N — четное, только решение уравнения Безу в явном виде придется поискать.

UPD 25 jan 2019:
Приведенное решение, оказывается, годится и для четного N. Для того, чтобы в этом убедится достаточно представить $f(x)$ в следующем виде:
$f(x) = 1 - x^{N+1}P(1-x) = (1-x)^{N+1}P(x)$
Конец апдейта

Приведу пару примеров. Для $N = 1$ получаем $f(x) = 2x^3 - 3x^2+1$.
Для $N = 3$ получаем $f(x) = 20x^7 - 70x^6 + 84x^5-35x^4+1$.
Вот график этих функций:
Чего тут нет — так это доказательства монотонности (т.е. нет локальных максимумов и минимумов внутри интервала). Ну и Б-г с ним.

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

Так же как и в первый раз даем два определения $f(\omega)$:
$f(\omega) = 1 - (\frac{1 - \cos \pi \omega}{2})^N Q_{1}(\cos \pi \omega)$ и
$f(\omega) = (\frac{1 + \cos \pi \omega}{2})^N Q_{2}(-\cos \pi \omega)$

Первое определение удовлетворяет условию слева ($x = 0$), второе ­— справа ($x = 1$).

Вспоминаем, что $\cos 2x = \cos^2 x - \sin^2 x = 1 - 2\sin^2 x = 2\cos^2 x -1$, откуда получаем $\sin^2 x = \frac{1-\cos 2x}{2}$, $\cos^2 x = \frac{1+\cos 2x}{2}$

Обозначим $x = \sin^2 \frac{\pi}{2} \omega$ и перепишем определения $f(\omega)$ в термах x, приравняв их друг к другу.

$1 - x^N Q_{1}(1-2x) = (1-x)^N Q_{2}(2x-1)$ или
$x^N Q_{1}(1-2x) + (1-x)^N Q_{2}(2x-1) = 1$

Пусть теперь $P(x) = Q_{2}(2x-1)$, а $P(1-x) = Q_{1}(1-2x)$. По сути, $Q_1$ и $Q_2$ — одни и те же многочлены, обозначены они по-разному из-за того, что мы не знаем этого заранее.

Получаем уравнение:
$x^N P(1-x) + (1-x)^N P(x) = 1$

Решение мы знаем:
 $P(x) = \sum_{n=0}^{N-1}C_{N+n-1}^{n}x^n$

Из $P(x)$ находим $Q_{1}(x)$ и $Q_{1}(\cos \pi \omega)$.

Пару примеров. При $N = 1$ получаем тривиальное $f(\omega) = \frac{1+\cos \pi \omega}{2}$.

При $N = 2$ получаем $P(x) = 1+2x$, $P(1-x) = 1+2(1-x) = 3 - 2x$, а $Q_{1}(1-2x)  = 2 + (1-2x) = 2 + \cos \pi \omega$.

Тогда $f(\omega) = 1 - (\frac{1- \cos \pi \omega}{2})^{2} (2 + \cos \pi \omega) = \frac{1}{16}(-\cos 3 \pi \omega + 9 \cos \pi \omega + 8)$

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

Comments

Popular posts from this blog

Ресемплинг аудио.

В этом посте расскажу о ресемплинге аудио (в частности, о даунсемплинге), под чем я подразумеваю понижение частоты дискретизации, с которой был записан звук. Известно, что в компьютере непрерывный аудио сигнал представлен в виде семплов (по-русски, измерений), сделанных в равно отстоящие друг от друга промежутки времени. Например, если звук записан с частотой дискретизации 44.1 кГц, то за одну секунду было сделано 44100 семплов. Здесь я расскажу о том, как восстановить исходный непрерывный сигнал по взятым семплам и как выбирать и менять частоту семплирования (или, что то же самое, частоту дискретизации). Преобразование Фурье. Преобразованим Фурье $F(\xi)$ функции $f(x)$ называется следующая функция: $F(\xi) = \int_{-\infty}^{\infty} f(x) e^{-ix\xi}dx$ обратное преобразование определено как $f(x) = \frac{1}{2\pi}\int_{-\infty}^{\infty} F(\xi) e^{ix\xi}d\xi$ Вопрос о том, когда существует преобразование Фурье и обратимо ли оно (т.е. работает ли вторая формула), очень сложе

На что влияет размер (глубина) семпла?

В этом посте я разберу ещё одну характеристику представления аудио на компьютере — размер или глубину семпла. В теореме Шеннона, с которой мы ознакомились в первом посте, сказано, что семплы нужно брать с определенной частотой, чтобы восстановить исходный сигнал. Однако, в ней ничего не сказано о том, сколько памяти компьютера нужно выделить на один семпл. По сути  ней семплы — это обычные вещественные числа. В компьютере же семплы чаще всего представляются целыми числами в диапазоне от $-2^{n-1}$ до $2^{n-1}-1$ (целые числа со знаком), где $n$ — количество бит на один семпл (bps, bits per sample). С количеством бит на семпл связана другая величина — динамический диапазон (или отношение сигнал-шум, что в данном случае то же самое). Он измеряется в децибелах определяется так: $DR = 20 \lg (2^{n}) \approx 6 n$ Динамический диапазон показывает логарифм отношения самого сильного сигнала к самому слабому. Теперь, что будет, если мы выберем $n$ маленьким? Если взять синусоидальную волн