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

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

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

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



Разбиения:







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

Подробнее


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



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



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

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


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








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



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



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


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

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

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

При использовании окна Ханна перекрытия в 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=\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 можно скачать здесь.
Share post

Comments 20

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

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

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

    И интересно, где вообще нужны нестандартные оконные функции, особенно несимметричной формы?

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

    И еще вопрос, может вы знаете ответ, раз уж занимаетесь обработкой сигналов. Если есть сигнал неизвестного вида, который может содержать частоту X + шумы и щелчки, и мы хотели бы получить на выходе амплитуду этой частоты X в сигнале, то есть увидеть, где она есть, где ее нет, где ее больше — есть какой-то готовый простой алгоритм, без вычисления FFT? Какой-то простой фильтр, может быть? Нужно для определения задержки прохождения сигнала через звуковую карту и ее драйвера.
      0
      пояснений немножечко не хватает. Вы пишете с расчетом на читателя, который прекрасно разбирается в оконных функциях, а наверно, стоило бы немного пояснить подробнее, зачем они нужны, для тех, кто не занимается ежедневно обработкой сигналов и просто заглянул из любопытства
      Пояснения к статье должны быть достаточно краткими для того, чтобы не отвлекать от сути. Я исходил из того, что если читатель не знаком с обработкой сигналов, то статью можно рассматривать просто как математическую задачу с решением.

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

      А то ваши оконные функции похожи на горки и сигнал без изменений остается только в серединке.
      Процент перекрытия может быть и меньше 50%. Я не стал заострять на этом внимание, потому что это казалось очевидным (каюсь):


        0
        И еще вопрос, может вы знаете ответ, раз уж занимаетесь обработкой сигналов. Если есть сигнал неизвестного вида, который может содержать частоту X + шумы и щелчки, и мы хотели бы получить на выходе амплитуду этой частоты X в сигнале, то есть увидеть, где она есть, где ее нет, где ее больше — есть какой-то готовый простой алгоритм, без вычисления FFT? Какой-то простой фильтр, может быть? Нужно для определения задержки прохождения сигнала через звуковую карту и ее драйвера.
        На мой взгляд, вопрос некорректно поставлен. Для определения задержки прохождения сигнала через звуковую карту и ее драйвера не нужно анализировать амплитуду частот, нужно обеспечить точку отсчёта и прогонять единичный импульс или ЛЧМ. Задача достаточно сложная, чтобы можно было дать на неё исчерпывающий ответ в комментарии.
          0
          Ну, это слишком просто и неинтересно. У нас могут быть посторонние шумы и искажения, так как вход и выход не соединены кабелем, а акустической связью через динамики/наушники и микрофон с неизвестными характеристиками.
            0
            А чем FFT не устраивает? Медленно считается? Если нужна только одна частота (т.е. одна точка на дискретном спектре), то можно считать дискретное преобразование Фурье для одной частоты. Сложность FFT — O(n * log (n)), а расчёта одной точки спектра «в лоб» — O(n), т.е. для очень узких фильтров (при больших n) будет в несколько раз быстрей.
            Но учтите, что чем точней вам надо получить частоту, тем более длинную реализацию сигнала (окно) надо брать. Это принципиальное ограничение, не важно используете FFT (хотя ресурсов на расчёт FFT тоже надо больше), «железные» фильтры или ещё что.
            Так что узкополосный сигнал может быть не самым удачным выбором, если надо измерять время. Как написал Refridgerator, лучше мощный импульсный сигнал или какой-то более сложный с широким спектром (ЛЧМ или ещё что).
          0
          И интересно, где вообще нужны нестандартные оконные функции, особенно несимметричной формы?
          При нелинейной обработке сигнала через FFT результат сильно зависит от оконной функции. Если после обработки возникает артефакт в виде реверберации или «звона», то несимметричная форма окна позволяет его сфокусировать в правой части импульса, что может дать более естественное звучание. Отличия те же, что и между фазолинейными и минимально-фазовыми фильтрами.
            +1
            Например, с их помощью я делал генератор шума дождя — из записи дождя ограниченной продолжительности извлекал небольшие куски со случайным смещением и соединял их заново (в 4-х каналах), чтобы не слушать бесконечно одно и то же и получить объёмный звук. Это позволяло также настраивать «плотность» звука, смешивая в разных пропорциях разные части записи.
              0
              По поводу определения неизвестной частоты в сигнале, возможно Generalized Harmonic Analysis то что вы ищете, например статья «Fast and Accurate Generalized Harmonic Analysis Using Newton’s Method». Ну и пользуясь случаем порекламирую библиотечку, которую как раз начал писать по этой статье, libgha
              Идея там такая, сначала через FFT грубо ищем максимум, потом используя метод Ньютона уточняем частоту и фазу.
              Работает это в предположении что в сигнале нет двух синусоид частоты которых расположены ближе частотного разрешения FFT с заданным окном.
                0
                И еще вопрос, может вы знаете ответ, раз уж занимаетесь обработкой сигналов. Если есть сигнал неизвестного вида, который может содержать частоту X + шумы и щелчки, и мы хотели бы получить на выходе амплитуду этой частоты X в сигнале, то есть увидеть, где она есть, где ее нет, где ее больше — есть какой-то готовый простой алгоритм, без вычисления FFT? Какой-то простой фильтр, может быть? Нужно для определения задержки прохождения сигнала через звуковую карту и ее драйвера.

                Посмотрите алгоритм Гёрцеля
                0

                Спасибо, интересная статья. Взял на заметку.
                А в чем вы рисовали графики?

                  0
                  Не только графики, но и все вычисления сделаны в Wolfram Mathematica, в конце статьи есть ссылка на исходник. Код для графиков там скрыт по умолчанию, для открытия нужно выделить ячейку серым цветом над графиком и в меню Cell->Cell Properties->Open поставить галочку.
                  0
                  За счёт симметричности функции косинуса при сложении они суммируются в единицу:

                  А какое сложение используется?
                  Или имеется ввиду, что будет единица от первого «бугра» до последнего?
                    +1
                    Совершенно случайно узнал, что в Vorbis (аудиокодек с потерями) для двойного наложения с 50%-ым перекрытием используется вот такое вот хитрое окно:


                    Для обнуления 1-ой производной они используют модуляцию аргумента. Круто.

                    Не удержался и добавил модуляцию модуляции — и обнулилась не только 2-ая производная, но и 3-я:



                    Правда, преобразование Фурье (аналитическое) от таких окон посчитать не получилось.
                    Wolfram с большим скрипом таки выдал результат на первое окно:

                    А вот второе уже никак.
                      0
                      В ATRAC3 (https://wiki.multimedia.cx/index.php/ATRAC3#Encoding_algorithm) используются разные окна при кодировании и декодировании, можете прокомментировать такое решение?
                        +1
                        Правильно понял, речь об этих окнах?
                          0
                          Да, они самые.
                          Наверно, нужно добавить, там так же используется mdct (как и в большенстве аудиокодеков).
                            +1
                            Для начала нужно привести всё к удобочитаемому виду, так как из кода от FFmpeg вообще ничего не понятно, особенно в строке w = 0.5 * (wi * wi + wj * wj). Вероятно, таким образом авторы пытались оптимизировать вычисления (хотя в этом нет никакого смысла, поскольку вызывается эта функция один раз и тормозить там вообще нечему).
                            код
                            static av_cold void init_imdct_window(void)
                            {
                            	int i, j;
                            	for (i = 0, j = 255; i < 128; i++, j--) 
                            	{
                            		float wi = sin(((i + 0.5) / 256.0 - 0.5) * M_PI) + 1.0;
                            		float wj = sin(((j + 0.5) / 256.0 - 0.5) * M_PI) + 1.0;
                            		float w  = 0.5 * (wi * wi + wj * wj);
                            		mdct_window[i] = mdct_window[511 - i] = wi / w;
                            		mdct_window[j] = mdct_window[511 - j] = wj / w;
                            	}
                            }


                            Первым делом я посмотрел, что рисует (wi * wi + wj * wj), и получилась всё та же синусоида, только удвоенной частоты.Таким образом, если нормировать окно к (-1,1), то получим формулу



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




                            Тут уже хорошо видно, что оно симметричное и с перекрытием в 50% суммируется в единицу.

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

                            Отдельный интересный вопрос, как эта формула выводилась — потому что деление на косинус удвоенной частоты вообще не очевидно. Смещение 3 в знаменателе могло быть найдено решением уравнения

                            Или вообще изначально это была совсем другая формула, которая в итоге сократилась до этой.
                              0
                              Спасибо! Буду обдумывать.

                              А подскажите еще, утверждение:
                              потому что при умножении сигнала на окно происходит свёртка их спектров
                              так же верно и для mdct?
                                0
                                Само собой — оно же происходит до mdct либо любых других преобразований.
                                0
                                На всякий случай хочу добавить, что косинусное окно вовсе не является наилучшим. При таком подходе (разных окнах на вход и выход) на входе можно использовать вообще практически любое окно — а окно на выходе получать делением желаемого результирующего на входное. Например, при использовании окна Блэкмана



                                и того же результирующего окна



                                на выходе получим окно

                      Only users with full accounts can post comments. Log in, please.