Скользящие ДПФ

[Алексей Климов]

Здравствуйте.

Требуется устройство контролирующие параметры сигнала (частоту, фазу, амплитуду).

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

Использую скользящие ДПФ (Модифицированный алгоритм Герцеля).
Формат данных - single(float).

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

В статье "Accurate, Guaranteed - Stable, Sliding DFT", Krzysztof Duda, приводится метод стабилизации по следующей формуле

21-09.PNG (4.68 КБ) 8569 просмотров

, где 0 << r < 1 (пробую r = 0.99999). Этот метод приносит в жертву точность (как написано в статье).

Вопрос в том, можно ли после преобразование домножить вычисленные отсчеты на какой-нибудь коэффициент > 1 для повышения точности и как его рассчитать? Пробовал подбирать коэффициенты для каждого отсчета, но это портит интерполяцию ДПФ, возможно нужно подбирать один коэффициент на все отсчеты.

И главное, нет ли более правильного пути решения проблемы?

[Бахурин Сергей]

Добрый день. Необходимо понять причину неустойчивости. Это проблема округления комплексного фазового коэффициента? Если да, то можно попробовать сместить полюс фильтра в зону устойчивости.
Было бы здорово, если вы покажите проблему более детально. В идеале - дадите данные на которых алгоритм расходится.

[Алексей Климов]

Насколько я понял статью, то да - ошибка округления фазовых коэффициентов приводит к неустойчивости и умножение их на коэффициент r < 1 (в статье используется r = 0.99999) перемещает их внутрь единичной окружности.

В архиве две картинки: dft_res - выход реальной части ДПФ (в виде ких фильтра) вблизи максимума, sdft_res - тоже для скользящего ДПФ и сигнал поданный на вход (фазовая модуляция 5 Гц на несущей 50 Гц, сохранил из матлаба переменную).

Фильтр для sdft:

Параметры фильтра:
Длина окна M = 1152, номер отсчета ДПФ k = 3, частота дискретизации 19200 Гц.
Соответственно коэффициенты фильтра:
a = cos(2*pi*3/1152)
b = sin(2*pi*3/1152)

После изменения фильтра по формуле из первого сообщения результат следующий:

Похоже, что результат стабилизируется, но погрешность существенна.

[Бахурин Сергей]

Вопрос такой. Насколько важно корректность фазы? Или амплитуды будет достаточно?

Ещё вопрос в схеме задержка , чему равно d?

Попробовал ваши данные. Если все в формате double, то все устойчиво. Неустойчивость возможна если сильно округлить поворотный коэффициент (до 4 знаков после запятой в большую сторону). Но замечу что задержка d =M+1, а не M (часто сам ошибаюсь).

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

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

[Алексей Климов]

Модель должна мерить фазу входного сигнала с точностью не менее 0.1 градус.

Задержка 1152 равна длине окна M. Вот схема из статьи:

Ошибка? Попробовал 1152+1 и результат совсем странный.

Данные в double, но на вход подается single. Я их так и преобразую single(x[k]) как и поворотные коэффициенты и после этого 3 отсчет становится неустойчивым.

Стабилизирую я вот такой схемой но точность падает:

rsdft.PNG (8.66 КБ) 8505 просмотров

Вы о ней говорите?

[Бахурин Сергей]

Данные могут быть любые. Главное чтобы коэффициент w был без округления.

Задержка правильная на длину окна 1152

Схема стабилизации иная. Умножитель стоит на входе и смещает бин в ноль. Потом стабильный фильтр с единичными к-тами.

Требуется оценить 0.1 градус по отношению к чему?

[Алексей Климов]

Я пробовал нижнюю модель. Для нее неустойчив 1 бин.

Первая ведь эквивалентна?

По отношению к фазе входного сигнала. Параметры входного гармоники нужно измерить в итоге.

[Бахурин Сергей]

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

Если у вас неустойчивость, то надо понять причину.

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

[Алексей Климов]

Вот так у меня реализовано MSDFT:

Counter считает от 0 до 1151 с частотой 19200 Гц.

Задержка d = 1152.

В таблицах записаны коэффициенты (k - номер бина):
m = 0 : 1151;
M = 1152;
ReW = cos(-2*pi * k * m ./ M);
ImW = sin(-2*pi * k * m ./ M);

ReWshift и ImWshift сдвинутые влево циклически ReW и -1*ImW.
Т.е. W^(-km) и W^(k(m+1)).

Комплексное умножение:

Можете оценить правильность?

Про фазу. На текущем примере у нас несущая на 50 Гц и модулированная фазы в 5 Гц . Вот такой же сигнал нужно получить на выходе устройства. Я это делаю через arctg отношения мнимой и реальной части(после прореживания в Fs/50 раз), а при анализе и сравнение применяю unwrap.
Вот эта картинка убедила меня, что я понимаю что происходит:

phase.PNG (27.85 КБ) 8468 просмотров

И, кажется, ответ на ваш вопрос - относительно гармоники 50 Гц.

[Алексей Климов]

О алгоритме MSDFT в статье авторы делают следующие замечание:

Использование формата single (плавающей точки), видимо, причина неустойчивости.
Но я не понял идею с фиксированной точкой...