Автор: Александр Калачев
Дата публикации: 2009.12.07
Публикуется впервые
Вычисление преобразования Фурье на процессоре SEAforth40
Введение
Рассмотрим задачу вычисления дискретного преобразования Фурье на процессоре SEAforth40. Ограничимся вычислениями в формате фиксированной точкой - т.е. только целочисленные операции. Для вычисления преобразования Фурье требуется введение комплексной арифметики. В данном случае ограничимся только действительными числами при помощи дополнительного преобразования - преобразования Хартли, которое является чисто действительным. Имея вычисленные отсчёты преобразованиям Хартли можно получить как действительные и мнимые коэффициенты Фурье преобразования, так и спектр мощности и фазовый спектр.
Дискретное преобразование Хартли вычисляется следующим образом:
H(v)=1/N*(Summ((t=0;N-1) f(t)*cas(2*pi*t*v/N)));
где
cas(y)=cos(y)+sin(y).
Одним из замечательных свойств этого преобразования является то, что обратное преобразование вычисляется по аналогичной формуле, за исключением масштабирующего коэффициента 1/N.
С преобразованием Фурье оно связано следующим образом.
гармоники спектра мощности:
Фр(v)=(H(v)^2+H(N-v)^2)/2;
фазовый спектр:
Фф(v)=arctg((H(v)-H(N-v))/(H(v)+H(N-v))).
Таким образом, задача вычисления преобразования Фурье разбивается на два больших этапа - вычисление преобразование Хартли и уже на его основе преобразования Фурье.
Трудоемкость задачи может быть оценена в N^2 операций сложения и умножения на вычисление H(v), плюс 3N умножений и сложений на получение отсчётов преобразования Фурье. Прямое вычисление Фурье преобразования даст порядка (2N)^2 операций сложения и умножения. Для быстрого преобразования Фурье имеем 2NlogN умножений и 3NlogN сложений. Быстрое Преобразование Хартли имеет порядка 3/2NlogN сложений и NlogN умножений.
Разбиение задачи
Учитывая скромные ресурсы процессора по встроенной оперативной памяти, ограничимся преобразованиями с небольшим числом N. Для ускорения вычислений БПХ распределим вычисления коэффициентов преобразования на различные ядра процессора - пусть каждое ядро вычисляет один или несколько коэффициентов преобразования.
Рассмотрим случай N=16. Для него возможен вариант, при котором одно ядро процессора вычисляет один коэффициент (аналогичная ситуация возможна и при количестве отсчётов сигнала порядка 32-х). Т.о количество операций выполняемых ядром можно оценить как N+3 операций сложения и умножения. Вычисления тригонометрических функций целесообразно выполнить табличным методом - каждое ядро будет иметь свой фиксированный набор значений функции cas() и arctg().
Для удобства в вычислениях будут задействованы центральные 16 ядер. Предполагается, что отсчёты сигнала поступают в процессор посредством одного из периферийных ядер. В рассмотренном ниже примере источником сигнала является ядро №10 (можно предполагать, что оно получает сигнал с внешнего АЦП), ядро №20 принимает сигнал после подачи его на преобразователь и может передать его на дальнейшую обработку другим ядрам. Ядро, вычисляющее коэффициент преобразования ждёт отсчёт входного сигнала, копирует его себе, передает следующему ядру, вычисляет произведение с накоплением. Результаты преобразования остаются в ядрах.
Диаграммы потов данных
В итоге имеем следующую диаграмму потоков данных при вычислении преобразования Хартли.
Код:
[20]<-f(t)-[21]<-f(t)-..<-[28]
^
|
f(t)
|
[10]-f(t)->[11]-f(t)->..->[18]
11е ядро вычисляет H(0), 12е - H(1), ... 28e H(8) и т.д.
На втором этапе - вычислении Фурье спектра диаграмма потов данных для рассматриваемого N выглядит следующим образом.
Код:
[20] [21] [28]
| |
H(15) H(8)
H(0) H(7)
| |
[10] [11] [18]
Где «|» соответствует взаимному обмену коэффициентами.
В силу симметрии ДПФ относительно N/2 для вычисления спектра мощности достаточно произвести вычисление только на половине ядер.
Алгоритмы работы ядерАлгоритм работы ядер, вычисляющих коэффициенты преобразований, практически одинаков, и может быть представлен в виде последовательности нескольких шагов.
1. Начальная инициализация, промежуточная сумма s=0;
2. Цикл от 0 до 15 из следующих ниже шагов;
3. прием входного отсчёта сигнала - f(t);
4. копирование сигнала себе;
5. передача отсчёта сигнала соседнему ядру;
6. выборка из таблицы значения cas(vt);
7. вычисление произведения p=p*cas(vt);
8. вычисление промежуточной суммы s=s+p;
9. если отсчёт не 16й, возврат на начало цикла - п.3.
10. H(v)=s;
11. прием H(N-v);
12. вычисление Фр(v)=(H(v)^2+H(N-v)^2)/2;
Ниже приводится один из вариантов реализации описанного алгоритма.
Основная программаhart2.vf реализует этапы с 1 по 10й.
fourier2.fv этапы 11, 12.
файл hartley_test.vfКод:
v.VF +include" c7Gr01/romconfig.f" \ подключение функций ROM нужной версии процессора
include <path>\Fpmath.f \ для вычисления таблиц функции cas() подключаем библиотеку плавающей точки компилятора SwiftForth
16 VALUE num \ =N количество отсчетов
0 VALUE v \ номер отсчет Хартли спектра
include coef0.vf \ определяем слова, вычисляющие cas в формате с фиксированной точкой
10 {node \ ядро выдает 16 последовательных отсчётов сигнала
0 org
here =p
'r--- # b!
1 # !b
1 # !b
1 # !b
1 # !b
0 # !b
0 # !b
0 # !b
0 # !b
1 # !b
1 # !b
1 # !b
1 # !b
0 # !b
0 # !b
0 # !b
0 # !b
node}
( код для ядер практически одинаков, за исключением портов приема и передачи данных, определяем несколько переменных, позволяющих настроить порты в коде для каждого из ядер
)
0 VALUE in_
0 VALUE out_
0 TO v \ - номер отсчёта Х
11 {node 0 org \ установка начального адреса компиляции
\ задание таблицы коэффициентов cas(vt)
coef , , , , , , , , , , , , , , , ,
here *cy =p \ включение режима расширенной арифметики
'r--- TO in_ '--l- TO out_ include hart2.vf \ настройка переменных и переход к вычислению коэффициента H(v)
\ -- hh hl
'---u TO in_ '---u TO out_ include fourier2.vf \ вычисление отсчёта мощностного спектра Фурье - Фр(v)
\ -- fh fl
node}
1 TO v
12 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy =p \ включение режима расширенной арифметики
'--l- TO in_ 'r--- TO out_ include hart2.vf
'---u TO in_ '---u TO out_ include fourier2.vf
node}
2 TO v
13 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy =p \ включение режима расширенной арифметики
'r--- TO in_ '--l- TO out_ include hart2.vf
'---u TO in_ '---u TO out_ include fourier2.vf
node}
3 TO v
14 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy =p \ включение режима расширенной арифметики
'--l- TO in_ 'r--- TO out_ include hart2.vf
'---u TO in_ '---u TO out_ include fourier2.vf
node}
4 TO v
15 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy =p \ включение режима расширенной арифметики
'r--- TO in_ '--l- TO out_ include hart2.vf
'---u TO in_ '---u TO out_ include fourier2.vf
node}
5 TO v
16 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy =p \ включение режима расширенной арифметики
'--l- TO in_ 'r--- TO out_ include hart2.vf
'---u TO in_ '---u TO out_ include fourier2.vf
node}
6 TO v
17 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy =p \ включение режима расширенной арифметики
'r--- TO in_ '--l- TO out_ include hart2.vf
'---u TO in_ '---u TO out_ include fourier2.vf
node}
7 TO v
18 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy =p \ включение режима расширенной арифметики
'--l- TO in_ '---u TO out_ include hart2.vf
'---u TO in_ '---u TO out_ include fourier2.vf
node}
8 TO v
28 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy =p \ включение режима расширенной арифметики
'---u TO in_ '--l- TO out_ include hart2.vf
$3fff0 # ~u/mod \ масштабирование коэффициентов
'---u TO out_ out_ # b! dup !b \ передача результата ядрам, вычисляющим отсчеты Фурье спектра
node}
9 TO v
27 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy =p \ включение режима расширенной арифметики
'--l- TO in_ 'r--- TO out_ include hart2.vf
$3fff0 # ~u/mod \ масштабирование коэффициентов
'---u TO out_ out_ # b! dup !b
node}
10 TO v
26 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy =p
'r--- TO in_ '--l- TO out_ include hart2.vf
$3fff0 # ~u/mod \ масштабирование коэффициентов
'---u TO out_ out_ # b! dup !b
node}
11 TO v
25 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy \ включение режима расширенной арифметики
'--l- TO in_ 'r--- TO out_ include hart2.vf
$3fff0 # ~u/mod \ масштабирование коэффициентов
'---u TO out_ out_ # b! dup !b
node}
12 TO v
24 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy =p \ включение режима расширенной арифметики
'r--- TO in_ '--l- TO out_ include hart2.vf
$3fff0 # ~u/mod \ масштабирование коэффициентов
'---u TO out_ out_ # b! dup !b
node}
13 TO v
23 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy =p \ включение режима расширенной арифметики
'--l- TO in_ 'r--- TO out_ include hart2.vf
$3fff0 # ~u/mod \ масштабирование коэффициентов
'---u TO out_ out_ # b! dup !b
node}
14 TO v
22 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy =p \ включение режима расширенной арифметики
'r--- TO in_ '--l- TO out_ include hart2.vf
$3fff0 # ~u/mod \ масштабирование коэффициентов
'---u TO out_ out_ # b! dup !b
node}
15 TO v
21 {node 0 org
coef , , , , , , , , , , , , , , , ,
here *cy =p \ включение режима расширенной арифметики
'--l- TO in_ 'r--- TO out_ include hart2.vf
$3fff0 # ~u/mod \ масштабирование коэффициентов
'---u TO out_ out_ # b! dup !b
node}
reset \ сброс процессора
11 watch1 10 setstep \ настройка детального просмотра состояния 11-го ядра, задали шаг симулята в 10 «тактов»
sim \ запуск симуляции
Вычисление коэффициентов ХартлиРассмотрим вычисление отсчетов Хартли
Ядро вычисляет сумму произведений отсчетов сигнала на функцию cas(vt):
H(v)=Summ((t=0;N-1) f(t)*cas(2*pi*t*v/N));
При необходимости масштабирование коэффициента выполняется отдельно.
Выполняются следующие действия:
1. Начальная инициализация, промежуточная сумма s=0;
2. Цикл от 0 до 15 из следующих ниже шагов;
3. прием входного отсчёта сигнала - f(t);
4. копирование сигнала себе;
5. передача отсчёта сигнала соседнему ядру;
6. выборка из таблицы значения cas(vt);
7. вычисление произведения p=p*cas(vt);
8. вычисление промежуточной суммы s=s+p;
9. если отсчёт не 16й, возврат на начало цикла - п.3.
10. H(v)=s.
Накопление суммы ведется с 36-ти битной точностью, входной сигнал и коэффициенты cas – 18-раздные.
Файл hart2.vfКод:
0 # 0 # 0 #
\ на стеке начальная сумма в формате двойной точности и номер отсчета входного сигнала -- sh sl t
15 # for \ начинаем цикл вычисления H(v)
\ -- sh sl t
dup dup xor dup . + drop \ очищаем бит переноса
in_ # b! @b \ -- sh sl t f(t) принимаем входной отсчёт
dup out_ # b! !b \ -- sh sl t f(t) скопировали на стек и передали следующему ядру
\ выбираем из таблицы коэффициентов cas(vt)
push \ сохраняем отсчет сигнала -- sh sl i r: -- f(t)
a! @a+ \ -- sh sl cas; a=t+1; r: -- f(t)
pop a@ push \ временно сохраняем номер отсчета на стеке возвратов -- sh sl cas f(t) ; r: -- t+1
( умножаем f(t) на cas(vt)
учитываем особенности реализации слова * - умноженная на 2 старшая часть на стеке; - сохранение множимого во втором элементе стека; в регистре а – младшая часть с инвертированным старшим битом
)
* 2/ \ умножаем и приводим к нормальному виду старшую часть произведения
\ -- sh sl cas cas*f(t) ; r: -- t+1
push drop pop \ избавляемся от множителя на второй позиции стека
a@ $20000 # xor \ помещаем на стек младшую часть
\ -- sh sl yh yl; r: -- t+1; где y=cas*f(t)
\ производим сложение чисел двойной точности
\ складываем младшие слова
push a! \ -- sh sl; r: -- yl t+1; a=yh
pop \ -- sh sl yl; r: -- t+1; a=yh
. + \ -- sh sl+yl; r: -- t+1; a=yh
\ складываем старшие слова
push \ -- sh ; r: -- sl+yl t+1; a=yh
a@ | . + pop \ -- sh+yh+c sl+yl t+1; r: -- t+1; a=yh
pop \ возвращаем на стек сохраненный номер отсчета -- sh+yh+c sl+yl t+1; a=yh
next \ -- hvh hvl t+1
drop \ сбрасываем ненужный теперь номер отсчета t -- hh hl - отсчет Хартли спектра
Вычисление преобразования ФурьеВычисление коэффициентов Фурье спектра ведется на основе вычисленных ранее коэффициентов Хартли спектра. Последовательность действий следующая:
1. получаем H(v);
2. принимаем от соответствующего ядра отсчет H(N-v);
3. вычисляем коэффициент Фурье спектра мощности Фр(v)=(H(v)^2+H(N-v)^2)/2;
Файл fourier2.vfКод:
\ на стеке вычисленный в hart2.vf отсчет Хартли (старшая и младшая часть) -- hh hl
$3fff0 # ~u/mod \ масштабируем коэффициент – делим на N (N=16) -- r q ( -- r hv )
( получившийся остаток можно не учитывать, а поскольку переполнения стека в данном процессоре не бывает, можно считать ячейки стека ниже hv пустыми)
in_ # b! @b \ принимаем масштабированный отсчет h(N-v) -- hv h(N-v)
\ вычисляем H(N-v)^2
dup \ дублируем
* \ умножаем
\ приводим произведение к нормальному виду
2/ push drop pop \ -- hv f2h ; a=f2l
a@ $20000 # xor \ -- hv f2h f2l ; a=f2l
\ временно сохраняем его на стеке возвратов
push push \ -- hv ; r: -- f2h f2l;
\ вычисляем H(v)^2
dup \ дублируем
* \ умножаем
\ приводим произведение к нормальному виду
2/ push drop pop \
a@ $20000 # xor \ -- f1h f1l ; r: -- f2h f2l;
\ складываем два числа с двойной точностью, одно из них находится на стеке возвратов
pop a! pop \ взяли младшую часть слагаемого со стека возвратов -- f1h f1l f2l; a=f2h
. + push \ сложили, результат запомнили на стеке возвратов -- f1h ; a=f2h ; r: -- f1l+f2l
a@ . + pop \ сложили старшие части, младшую часть суммы вернули на стек -- Фрh Фрl
\ на стеке удвоенный отсчет Фурье спектра мощности -- Фрh Фрl
……..
\ далее может идти код, нацеленный на последующую обработку или передачу результатов
'-d-- # b! @b \ останов ядра – исключительно для просмотра результата в симуляторе
Вспомогательный файл coef0.vf содержит слова, помогающие автоматизировать заполнение таблиц функции cas(vt).
Файл coef0.vfКод:
: cas ( f: x -- ; -- cas ) \ вычисляет целочисленный вариант cas(x)
fdup fcos fswap fsin f+ $01000 s>f f* f>s
;
: 2pivt/N ( N v t -- ; f: -- 2pivt/N ) \ формирует аргумент для cas
s>f s>f f* s>f f/ pi 2.e f* f*
;
: coef ( -- cas[N] cas[N-1] … cas[0] )
\ оставляет на стеке N отсчетов в обратном порядке для заполнения таблицы cas(vt)
0 num ?do
num v i 2pivt/N cas
-1 +loop
;
РезультатыОперативная память ядер вычисляющих и H(v) и Фр(v) загружена на 92% из которых 27% за-нято таблицей cas(vt).
Дамп памяти одного из ядер приведен ниже.
Код:
RAM Node 12
addr data mnemonics/code
000 01000 --
001 014E8 ¦
002 016A1 ¦
003 014E8 ¦
004 01000 ¦
005 008A9 ¦
006 00000 ¦
007 3F757 ¦
008 3F000 } таблица коэффициентов cas(vt)
009 3EB18 ¦
00A 3E95F ¦
00B 3EB18 ¦
00C 3F000 ¦
00D 3F757 ¦
00E 00000 ¦
00F 008A9 --
010 05D17 @p+ @p+ @p+ @p+ ---------------------------------
011 00000
012 00000
013 00000
014 0000F
015 2E9B2 push . . .
016 24DE3 dup dup xor dup
017 2C1EF . + drop @p+
018 00175
019 29F97 b! @b dup @p+
01A 001D5
01B 29BBA b! !b push .
01C 2BC9A a! @a+ pop . вычисление H(v)
01D 228B2 a@ push . .
01E 134CA call CA *
01F 308EA 2/ push drop .
020 26E12 pop a@ @p+ .
021 20000
022 388AA xor push a! .
023 269F2 pop . + .
024 2EEB2 push a@ . .
025 2C19A . + pop .
026 27016 pop next 16
027 3BDB2 drop @p+ . . ---------------------------------
028 3FFF0
029 136AE call 2AE ~u/mod
02A 04B03 @p+ b! @b dup
02B 00145
02C 134CA call CA *
02D 308EA 2/ push drop .
02E 26E12 pop a@ @p+ .
02F 20000
030 388BB xor push push dup вычисление Фр(v)
031 134CA call CA *
032 308EA 2/ push drop .
033 26E12 pop a@ @p+ .
034 20000
035 38CAA xor pop a! .
036 269F2 pop . + .
037 2EEB0 push a@ . +
038 27DA2 pop @p+ b! . ---------------------------------
039 00115
03A 00000 @b and @b +
03B
03C
03D
03E
03F
Временные параметры вычислений следующие:
- вычисление коэффициента Хартли спектра ~ 2940 тактов;
- вычисление Фурье ~ 350 тактов;
- общее время выполнения ~ 3320 тактов.
С учетом того, что такт примерно соответствует 1.4нс, возможно выполнение примерно 215000 преобразований в секунду.
ОптимизацияПриведенный выше код для вычисления H(v) и Фр(v) не вполне соответствует форт стилю программирования и более напоминает код на ассемблере, что впрочем, для данного процессора не далеко от истины. Более оптимальным, с точки зрения структурного программирования и экономии места было бы определение слов M*, D+ с вызовом их в нужных местах.
Файл math.vfКод:
: M* ( x y -- ph pl )
* 2/ \ -- x y ;
push drop pop
a@ $20000 # xor \ -- ph pl
;
: D+ ( ah al bh bl -- sh sl )
push a! pop \ -- ah al bl; a=bh
. + push \ -- ah ; a=bh ; r: -- sl
a@ . + pop \ -- sh sl
;
Код для ядер изменится следующим образом.
К примеру код для 11-го ядра
Код:
11 {node 0 org
coef , , , , , , , , , , , , , , , ,
include math.vf
here *cy =p
'r--- TO in_ '--l- TO out_ include hart2_.vf \ -- hh hl
'---u TO in_ '---u TO out_ include fourier2_.vf \ -- fh fl
node}
Аналогично, для остальных задействованных ядер после формирования таблицы, подключа-ются слова M* , D+.
Код для вычисления коэффициентов Хартли и Фурье примет вид
Файл hart2_.vf Код:
0 # dup dup
15 # for
dup dup xor dup . + drop \ очищаем бит переноса
\ -- sh sl i
in_ # b! @b \ -- sh sl i x приняли входной отсчёт
dup out_ # b! !b \ -- sh sl i x передали следующему
push \ -- sh sl i r: -- x
a! @a+ \ -- sh sl cas; a=i+1; r: -- x
pop a@ push \ -- sh sl cas x ; r: -- i+1
M* \ -- sh sl yh yl; r: -- i+1; где y=cas*x
D+
pop
next \ -- hvh hvl i+1 - отсчет Хартли спектра
drop \ -- hh hl
\ '-d-- # b! @b \ останов ядра
Файл fourier2_.vfКод:
\ -- hh hl \ - вычисляем Фурье спектр на основе Хартли отсчетов
$3fff0 # ~u/mod \ -- r q ( -- r hv )
in_ # b! @b \ -- hv h(N-v)
dup
M*
push push \ -- hv ; r: -- f2h f2l;
dup
M*
pop pop
D+
'-d-- # b! @b \ останов ядра
При таком подходе оперативная память ядер вычисляющих и H(v) и Фр(v) загружена на 80%.
Дамп памяти одного из ядер приведен ниже.
Код:
RAM Node 12
000 01000 ---
001 014E8
002 016A1
003 014E8
004 01000
005 008A9
006 00000
007 3F757
008 3F000 }коэффициенты
009 3EB18
00A 3E95F
00B 3EB18
00C 3F000
00D 3F757
00E 00000
00F 008A9 ---
: M*
010 134CA call CA *
011 308EA 2/ push drop .
012 26E12 pop a@ @p+ .
013 20000
014 39555 xor ;
: D+
015 2EA9A push a! pop .
016 3C88A + push a@ .
017 3CC55 + pop ;
018 04D97 @p+ dup dup @p+ -----------------------+
019 00000
01A 0000F
01B 2E9B2 push . . .
01C 24DE3 dup dup xor dup
01D 2C1EF . + drop @p+
01E 00175
01F 29F97 b! @b dup @p+
020 001D5 вычисление H(v)
021 29BBA b! !b push .
022 2BC9A a! @a+ pop .
023 228B2 a@ push . .
024 13410 call 10 M*
025 13415 call 15 D+
026 2701C pop next 1C
027 3BDB2 drop @p+ . . ----------------------------+
028 3FFF0
029 136AE call 2AE ~u/mod
02A 04B03 @p+ b! @b dup
02B 00145
02C 13410 call 10 M*
02D 2E892 push push dup . вычисление Фр(v)
02E 13410 call 10 M*
02F 26CB2 pop pop . .
030 13415 call 15 D+
031 04B00 @p+ b! @b +
032 00115 -----------------------------+
033
034
035
036
037
038
039
03A
03B
03C
03D
03E
03F
Временные параметры вычислений следующие:
- вычисление коэффициента Хартли спектра ~ 3110 тактов;
- вычисление Фурье ~ 380 тактов;
- общее время выполнения ~ 3490 тактов.
Скорость преобразования упала до 200000 преобразований в секунду.
Литература:
1. Брейсуэлл Р. Преобразование Хартли: пер. с англ. - М.: Мир, 1990.-175с.
2. SEAforth 40C18 Device Data Sheet. [Электронный ресурс] – Режим доступа: http://www.intellasys.net/templates/tri ... aSheet.pdf, свободный. – Яз. Англ.
3. IntellaSys - SEAforth 40C18. / http://www.intellasys.net/index.php?opt ... &Itemid=75.
4. New VentureForth Programmers Guide. / http://www.intellasys.net/index.php?opt ... &Itemid=68.
5. Александр Калачев. Процессоры семейства SEAforth. // Журн. Компоненты и технологии. – 2009. – №4. – С. 66-73.