Можно ли ускорить дискретное преобразование Фурье?

[TripleKill]

Всем привет.
Вопрос скорее из области теории.
Задача такова: есть сигнал от некого датчика, в котором полезная составляющая находится на определенной частоте, порядка десятков-сотни с небольшим кГц. Нужно отфильтровать эту частоту и вычислить ее амплитуду.
Применяю ДПФ, хочу ускорить расчет, но есть одна загвоздка: для быстрого преобразования требуется, чтобы количество точек замеров (N) сигнала выражалось степенью двойки.
Замеры осуществляются АЦП по таймеру, частота работы которого (Fs) со всеми предделителями и прочим не подгоняется под степень двойки. Т.е. сделав число замеров степенью двойки, я получаю шаги по частотам гармоник df=Fs/N, которые целое число раз не укладываются в фильтруемую мной частоту. Т.е. вместо одного пика на нужной частоте я получу два пика на соседних гармониках, из которых не особо-то и вычислишь искомую амплитуду.
Собственно, вопрос. Могу ли я в такой ситуации ускорить ДПФ какими-нибудь другими методами, не подгоняя специально количество замеров или частоту их взятия, но используя тот факт, что искомая частота находится на точно известной гармонике? Или, может, тут вообще лучше применить другую математику (не ДПФ, или его частный случай)? Спектр рассчитывать не требуется, только вытащить амплитуду одной заданной частоты.

[BlackKilkennyCat]

может, раз многое известно, вообще не вычислять, а сравнивать с шаблонами?

[jcxz]

TripleKill писал(а):

Собственно, вопрос. Могу ли я в такой ситуации ускорить ДПФ какими-нибудь другими методами, не подгоняя специально количество замеров или частоту их взятия, но используя тот факт, что искомая частота находится на точно известной гармонике? Или, может, тут вообще лучше применить другую математику (не ДПФ, или его частный случай)? Спектр рассчитывать не требуется, только вытащить амплитуду одной заданной частоты.

Если спектр не нужен - зачем тогда ДПФ? Сделайте полосовой фильтр и измеряйте амплитуду на его выходе. Фильтр требует меньше вычислений чем ДПФ/БПФ.

[КРАМ]

TripleKill писал(а):

Всем привет.
Вопрос скорее из области теории.

Во-первых, количество фильтров в ДПФ может быть произвольным. Можно посчитать только один, например.
Во-вторых, каждый фильтр ДПФ является квадратурным приемником прямого преобразования с полосой пропускания определяемой длиной накопления сигнала. То есть вы ничего не сказали о полосе сигнала. Лучше всего, если вы расскажите что это за сигнал (не его параметры, а откуда вы его взяли для измерений).
В третьих, если сигнал попал в два фильтра, то можно посчитать его амплитуду как корень квадратный из суммы квадратов (энергий) сигнала в каждом фильтре.
Итого, фильтры ДПФ должны перекрывать весь диапазон изменения, либо вы должны считать только один или два по известной частоте сигнала.

jcxz писал(а):

Сделайте полосовой фильтр и измеряйте амплитуду на его выходе.

А что это как не ДПФ? Это и есть он самый. Только его один бин.

[TripleKill]

Во-первых, количество фильтров в ДПФ может быть произвольным. Можно посчитать только один, например.

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

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

считать только один или два по известной частоте сигнала

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


Здесь Z0 и Z1 - действительная и мнимая части (если ничего не путаю), Sig[] - массив с замерами сигнала, u - номер гармоники, соответствующий искомой частоте, Npoints - количество замеров сигнала (= длина накопления).
Компилятор, конечно, использует возможности контроллера в части вычисления квадратного корня, но все равно достаточно долго. Один из вариантов, не относящихся чисто к математике - это затабулировать синусы и косинусы, что избавит контроллер от многократного их вычисления. Но есть минус - тогда я не смогу вычислить амплитуду произвольной частоты.

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

Сигнал взят с датчика. Про сам этот датчик мне немногое известно, но форма сигнала такая:

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

[Eddy_Em]

Ну вот если бы сразу поставили вопрос правильно, вам бы сказали, что Фурье здесь не нужно! Вас не волнует фазовая составляющая, поэтому можно было бы использовать, например, дискретное косинусное преобразование. Если частота известна достаточно точно, то вообще сверткой можно ограничиться, если же известен диапазон — тоже можно найти оптимальнй метод.
А преобразования Фурье пихать куда попало не следует!

[КРАМ]

Eddy_Em писал(а):

Вас не волнует фазовая составляющая, поэтому можно было бы использовать, например, дискретное косинусное преобразование.

Даладна!!!
Амплитуда и фаза требуют одинаковых расчетов , поскольку им обеим требуется результат в комплексной алгебраической форме. Получение одной проекции сигнала возможно лишь при условии КОГЕРЕНТНОСТИ преобразования. То есть когда фаза равна 0 или 180 градусам. Что само по себе будет технической проблемой.
Вычисление ДПФ и есть по сути интеграл свертки, есличо...

[tonyk]

Полоса пропускания фильтра должна быть 125 Гц, достаточно узко

125 Герц или кило Герц? Чё-т путаешься в показаниях. Или я чё-то не так понял.

При таких исходных данных получается, что тебе нужен просто узкополосный полосовой фильтр, а не БПФ. А если ещё и 125кГц, то вообще можно даже подумать и сделать его на ОУ. Если же 125Гц (не кГц!), то лучше в "цифре", тут аналоговый будет непростым.

[jcxz]

tonyk писал(а):

125 Герц или кило Герц? Чё-т путаешься в показаниях. Или я чё-то не так понял.

Видимо всё-таки - 2-ое.
125 кГц - это видимо центральная частота ППФ;
125 Гц - это видимо ширина полосы пропускания ППФ (определяется по максимально допустимому затуханию сигнала в полосе пропускания).
Никакой путаницы нет.

[КРАМ]

TripleKill писал(а):

Компилятор, конечно, использует возможности контроллера в части вычисления квадратного корня, но все равно достаточно долго. Один из вариантов, не относящихся чисто к математике - это затабулировать синусы и косинусы

Это "шуткатакой"?
Вы всерьез считаете синусы-косинусы аналитически? Милостивый государь, забудьте как дурной сон аналитические выражения ДПФ из учебника, патамушта при реализации в реальной программе нужно понимать откуда и как процессор извлекает данные и куда он их потом размещает. Нужно понимать размерность операндов и результата, чтобы потом не искать ресурсы на простейшие вычисления.
Итак.
1. Расчет фильтров ДПФ лучше (и нужно) осуществлять в формате fixed-point.
2. Коэффициенты ДПФ (те самые синусы и косинусы) должны быть константами во флеше.
3. Извлечение корня придется написать самому, иначе компилятор потребует флоаты. Тем более, что извлечение целочисленного корня совсем не ресурсоемкая функция. Могу предложить метод последовательных приближений, когда число итераций будет равно разрядности результата. Метод дает максимально точный результат ограниченный только заданной разрядностью.
4. К сожалению, АРМы не имеют trueDSP инструкций, поэтому адресация операндов при накоплении произведений потребует отдельных инструкций цикла.
Полоса 125 Гц потребует накопления массива сигнала в течении 8 мс. То есть с учетом частоты Найквиста и возможностей реализации антиалиасинга это будет не менее 4096 отсчетов.

Добавлено after 14 minutes 45 seconds:

tonyk писал(а):

а не БПФ

Автор темы ни словом не заикнулся про БПФ. Разговор шел о ДПФ.

Добавлено after 19 minutes 27 seconds:

TripleKill писал(а):

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

А вот тут есть некоторая проблема.
Дело в том, что фильтр ДПФ можно вычислять только один, если частота дискретизации сигнала кратна центральной частоте фильтра. Это значит, что все будет совсем просто, если частота сигнала будет известна ДО ТОГО, как будет запущено накопление массива. Так можно запустить преобразование на х4 или х8 частоте сигнала (фильтра) и вся "таблица" синусов-косинусов станет предельно примитивной. Для х4 эта таблица будет иметь всего 4 значения для каждой функции. Например в формате int16_t - cos = {0x7FFF, 0, 0x8000, 0}. Это и есть fixed-point: 0x7FFF = 0,9999695 0x8000 = -1

Добавлено after 16 minutes 33 seconds:

TripleKill писал(а):

В текущем варианте она так и вычисляется.

Только что заметил. Корень квадратный из суммы квадратов Real и Imaginary - это У ВАС вычисление амплитуды В ОДНОМ фильтре. Чтобы найти амплитуду сигнала попавшего между фильтрами, нужно сложить квадраты амплитуд ДВУХ фильтров и только потом извлечь корень.

[jcxz]

КРАМ писал(а):

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

Ну - если у автора ядро поддерживает аппаратные флоаты, то быстрее будет перевести fixed->float и воспользоваться соответствующей инструкцией CPU. Например VSQRT у Cortex-M4F выполняется за какие-то считанные такты (у меня отладчик насчитал 4 такта, может и врёт но около того). Если по таймеру - 6 тактов. Это на аргументе 2^31-1.
Целочисленное же вычисление корня от того же аргумента на этом же самом CPU методом Ньютона-Рафсона при полной оптимизации == ~89 тактов (и отладчиком и по таймеру).
Так что - глупо не воспользоваться возможностями FPU, коли оно и так имеется. Другое дело если его нет....

КРАМ писал(а):

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

А сколько это будет в тактах (для аргумента = 2^31-1)? Например для Cortex-M.

[КРАМ]

jcxz писал(а):

если у автора ядро поддерживает аппаратные флоаты

Если поддерживает, то можно. Но на самом деле это самая непринципиальная операция. Она ОДНА на весь фильтр. Накопление 2*4096 произведений с адресациями даст более 32000 машинных циклов В ЛУЧШЕМ случае...

jcxz писал(а):

А сколько это будет в тактах (для аргумента = 2^31-1)? Например для Cortex-M.

Не могу сказать. Не считал. Могу дать ассемблерный код для dsPIC33 с комментариями и Вы легко переведете его в любую систему команд.

[jcxz]

КРАМ писал(а):

Не могу сказать. Не считал. Могу дать ассемблерный код для dsPIC33 с комментариями и Вы легко переведете его в любую систему команд.

Лучше - на си. Так как я не знаком с системой команд PIC.
Да и подозреваю, что там будет всё тот же Ньютон-Рафсон.

[КРАМ]

jcxz писал(а):

Так как я не знаком с системой команд PIC.

Система команд PIC24/dsPIC33 читабельна настолько, что даже комментарии будут лишними, но я их приведу. Готового варианта на Си у меня нет, а писать лениво, извините.

Код:

;--- функция извлечения квадратного корня
;--- W0=SQRT[W1:W0]
SQRT:
   mov      # 0x8000, W5       ; 1 в старший разряд W5
   clr      W6                       ; W6=0 сбрасываем регистр результата
LoopSQRT:
   ior      W6, W5, W6          ; W6=W6 | W5
   mul.uu      W6, W6, W2  ; W2:W3=W6*W6 возводим в квадрат
   sub      W0, W2, W4           ; W4=W0:W1-W2:W3 вычитаем 32-разрядные операнды
   subb      W1, W3, W4            ; тоже самое
   bra      C, $+4                    ; пропускаем одну инструкцию, если есть перенос
   sub      W6, W5, W6            ; W6=W6-W5
   lsr      W5, W5                   ; правый логический сдвиг W5
   bra      NC, LoopSQRT          ;переход на LoopSQRT если нет переноса
   mov      W6, W0                    ; W0=W6 возвращаем результат в W0 (соглашение Си о возврате значения)
   return
;----------


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

[TripleKill]

Если поддерживает, то можно

Там и больше можно) Прошу прощения, забыл указать, что это все считается на STM32F407, а там этот самый FPU есть.

Мне доводилось писать аналог такого фильтра на arm-asm, и с использованием FPU все это весьма лаконично выглядит, корень одной строкой. Плюс MAC-инструкции.

2. Коэффициенты ДПФ (те самые синусы и косинусы) должны быть константами во флеше.

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

фильтр ДПФ можно вычислять только один, если частота дискретизации сигнала кратна центральной частоте фильтра

А вот это очень интересно. Я и решил написать первый пост, думая о чем-то таком: например, если частоты кратны, то можно брать не все отсчеты, а каждый n-ый, или, как Вы говорите, таблица синусов/косинусов упростится.

Чтобы найти амплитуду сигнала попавшего между фильтрами

Полагаю, это верно только в том случае, если искомая частота попадает точно посередине двух гармоник. А если она смещена к одной из них?

вообще можно даже подумать и сделать его на ОУ

Реализация на ОУ уже существует. Она была сделана до меня и не устраивает, скажем так, "заказчика", по точности. И нужно в любом случае получить оцифрованное значение амплитуды сигнала.

Автор темы ни словом не заикнулся про БПФ

Я не упоминал сам термин БПФ, но в первом посте написал, что читал про методы ускорения расчета. БПФ в плане дискретности, если не ошибаюсь, - это как раз про степени двойки.

[КРАМ]

TripleKill писал(а):

Это очень удобно, если частота, которую требуется отфильтровать, не будет меняться. Но если вдруг потребуется перестроить фильтр на другую частоту, то будет проблема.

Это чушь, извините. Таблица синуса/косинуса будет ОДНА. Частоты бинов ДПФ определяются частотой дискретизации. Выбор коэффициентов для фильтров - это вопрос расчета шага адресации таблицы. И все.

TripleKill писал(а):

Я не упоминал сам термин БПФ, но в первом посте написал, что читал про методы ускорения расчета. БПФ в плане дискретности, если не ошибаюсь, - это как раз про степени двойки.

Причем тут степень двойки? БПФ использует симметрию коэффициентов. Но расчет БПФ возможен только для ВСЕХ фильтров. Поэтому существует точка равенства времени расчета для БПФ и ДПФ и она определяется потребным числом фильтров из полного их количества. Кроме того, расчет элементарной операции БПФ - "бабочки" достаточно объемен и даже при наличии trueDSP занимает более, чем пару десятков инструкций. А таких "бабочек" будет N*logN. Для ДПФ расчет одного фильтра - это 2*N MAC-инструкций (при условии trueDSP, чего у ARM-ов нет, что приводит к фактической потере производительности примерно в 4 раза)
ЗЫ. trueDSP MAC - это не просто умножение-сложение за 1 цикл. Это еще и предвыборка одновременно двух операндов и перевод указателей для следующей предвыборки. И все за 1 машинный цикл. Этого у ARM-ов нет, поскольку требует иной архитектуры шин.

Добавлено after 10 minutes 32 seconds:

TripleKill писал(а):

Полагаю, это верно только в том случае, если искомая частота попадает точно посередине двух гармоник. А если она смещена к одной из них?

Полагаете неправильно. Сумма квадратов амплитуд - это сумма энергий. Корень квадратный из энергии - это амплитуда. Соотношение энергий не имеет значения.

Добавлено after 3 minutes 38 seconds:

TripleKill писал(а):

А вот это очень интересно. Я и решил написать первый пост, думая о чем-то таком: например, если частоты кратны, то можно брать не все отсчеты, а каждый n-ый, или, как Вы говорите, таблица синусов/косинусов упростится.

К кратности частот имеет отношение попадание сигнала точно в один фильтр. Если вы будете одновременно измерять частоту сигнала и формировать частоту дискретизации кратно 4 или 8, то вы фактически сделаете следящий приемник прямого преобразования с комплексным (квадратурным) выходом. И останется только найти модуль по упомянутой формуле Эйлера.

[Andrey_B]

TripleKill писал(а):

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

Гляньте алгоритм Гёрцеля https://ru.wikipedia.org/wiki/%D0%90%D0 ... 0%BB%D1%8F . В свое время применяли для АОНизма аж на 8080. В сети есть примеры на Си, например это http://we.easyelectronics.ru/samnetmon/ ... celya.html.

[TripleKill]

Таблица синуса/косинуса будет ОДНА

Для ясности приведу скрин из программы, в которой это моделировалось:

Я имел в виду вот этот расчет. При вычислении p используется номер гармоники u, который зависит от искомой частоты. Если при изменении искомой частоты я изменю частоту выборки таким образом, что u не изменится, то и p будут теми же (= одна и та же таблица синусов/косинусов). Но я не всегда могу обеспечить такую подгонку частоты, стало быть, в некоторых случаях p все же изменится.
Ну а если искомая частота неизменна, то, конечно, нет никакого смысла их пересчитывать.

Полагаете неправильно

Поясню предположение. Допустим, есть спектр частотного анализа:

Желтым изображены величины гармоник, полученные при расчете. Нарисуем огибающую (черную).
И, условно, есть некая искомая частота, не попавшая точно в отсчет. Когда я вычисляю корень из суммы квадратов амплитуд, амплитуду какой гармоники я получу? А или В? Математически она, конечно, больше похожа на А, но что, если искомая частота находится на месте В?

Добавлено after 18 minutes 34 seconds:

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

Да, спасибо, тоже на него в поиске натыкался. Интересно бы их точность сравнить.

[КРАМ]

TripleKill писал(а):

Когда я вычисляю корень из суммы квадратов амплитуд, амплитуду какой гармоники я получу? А или В? Математически она, конечно, больше похожа на А, но что, если искомая частота находится на месте В?

Наличие откликов в двух (и более) соседних фильтрах ничего не говорит о форме реального сигнала попавшего в эти фильтры. То есть разговоры об АМПЛИТУДЕ имеют смысл только если мы знаем, что это ОДИН ГАРМОНИЧЕСКИЙ сигнал (или одна гармоника более сложного сигнала). В таком случае нет никаких А и В гармоник. Есть один сигнал, энергия которого распределилась между двумя фильтрами в соотношении расстроек частоты сигнала относительно центральной частоты фильтров. Поэтому, сложив энергии ЭТОГО СИГНАЛА, попавшего в оба фильтра мы получим энергию этого сигнала и зная, что это синусоида, легко восстановим амплитуду путем извлечения корня. Если сигнал несинусоидальный и его форма может быть любой, то спектральный расчет может быть произведен только относительно ДЕЙСТВУЮЩЕГО значения сигнала. Разговоры об амплитуде перестают иметь смысл.
Еще раз повторю. ДПФ не выделяет гармоники сигнала. ДПФ - это набор фильтров с определенной полосой пропускания. Что попадает в фильтр, то мы и увидим. Аналитическое Фурье-преобразование имеет дело с БЕСКОНЕЧНЫМИ ПЕРИОДИЧЕСКИМИ сигналами. Поэтому в аналитической фурье-спектрограмме гармоника отображается дельта-функцией (бесконечно узкой палкой).
Вы привели сигнал, который подлежит анализу. В этом сигнале НЕ МОЖЕТ быть гармоник не кратных частоте полного сигнала. То есть при полосе анализа В ТЫСЯЧУ РАЗ МЕНЬШЕ частоты анализируемой гармоники, проблема может возникнуть только если эта гармоника является какой нибудь пятисотой или около того... Но в таком случае нужно будет увеличить время накопления...

[TripleKill]

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