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

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

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



Разбиения:







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

Подробнее


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



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



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

Перекрытие в 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 можно скачать здесь.

Комментарии (12)


  1. Sabubu
    23.11.2018 11:38

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

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

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

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

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

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


    1. Refridgerator Автор
      23.11.2018 12:45

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

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

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



    1. Refridgerator Автор
      23.11.2018 12:57

      И еще вопрос, может вы знаете ответ, раз уж занимаетесь обработкой сигналов. Если есть сигнал неизвестного вида, который может содержать частоту X + шумы и щелчки, и мы хотели бы получить на выходе амплитуду этой частоты X в сигнале, то есть увидеть, где она есть, где ее нет, где ее больше — есть какой-то готовый простой алгоритм, без вычисления FFT? Какой-то простой фильтр, может быть? Нужно для определения задержки прохождения сигнала через звуковую карту и ее драйвера.
      На мой взгляд, вопрос некорректно поставлен. Для определения задержки прохождения сигнала через звуковую карту и ее драйвера не нужно анализировать амплитуду частот, нужно обеспечить точку отсчёта и прогонять единичный импульс или ЛЧМ. Задача достаточно сложная, чтобы можно было дать на неё исчерпывающий ответ в комментарии.


      1. Sabubu
        23.11.2018 17:46

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


        1. Nokta_strigo
          23.11.2018 22:20

          А чем FFT не устраивает? Медленно считается? Если нужна только одна частота (т.е. одна точка на дискретном спектре), то можно считать дискретное преобразование Фурье для одной частоты. Сложность FFT — O(n * log (n)), а расчёта одной точки спектра «в лоб» — O(n), т.е. для очень узких фильтров (при больших n) будет в несколько раз быстрей.
          Но учтите, что чем точней вам надо получить частоту, тем более длинную реализацию сигнала (окно) надо брать. Это принципиальное ограничение, не важно используете FFT (хотя ресурсов на расчёт FFT тоже надо больше), «железные» фильтры или ещё что.
          Так что узкополосный сигнал может быть не самым удачным выбором, если надо измерять время. Как написал Refridgerator, лучше мощный импульсный сигнал или какой-то более сложный с широким спектром (ЛЧМ или ещё что).


    1. Refridgerator Автор
      23.11.2018 13:08

      И интересно, где вообще нужны нестандартные оконные функции, особенно несимметричной формы?
      При нелинейной обработке сигнала через FFT результат сильно зависит от оконной функции. Если после обработки возникает артефакт в виде реверберации или «звона», то несимметричная форма окна позволяет его сфокусировать в правой части импульса, что может дать более естественное звучание. Отличия те же, что и между фазолинейными и минимально-фазовыми фильтрами.


    1. Refridgerator Автор
      23.11.2018 13:40
      +1

      Например, с их помощью я делал генератор шума дождя — из записи дождя ограниченной продолжительности извлекал небольшие куски со случайным смещением и соединял их заново (в 4-х каналах), чтобы не слушать бесконечно одно и то же и получить объёмный звук. Это позволяло также настраивать «плотность» звука, смешивая в разных пропорциях разные части записи.


    1. interrupt
      24.11.2018 00:50

      По поводу определения неизвестной частоты в сигнале, возможно Generalized Harmonic Analysis то что вы ищете, например статья «Fast and Accurate Generalized Harmonic Analysis Using Newton’s Method». Ну и пользуясь случаем порекламирую библиотечку, которую как раз начал писать по этой статье, libgha
      Идея там такая, сначала через FFT грубо ищем максимум, потом используя метод Ньютона уточняем частоту и фазу.
      Работает это в предположении что в сигнале нет двух синусоид частоты которых расположены ближе частотного разрешения FFT с заданным окном.


    1. blackghost56
      24.11.2018 15:21

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

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


  1. maybe_im_a_leo
    23.11.2018 22:19

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


    1. Refridgerator Автор
      24.11.2018 15:17

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


  1. ser-mk
    23.11.2018 23:27

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

    А какое сложение используется?
    Или имеется ввиду, что будет единица от первого «бугра» до последнего?