Альтернатива оконному Фурье?

[Cos(x-pi/2)]

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

Да, на практике растекания спектра не избежать (ведь крайне маловероятно, чтобы частоты исследуемых сигналов всегда попадали точно в бины FFT), но стоит ли бороться с боковыми лепестками, остающимися после сглаживания окном (или даже без окна)? Так, если сигнал представляет собой сумму стационарных синусоид, заметно различающихся по частоте, как в исходном примере, а задача состоит в нахождении этих частот, то может быть достаточным обычное дискретное пр-е Фурье даже на небольшом количестве отсчётов N и даже без окна (т.е. в прямоугольном окне, выделяющем N отсчётов).

На одном из радиофорумов, в разговоре о точном измерении частоты, когда "знатоки"-любители (к ним по невежеству примкнул было и я :-) уже сошлись во мнении о необходимости долгого накопления отсчётов, один специалист по цифровой обработке сигналов в качестве опровержения дал ссылку на вот такую статью 2007 года:

"Fast, Accurate Frequency Estimators"
E. Jacobsen; P. Kootsookos
https://ieeexplore.ieee.org/document/4205098

Идея там простая: частота тона вычисляется по трём фурье-амплитудам $X_k$ в соседних бинах - около максимума спектрального пика и по обе стороны от него - по формулам (3) и (1) в статье.

Вот, что выходит для примера B@R5uk с частотами $cyc1 = 5.25,$ $cyc2 = 27.75$ с частотой дискретизации $1/dt = num = 280.$ Беру (в Маткаде) в роли сигнала $s$ первые $2^8=256$ отсчётов. В массиве $X=\text{FFT}(s)$ смотрим амплитуды $|X_k|.$ Там виден локальный максимум при $k=5$ и ещё один максимум при $k=25.$ Вычисляем для них значения $\delta_1$ и $\delta_2$ по формуле (3) из статьи:

$\delta =-\text{Re} \left( \dfrac{X_{k+1}-X_{k-1}}{2X_k-X_{k+1}-X_{k-1}} \right).$

У нас:

$X_4=.0475868602+0.055390692i$
$X_5=-0.1762312527-0.2516369141i$
$X_6=-0.027672292-0.0463420793i$

и

$X_{24}=-0.0133080697-0.0328703847i$
$X_{25}=-0.0471254727-0.1119558947i$
$X_{26}=0.0268346206+0.0611161691i$

Получается: $\delta_1 = -0.199793$ и $\delta_2 = 0.372688.$ Вычисляем значения частот по формуле, соответствующей формуле (1) в статье:

$tone=\dfrac{k+\delta}{dt \,256}$

Результат: $tone1 = 5.25023$ и $tone2 = 27.75138.$ В обеих частотах таким простым путём удалось найти четыре верных знака!

(ещё пример с оконным кратковременным Фурье)

Одна из задач, где нет альтернативы оконному Фурье-преобразованию небольших блоков сигнала, притом в реальном времени, - радиоприёмники типа SDR. (Software-defined radio. Поиском по форуму я почему-то не нашёл упоминаний об SDR, поэтому приведу немного подробностей; извините, пожалуйста, если это общеизвестные банальности).

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

По ссылке http://www.websdr.org/ даётся некий список sdr-приёмников, раздающих принятые сигналы в интернет. Первый пункт в списке самый известный - это широкополосный приёмник в университете Twente:

http://websdr.ewi.utwente.nl:8901/

Пояснение тем, кто знакомится с этим SDR впервые: для примера введите там в левом поле частоту 1008 kHz, тип модуляции нажмите АМ, в правом поле Waterfall zoom жмите +, чтобы шкала частот растянулась и появились названия станций. Мышью можно двигать серединку изображения АЧХ демодулятора под шкалой частот - тем самым перестраиваться со станции на станцию.

Без оконного Фурье, в смысле - без обновляемой в реальном времени спектрограммы, - такая замечательная штука была бы невозможна!

[B@R5uk]

Cos(x-pi/2) в сообщении #1566682 писал(а):

В обеих частотах таким простым путём удалось найти четыре верных знака!

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

Cos(x-pi/2) в сообщении #1566682 писал(а):

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

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

Разумеется, если заранее известна модель сигнала, то будет разумнее воспользоваться ею для распознавания компонент. Например, если известно, что сигнал представляет собой сумму небольшого числа синусов, то можно просто зафитить сигнал с помощью МНК, подгоняя частоты численно с помощью, например, алгоритма Нелдера-Мида, взяв за отправную точку набор частот, рассчитанные через максимумы БПФ (в том числе предложенным вами способом). Немного трудоёмко, за то абсолютно корректно.

[Евгений Машеров]

Это решение задачи, но другой. Как получить оценку частоты пика с точностью выше, чем шаг частотной шкалы. Накопление по эпохам и windowing для другого, они нужны потому, что оценка спектра посредством DFT не является состоятельной статистически, то есть её дисперсия не падает с увеличением длины отрезка. И чтобы повысить точность, либо разбивают длинный сигнал на отрезки и полученные по ним оценки усредняют, дисперсия падает обратно числу усреднений, либо усредняют сам спектр, а это мероприятие в частотной области оказывается во временной области эквивалентно наложению окон. И то, и то повышает точность оценок мощности спектра ценой загрубления по частоте.
Предложенный же подход основан на предположении, что истинный спектр вблизи пика можно аппроксимировать параболой, а положение её максимума оценить, подогнав по трём точкам. Если пиков много - эффект может быть и обратен желаемому.
А если мало - можно подогнать авторегрессию порядка $2n$, где n - известное число синусоидальных составляющих, от параметров авторегрессии к их частотам перейти легко. Ну, или подгонять синусоиду. Только у нас вместо Нелдера-Мида Левенберг-Марквардт был, там линеаризуется хорошо.

[B@R5uk]

Евгений Машеров в сообщении #1566690 писал(а):

Только у нас вместо Нелдера-Мида Левенберг-Марквардт был, там линеаризуется хорошо.

Амплитуды синусов/косинусов тоже входили в вектор-параметр оптимизации вместе с частотами?

[ilghiz]

Подпишусь под комментарием Евгения Машерова по поводу минимизации методом Левенберга-Марквардта и от себя добавлю (выше уже писал) что и Ньютон и квази-Ньютон тут хорошо работает.

B@R5uk в сообщении #1566692 писал(а):

Амплитуды синусов/косинусов тоже входили в вектор-параметр оптимизации вместе с частотами?

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

[B@R5uk]

ilghiz в сообщении #1566693 писал(а):

фактически одномерная минимизация

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

$$\omega=2\arcsin\left(\frac{1}{2}\sqrt{-\frac{D^2Y}{Y}}\right)$$

Наличие шума, правда, всё портит, так как вектор производной содержит исходный несмещённый вектор:

$$D^2Y=Y_{-1}+Y_{+1}-2Y$$

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

$$Y\to\frac{1}{2}(Y_{-2}+Y_{+2})$$

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

[ilghiz]

B@R5uk в сообщении #1566699 писал(а):

Наличие шума, правда, всё портит

вы сами и сказали, тут даже и комментировать нечего, поэтому и фиттят наименьшими квадратами.

[B@R5uk]

ilghiz в сообщении #1566700 писал(а):

вы сами и сказали

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

[B@R5uk]

B@R5uk в сообщении #1566699 писал(а):

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

Как-то так:

$$\frac{D^2Y}{Y_{-2}+Y_{+2}}=-\frac{1}{k}$$

$$\omega=\arccos\left(\frac{2k+2}{k+\sqrt{k^2+8k+8}}\right)$$

Для надёжности ещё вектора в частном через 5 сэмплов проредить надо, тогда оценка для частоты (в пренебрежении нелинейностями формул) будет точно несмещённой.

[B@R5uk]

Вместо последней лучше вот такой формулой воспользоваться:

$$\omega=2\arcsin\sqrt{\frac{-k-2+\sqrt{k^2+8k+8}}{2k+2\sqrt{k^2+8k+8}}}$$

А то при малых частотах арккосинус близкой к единице величине не очень точно считается (особенно для какого-нибудь float'а).