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

КРАМ

[uquote="TripleKill",url="/forum/viewtopic.php?p=3900708#p3900708"]Можете прокомментировать...?[/uquote]
Чтобы увести разговор из формализованной математики, переведу стрелки на физический смысл.
Итак.
Каждый бин ДПФ - это квадратурный приемник прямого преобразования. И гетеродин (гетеродины) у него синусно/косинусная последовательность кратной первому фильтру частоты. Значит в полный массив сигнала всегда укладывается полное число периодов гетеродинов ЛЮБОГО из фильтров. Сиречь, имея ОДНУ таблицу синуса или косинуса с дискретностью по фазе равной числу отсчетов входного массива, мы сможем брать готовые отсчеты синуса и косинуса, просто рассчитывая адресное смещение и адресный шаг. Ну или создав для каждого фильтра две константы - смещение и шаг адреса.
Естественно, что если фильтр один, то полная таблица превращается в прореженный огрызок ОДНОГО периода гетеродина, который на полном массиве просто повторяется некоторое количество раз.

Реклама

jcxz

[uquote="КРАМ",url="/forum/viewtopic.php?p=3900471#p3900471"]Все достаточно просто. Поразрядно, начиная со старшего разряда добавляем 1 и возводим в квадрат. После чего сравниваем с исходным аргументом. Оставляем единицу, если меньше или равно.[/uquote]Спасибо! Жаль не могу почему-то поставить плюс.

Перевёл на си (если не ошибся в чём-то):

Код: Выделить всё

u32 Sqrt(u32 x)
{
  u32 w5 = B15, w6 = 0;
  do {
    w6 += w5;
    if (w6 * w6 > x) w6 -= w5;
  } while (w5 >>= 1);
  return w6;
}

скомпилил с макс. оптимизацией по скорости:

Код: Выделить всё

         MOV      R6,#+32768
         MOVS     R2,#+0
DteSqrt: ADDS     R2,R6,R2
         MUL      R7,R2,R2
         CMP      R1,R7
         IT       CC
         SUBCC    R2,R2,R6
         LSRS     R6,R6,#+1
         BNE      DteSqrt

получилось ооооочень медленно.

Результат для sqrt(127):
Ньютон_Рафсон = 63 тактов
посл_приближение = 192 тактов
Результат для sqrt(2^31-1):
Ньютон_Рафсон = 89 тактов
посл_приближение = 192 тактов
Выполнение из ОЗУ.
Видимо - сказывается большое число итераций в вашем методе и высокая цена за переход: (1+1+1+1+1+1+5)*16+k=192 такта.
К тому же у НР с уменьшением величины числа уменьшается время вычисления.

PS: Так что для ARM-ов последовательное приближение - не самый лучший способ.

Реклама

КРАМ

Я писал для 16-битных Микрочипов, причем без применения инструкции do, чтобы переваривалось PIC24, в котором этой инструкции нет. Плюс к этому конвейер у Микрочипа короче и переход выполняется за 2 или 3 маш. цикла.
Итого в предложенном мной варианте в PIC-ах получается 9*16+2=146.
Опять же нужно посмотреть результат компиляции по Н-Р. Он может плохо лечь на систему команд Микрочипа. Переменная длительность исполнения ничего не дает по факту, как вы сами понимаете, важен самый длинный вариант.
Ну и слишком малый вес корня в общей производительности Фурье не слишком мотивирует к поискам оптимума по скорости.

12val12

Не нужен тут ДПФ
НУЖЕН просто фильтр настроенный на одну частоту
FIR или IIR

Добавлено after 3 minutes 49 seconds:
автор какая (f/FS ) ?
работы на 3 мин
http://www.winfilter.20m.com/

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

Реклама

КРАМ

[uquote="12val12",url="/forum/viewtopic.php?p=3910335#p3910335"]Не нужен тут ДПФ[/uquote]
Вообще то автору не нужен сигнал на выходе фильтра. Ему нужна АМПЛИТУДА.
И что ему делать с полоснопропускающим фильтром для этого?

Реклама

12val12

Вычислить амплитуду через 3 тау установления фильтра
Не есть проблема .

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

Реклама

КРАМ

[uquote="12val12",url="/forum/viewtopic.php?p=3910656#p3910656"]Вычислить амплитуду через 3 тау установления фильтра[/uquote]
Я не спрашивал КОГДА, я спрашивал КАК (впрочем, "когда" является отдельной проблемой - от какого момента считать время?).
На выходе бина ДПФ мы получим СРАЗУ комплексное значение сигнала. Из которого элементарно находится амплитуда.
А два бина позволяют легко решить проблему дробной частоты.
На выходе фильтра мы получим просто мгновенное значение сигнала. И что с ним делать? Особенно если нет точного значения частоты....

12val12

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

12 бит Ацп центруется делителем R=R ref+ref- на коде 2047 (- + 50 не критично ! )
в начале программы дополнительно вычисляем усреднением 0
ловится перепад через ноль на фронте ->+
находится максимальное значение
находится минимальное значение
A= (max - min)/2
пока не следующий перепад через ->+
амплитуды соседних колебаний складыааются последовательно в массив

КРАМ

[uquote="12val12",url="/forum/viewtopic.php?p=3911304#p3911304"]ловится перепад через ноль на фронте[/uquote]

Помилосердствуйте, отец родной!!! Какой перепад через ноль? Ничего, что при дискретизации 4...5 семплов на период там никаких максимумов ВООБЩЕ НЕТ? Точнее попадание мгновенного значения на максимум можно ждать очень долго, поскольку фаза дискретизации относительно измеряемого сигнала может оказаться любой и асинхронной.
Извините, но все это эпичный бред. В узкополосном фильтре сигнал практически синусоидальный. Его амплитуда выясняется в полтычка ДВУМЯ квадратурными отсчетами. И нах не облокотилось чего то там ловить...
ЗЫ. К чему вы написали о нуле АЦП я так и не понял... Мы обсуждали как беззнаковый формат превратить в знаковый? Чего то я не припомню...

12val12

Согласен с 1 абзацем раз всего 500 килосамплов .ну пусть так .посему гибко меняю алгоритм

Если частота дискретизации всего в 4 раза выше частоты сигнала то
суммируются все абсолютные значения октлонения отсчетов за 8-32 периодав и делятся на количество самплов за это время

автор молчит .может уже неактуально
Но если появится я ему предложу решение на фильтре .

КРАМ

[uquote="12val12",url="/forum/viewtopic.php?p=3911612#p3911612"]суммируются все абсолютные значения октлонения отсчетов за 8-32 периодав[/uquote]
Может стоит наконец прекратить рандомные фантазии и понять, что каждый бин ДПФ (кстати, с окном, есличо) и есть искомый FIR, причем сразу с квадратурным (комплексным) результатом. И вычисление амплитуды по этому результату займет около 100 машинных циклов. И все...
А пара бинов даст амплитуду для дробной частоты находящейся между этими бинами.
Для снятия экспоненциальной огибающей затухающих колебаний вообще не требуется никаких нелинейных мероприятий типа поисков-сортировок-отбора. Оные мероприятия являются самими ресурсоемкими действиями на ЛЮБОЙ аппаратной платформе. Вместо них в ДПФ используются MAC-инструкции.

КРАМ

[uquote="jcxz",url="/forum/viewtopic.php?p=3900892#p3900892"]PS: Так что для ARM-ов последовательное приближение - не самый лучший способ.

[/uquote]
Вернулся к вопросу.
Основная проблема метода Ньютона для dsPIC33/PIC24 - деление 32/16 (19 машинных циклов), а не 32/32. Из-за этого возникает ограничение входного аргумента величиной 32 762^2=1 073 348 644. Для сигнальных применений в 12 разрядов (внутренний АЦП) непринципиально, но нужно иметь ввиду.
Написал на АСМе метод Ньютона для dsPIC33 и получил время извлечения корня 117 машинных циклов. Для dsPIC33С - 65 машинных циклов, поскольку деление у них не 19, а 7 циклов.
Последовательное приближение дает полный диапазон (0...65 535^2=4 294 836 225) за 167 циклов и диапазон до 32 762^2=1 073 348 644 за 157 циклов. То есть у Ньютона выигрыш в 40 циклов.
Резюме. Основная экономия метода Ньютона возникает на длительности деления и расчете первого приближения путем нахождения старшего значащего разряда (логарифмирование по основанию 2). Для чего, например, в системе команд dsPIC33/PIC24 есть инструкция FF1L.

bad2cat

[uquote="TripleKill",url="/forum/viewtopic.php?p=3900336#p3900336"]Всем привет.
Вопрос скорее из области теории.
Задача такова: есть сигнал от некого датчика, в котором полезная составляющая находится на определенной частоте, порядка десятков-сотни с небольшим кГц. Нужно отфильтровать эту частоту и вычислить ее амплитуду.
Применяю ДПФ, хочу ускорить расчет, но есть одна загвоздка: для быстрого преобразования требуется, чтобы количество точек замеров (N) сигнала выражалось степенью двойки. .[/uquote]попробуйте Герцеля
https://livepcwiki.ru/wiki/Goertzel_algorithm
также есть несколько алгоритмов бпф переменной длины (сложной или комплексной, Винограда, Адамара)
таакже можно забить нулями отчёты до длины кратной степени двойки.
у меня такая же проблема стояла, когда делал аон на 580ик80 в 1992. сами понимаете выч.мощность аховая ) но зато нужно вычислять всего 6 частот и ьыстро через каждые 200 Гц. я взял 64 отчёта и окно Хэмминга на подставке (это уже на z80), а друг делал на герцеля 40 отчётов, но затем всётаки перешёл на ких фильтры (по мере изучения в вузе ))
кстати, он также как и вашей задаче, определчл только наличие частот (я же считал амплитуду для большей чёткости ), поэтому у него входное ацп это просто вход компаратора, а умножение 1 разряда на команде xor.
и ничего, работало норм.

zelya

имхо куда-то не туда ушли - если нужно диагностировать появление одной частоты из выборки сигналов - адаптивный алгоритм герцеля.
https://ru.dsplib.org/content/goertzel/goertzel.html - далее есть вариант на последовательное дополнение на каждой новой выборке
- http://www.dsplib.ru/content/goertzelmo ... elmod.html (недавно искал такую задачу)...
на произвольной частоте выборок, выше как минимум в 2 раза от требуемой частоты, см. теорему котельникова, а на необходимые частоты выставляет коэффициенты и вуаля...