Проектирование оконных функций, суммирующихся в единицу с заданным уровнем перекрытия23.11.2018 11:32

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

Сам процесс разбиения математически представляется умножением на некоторую весовую (оконную) функцию со смещением. Для самого простого окна — прямоугольного — это может выглядеть так:

Исходный сигнал:

Разбиения:

Можно восстановить исходный сигнал, просто просуммировав их.

Подробнее

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

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

Для устранения этих недостатков используется перекрытие — когда каждое следующее окно захватывает часть данных из предыдущего;, а весовое окно, соответственно, плавно спадает к краям.

Перекрытие в 50%

Наиболее часто для этого используют окно Ханна (известное также под названием «приподнятый косинус») с перекрытием в 50%:

За счёт симметричности функции косинуса при сложении они суммируются в единицу:

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

Обработка с двойным наложением окон

Рассмотрим более подробно алгоритм для обработки сигнала с использованием быстрого преобразования Фурье (FFT):

разбиение сигнала на сегменты с наложением окна;
прямое FFT;
обработка спектра;
обратное FFT;
повторное наложение окна (поскольку после обратного FFT границы не обязательно будут стыковаться с нулём без разрыва);
суммирование результирующих сегментов.

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

При использовании окна Ханна перекрытия в 50% уже недостаточно, так как будут возникать провалы:

Поскольку двойное наложение равносильно возведению в квадрат, то очевидным решением будет использовать корень из окна Ханна, чтобы компенсировать возведение в квадрат:

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

Примечание

Интересно, что в данном случае у нас получилась половина периода синусоиды.

Можно пойти и другим путём — использовать перекрытие в 75%, и тогда окна также будут суммироваться в константу — только уже не в $1$ , а в $\frac{3}{2}$ :

В данном случае нам просто повезло. Если разложить квадрат на сумму — то можно увидеть, что оно представляет собой композицию из двух окон Ханна, но в разных масштабах, что и позволяет выполнять необходимые нам требования:

$\left(\frac{\cos (2 \pi x)+1}{2} \right)^2 = \frac{\cos (2 \pi x)+1}{2} + \frac{\cos (4 \pi x)-1}{8}$

С другими оконными функциями такой фокус не пройдёт; также далеко не все стандартные оконные функции могут обеспечить суммирование в константу и при 50% перекрытии.

Идея

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

$\left\{ \begin{array}{ll} -1 & x\leqslant -1 \\ 1 & x\geqslant 1 \\ \sin \left(\frac{\pi x}{2}\right) & -1 < x < 1 \\ \end{array} \right.$

можно записать

$f(x+1)-f(x-1)$

Рассмотрев компоненты суммы двух таких окон, их способность суммироваться в константу станет более наглядным:

На отрезке $[0,2]$ мы получили взаимную компенсацию при сложении функций $-f(x-1)$ и $f((x+1)-2)$ , поскольку $f((x+1)-2)=f (x-1)$

Можно выбрать и другое смещение, с большим перекрытием, например:

$f\left(x+\frac{1}{2}\right)-f\left(x-\frac{1}{2}\right)$

И тогда, при смещении с шагом $\frac{1}{2}$ , они также будут суммироваться в константу:

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

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

Итоговая формула

Теперь осталось только задать масштаб, чтобы области определения и значения не зависели от величины перекрытия. Определив окно на интервале $[0,1]$ получим

$\frac{f\left(\frac{2 t x}{t-1}-1\right)-f\left(\frac{2 t (x-1)}{t-1}+1\right)}{2}$

где $f$ — кусочно-непрерывная функция вида

$\left\{ \begin{array}{ll} -1 & x\leqslant -1 \\ 1 & x\geqslant 1 \\ g(x) & -1 < x < 1 \\ \end{array} \right.$

а $g$ — некоторая интерполирующая функция между точками $(-1,-1)$ и $(1,1) $ .

Параметр $t$ определяет уровень перекрытия — делитель, определяющий шаг, на который необходимо смещать каждое следующее окно, $x_{n+1}=x_n+\frac{1}{t}$ , и должен быть больше единицы, $t > 1$» />.</p>

<p>Процент перекрытия <img src= можно посчитать по формуле

$p=\frac{100 (t-1)}{t}$

и соответственно наоборот

$t=\frac{100}{100-p}$

Примечание

При работе с реальными может потребоваться пересчёт уровня перекрытия в зависимости от конкретного размера FFT. Например, при FFT на 2048 точек и уровнем перекрытия 3 получаем шаг 2048/3=682.666…, что на практике, естественно, нереализуемо. Поэтому округлим его до целого, а $t$

пересчитаем как 2048/683=2.998535871156662.

Либо же можно наоборот — использовать размер окна, заведомо делящийся на 3 (скажем, 999), а остаток до необходимого размера для FFT дополнить нулями (25-ю).

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

Полиномиальные окна

Являются наименее вычислительно затратными. Используя ранее выведенную формулу

$\frac{2 x \Gamma \left(n+\frac{1}{2}\right) \, _2F_1\left(\frac{1}{2},1-n;\frac{3}{2};x^2\right)}{\sqrt{\pi } \Gamma (n)}$

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

первые 10 полиномов

$\begin{array}{c} x \\ \frac{1}{2} x \left(3-x^2\right) \\ \frac{1}{8} x \left(3 x^4-10 x^2+15\right) \\ \frac{1}{16} x \left(-5 x^6+21 x^4-35 x^2+35\right) \\ \frac{1}{128} x \left(35 x^8-180 x^6+378 x^4-420 x^2+315\right) \\ \frac{1}{256} x \left(-63 x^{10}+385 x^8-990 x^6+1386 x^4-1155 x^2+693\right) \\ \frac{x \left(231 x^{12}-1638 x^{10}+5005 x^8-8580 x^6+9009 x^4-6006 x^2+3003\right)}{1024} \\ \frac{x \left(-429 x^{14}+3465 x^{12}-12285 x^{10}+25025 x^8-32175 x^6+27027 x^4-15015 x^2+6435\right)}{2048} \\ \frac{x \left(6435 x^{16}-58344 x^{14}+235620 x^{12}-556920 x^{10}+850850 x^8-875160 x^6+612612 x^4-291720 x^2+109395\right)}{32768} \\ \frac{x \left(-12155 x^{18}+122265 x^{16}-554268 x^{14}+1492260 x^{12}-2645370 x^{10}+3233230 x^8-2771340 x^6+1662804 x^4-692835 x^2+230945\right)}{65536} \\ \end{array}$

С перекрытием в 75% окна с различными значениями n окна будут выглядеть так:

А в случае с извлечением корня — так (при использовании 75% перекрытия):

или так (при использовании 50% перекрытия):

Максимально гладкое окно

Функция

$\tanh \left(\frac{k x}{\sqrt{1-x^2}}\right)$

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

Окно вида «юбочка»

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

Сначала определим желаемый вид оконной функции — например, следующим образом:

$-\log \left(k^2 x^2+1\right)+\log \left(k^2+1\right)-\frac{k^2 \left(1-x^2\right)}{k^2+1}$

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

$\frac{k x \left(k^2 \left(x^2+3\right)+6\right)+3 \left(k^2+1\right) \left(k x \left(\log \left(k^2+1\right)-\log \left(k^2 x^2+1\right)\right)-2 \tan ^{-1}(k x)\right)}{4 k^3-6 \left(k^2+1\right) \tan ^{-1}(k)+6 k}$

Для удобства настройки можно привязать параметр $k$ к уровню перекрытия $t$ — например, таким образом, чтобы 4-ая производная в центре окна была равна 0 — и тогда $k$ будет считаться как $\sqrt{3} (t-1)$ :

Здесь для наглядности все окна приведены к одному масштабу.

Окно вида «иголка»

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

$\frac{1}{x}\to \frac{1}{\sqrt{x^2}}\to \frac{1}{\sqrt{x^2+1}}\to \frac{1}{\sqrt{k^2 x^2+1}}\to \frac{\left(1-x^2\right)^2}{\sqrt{k^2 x^2+1}}$

и использования тех же шагов в виде интегрирования и масштабирования получили формулу

$\frac{k x \left(2 k^2 \left(x^2-4\right)-3\right) \sqrt{k^2 x^2+1}+\left(8 \left(k^4+k^2\right)+3\right) \sinh ^{-1}(k x)}{\left(8 \left(k^4+k^2\right)+3\right) \sinh ^{-1}(k)-3 k \sqrt{k^2+1} \left(2 k^2+1\right)}$

Здесь также можно привязать параметр $k$ к уровню перекрытия. Непосредственное решение уравнения 4-ой производной даёт слишком громоздкий результат, поэтому просто сделаем по образу и подобию из решения для предыдущего окна, определив $k$ как $k(t-1) $ , обеспечив таким образом роль параметра $k$ как «тонкая настройка». При $k=2.22$ окна будут выглядеть следующим образом:

Асимметричное окно

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

$(1-x)^2 \left(1-x^{10}\right)^2$

Здесь степень при x в веcовом окне (а именно 10) определяет «резкость» атаки. Мы использовали конкретное значение, а не символьный параметр, чтобы упростить формулы для наглядности — при желании, в дальнейшем его можно пересчитать.

После интегрирования просто масштабирования уже недостаточно — нужно ещё выровнять края:

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

$\frac{8775 \left(\frac{x^{27}}{27}-\frac{2 x^{26}}{26}+\frac{x^{25}}{25}-\frac{2 x^{15}}{15}+\frac{4 x^{14}}{14}-\frac{2 x^{13}}{13}+\frac{x^3}{3}-x^2+x+\frac{117592}{61425}\right)}{9856}-1$

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

Заключение

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

Вне рассмотрения остался спектральный состав таких оконных функций — этому нужно посвящать отдельные исследования.

Чуть более расширенный вариант статьи (с возможностью динамического изменения параметров в рассматриваемых окнах и скрытыми формулами) в виде документа Wolfram Mathematica можно скачать здесь.