Міністерство освіти і науки України
Національний університет «Львівська політехніка»
ІКТА
Кафедра ЕОМ
Курсовий проект
з курсу
«Методи, алгоритми, засоби цифрової обробки сигналів та зображень»
на тему:
«Розробка процессора ШПФ»
Завдання
Варіант № 27
Розробити процесор ШПФ з такими вхідними даними:
Тип процесора
ADSP-2181
Кількість точок
64
Основа ШПФ
4
Прорідження
часове
Розрядність вхідних даних
16
Такт поступлення вхідних даних
50 нс
Час обробки
1.8 мс
Анотація
В даному курсовому проекті розглянуто спосіб реалізації алгоритму ШПФ за основою 4 для сигнального процессора ADSP-2181 для 16-розрядних вхідних даних з часовим прорідженням, детально описано механізми обчислення швидкого перетворення Фур`є за заданною основою, властивості та основні характеристики процесора, на якому планується реалізація, підраховано часові ресурси для виконання обчислення, створена функціональна схема системи та написана програма, що реалізує вказаний алгоритм ШПФ на заданному процесорі.
Зміст
Вступ _______________________________________________5
Теоритичний розділ ____________________________________6
2.1. Характеристики сигнальногоо процесора ADSP2181___6
2.2. Опис ШПФ ______________________________________9
Аналіз блок-схеми виконання заданої
функції обробки сигналів та зображень
на заданому типі процесора ___________________________16
Розрахунковий розділ ________________________________19
Розробка функціональної схеми _______________________21
Розробка програми виконання алгоритму ШПФ __________24
Висновки ____________________________________________27
1.Вступ
Аналіз Фур'є закладає основи багатьох методів, що застосовуються в області цифрової обробки сигналів (ЦОС). По суті справи, перетворення Фур'є (фактично існує кілька варіантів таких перетворень) дозволяє співставити сигналу, заданому в часовій області, його еквівалентне представлення в частотній області, і навпаки, якщо відома частотна характеристика сигналу, то зворотне перетворення Фур'є дозволяє визначити відповідний сигнал у часовій області.
Крім того, ці перетворення корисні при проектуванні фільтрів. Частотна характеристика фільтра може бути отримана за допомогою перетворення Фур'є його імпульсної реакції. І навпаки, якщо визначена частотна характеристика сигналу, то необхідна імпульсна реакція може бути отримана за допомогою зворотнього перетворення Фур'є над його частотною характеристикою. Цифрові фільтри можуть бути створені на основі їхньої імпульсної реакції, оскільки коефіцієнти фільтра з кінцевою імпульсною характеристикою (КІХ) ідентичні дискретній імпульсній реакції фільтра.
2.Теоретичний розділ
2.1. Характеристики сигнального процесора ADSP-2181
2.1.1.Призначення ADSP-2181.
ADSP-2181 є однокристальним мікрокомп'ютером оптимізованим для обробка цифрового сигналу даних (DSP) і інших високошвидкісних задач обробки даних.ADSP-2181 об'єднує базову архітектуру сімейства ADSP-2100 (три обчислювальні блоки, генератори адреси даних ) з двома послідовними портами, внутрішнім 16-бітним DMA портом, байтовим DMA портом, програмованим таймером, ознакою I/O, пам'ять даних.ADSP-2181 інтегрує 80K байтів пам'яті розміщеної на кристалі, яка конфігурується як 16K слів (24-біти) програмної RAM, і 16K слів(16-біт) RAM даних.
2.1.2.Властивості ADSP-2181.
40 MІPS, 25 ns цикл виконання інструкцій, кожна інструкція виконується за
один цикл.
3-шинна архітектура дозволяє подвійну вибірку операндів вкожному
командному циклі
Багатофункціональні Інструкції
Код сімейства ADSP-2100, сумісний з розширенням множини інструкцій
Незалежний ALU, перемножувач/акумулятор, і зсувні обчислювальні блоки
Два незалежні генератори адреси даних
Зациклення по ознаці нуля
Умовне виконання команд
4 Мбайтний Інтерфейс Пам'яті для Зберігання Таблиць Даних
На рис.1 зображена блок-діаграма процесора ADSP-2181, на якій показані всі його основні вузли та зв’язки між ними. Основними вузлами є:
16-розрядні обчислювальні блоки – АЛП, вузол множення, зсувач;
генератори адреси даних (DAG1,DAG2);
програмний автомат (Program Sequencer) та регістр інструкцій;
таймер;
SRAM програм і SRAM даних;
зовнішні порти для зв'язку із зовнішньою пам'яттю та периферійними пристроями;
контролер DMA;
послідовні порти;
Логіка контролю за потужністю
Флаги
Програмовані I/O
Блок переривань
Рис.1.1. блок – діаграма процесора ADSP-2181
2.1.3.Опис виводів
Рис.1.2. Процесор ADSP-2181
В табл.1.1. використовуються такі позначення:
G – земля;
І – вхід;
О -- вихід;
P -- електроживлення;
* -- Ці виводи можуть під’єднуватись тільки до EZ-ICE коннекторів і використовуються тільки при емуляції;
Таблиця 1. 1 Опис виводів
Назва сигналу
I/O
Призначення
Address13-0
O
Адресні вихідні виводи для програмного простору, простору даних, байтового простору і простору вводу виводу
Data23-0
I/O
Вхідні/вихідні виводи даних для програмного простору і простору даних
I
Вхід скиду прочесора
I
Рівнево чутливий запит на переривання
I
Рівнево чутливі запити на переривання
I
Запит на шину
O
Дозвіл шини
O
Дозвіл шини
O
Вибір пам’яті програм
O
Вибір пам’яті даних
O
Вибір байтової пам’яті
O
Вибір простору в/в
O
Вибір комбінованої пам’яті
O
Дозвіл читання з пам’яті
O
Дозвіл запису в пам’ять
MMAP
I
Вибір карти пам’яті
BMODE
I
Завантажити вибіркові настройки
CLKIN,
XTAL
I
Тактова частота або кварц кристалу
CLKOUT
O
Тактова частота процесора
SPORT04-0
I/O
Послідовний порт в/в
SPORT14-0
I/O
Послідовний порт 1 в/в або два зовнішніх , ознака входу і ознака виходу
I
IDMA порт читання/запису
I
IDMA порт вибору
IAL
I
IDMA порт дозволу латчу адрес
IAD15-0
I/O
IDMA порт шини адресу/даних
O
IDMA порт підтвердження дозволу готовності
I
Контроль зниження потужності
PWDACK
O
Контроль зниження потужності
FL0,FL1,
FL2
O
Ознаки
PF7-0
I/O
Програмовані в/в
EE
*
Емулятор
*
Емулятор
*
Емулятор
*
Емулятор
*
Емулятор
*
Емулятор
ECLK
*
Емулятор
ELIN
*
Емулятор
ELOUT
*
Емулятор
GND10-0
G
Земля
VDD5-0
P
Живлення
2.2. Опис ШПФ
2.2.1.Опис швидкого перетворення Фур’є з прорідженням в часі
Дискретний матеріальний сигнал у вигляді кінцевої часової послідовності x(nТ) запишемо як x(nТ), де - число відліків, N – число відліків, T - період дискретизації.
N - точкове дискретне перетворення Фур'є (ДПФ) задається формулою:
де X(k) - частотний k-ий відлік чи k-а спектральна складова сигналу (визначає вихідну частотну послідовність, спектр сигналу),
комплексна експонента, W- ядро перетворення.
При зміні значення n*k на величину кратну N ядро не змінюється (у силу періодичності синуса і косинуса). Тобто ядро по верхньому індексу є періодичною функцією з періодом N. Тому замість добутку n*k можна вставити залишок від ділення його на N , тобто (n*k) mod N. Cпектральна функція X(k) також має період N по аргументу k.
Число множень дійсних відліків сигналу на комплексне ядро в (1) дорівнює N2, а число додавань комплексних чисел - (N -1)N. Кількість цих операцій різко зростає із збільшенням N і приводить до занадто великого часу перетворення.
ДПФ стало широко застосовуватися після винаходу швидких алгоритмів, в основі яких лежить принцип зведення багатоточкового перетворення до малоточкового. Один з них (що став уже класичним) називається ШПФ із проріджу-ванням за часом. Цей алгоритм отриманий за умови, якщо N є ступенем числа 2, тобто , де ν - ціле число, ν≥0.
Основні формули перетворення, до яких ми прийдемо, виходять при розбивці суми в (1) на дві суми, що містять відліки з парними і непарними номерами
Власне кажучи ця операція являє собою зміну порядку сумування (перегрупу-вання доданків), від якого сума не міняється.
Можна записати , . З врахуванням цього (2) прийме вигляд:
Зауважимо, що суми в (3) являють собою N/2 - точкові ДПФ часових відліків з парними номерами в першій сумі, і непарними номерами в другій сумі.
Позначимо, згідно з [1],
X(k) = Xν(k) - ДПФ усіх вхідних відліків дискретного сигналу,
вхідних відліків з парними номерами (перше число в нижньому індексі дорівнює ν - 1, а друге - парному числу (нулю)) ,
вхідних відліків з непарними номерами (друге число в нижньому індексі дорівнює непарному числу (одиниці)).
З урахуванням введених позначень одержимо коротку форму запису (3)
Спектри і є періодичними функціями з періодом N/2 (у ядрі перетворення замість N стоїть N/2). Величина називається повертаючим множником і володіє наступною цікавою властивістю
Ця властивість надає при обчисленні з номерами k від N/2 до (N -1) використовувати відповідні значення раніше обчислених з номерами від 0 до (N/2 -1) (потрібно тільки змінити знак).
За звичай формулу (4) розбивають на два співвідношення для зазначених вище номерів і одержують основні формули (базові співвідношення) перетворення Фур'є у вигляді
Базові співвідношення вже можна назвати швидким перетворенням Фур'є (ШПФ), тому що число операцій зменшилося. Число множень матеріальних відліків сигналу на комплексне ядро дорівнює . До цього потрібно додати множень комплексних чисел. Число додавань дорівнює
Процес розбиття за схемою базових співвідношень можна продовжувати до тих пір, доки не прийдемо до перетворень Фур'є одиничних відліків (що збігаються із самими відліками). Необхідне число операцій при цьому буде значно менше. У такому вигляді це ШПФ називають ШПФ із проріджуванням за часом.
Для пояснення процесу розбиття розглянемо більш детально, приклад ШПФ при .
Позначимо оператор ДПФ визначених вхідних відліків через F і запишемо відповідності між символами ДПФ, введеними вище, і цим оператором.
Зв'язок між ДПФ із великим і меншим числом точок відповідно до базових співвідношень записується в такий спосіб:
Роботу ШПФ можна пояснити графічно, якщо базові співвідношення, записані в дуже короткій формі
або зобразити графом
Верхній стрілці відповідає додавання, а нижній - віднімання. Попередньо bm-1 домножається на множник повороту WN.
Використовуючи вищенаведене і розташовуючи символи ДПФ у впорядковано-му вигляді, зобразимо ШПФ геометрично за допомогою графа.
Рис.1.2.1. Граф ШПФ із проріджуванням за часом при N=8
Відзначимо, що відліки вхідної послідовності переходять у відповідні ДПФ нульового рівня відповідно до інверсії їхній двійкових номерів (операція називається перестановкою вхідних відліків). Наприклад, десятковий номер 4|10 у двійковому вигляді запишеться як 100|2. Інверсія числа 100|2 (у прочитанні з права на ліво) запишеться як 001|2 = 1|10. Таким чином, вхідний відлік під номером 4 співпаде з першим ДПФ X0,1(0). Перестановку для всіх відліків можна показати стрілками переходу: 0 → 0, 1 → 4, 2 → 2, 3 → 6, 4 → 1, 5→ 5, 6 → 3, 7 → 7.
2.2.2.Алгоритм перетворення
Алгоритм ШПФ можна скласти, спираючи на граф ШПФ при N=8. Зауважимо, що ДПФ з однаковим числом точок на графі розташовані у вигляді стовпців. Пронумеруємо ДПФ у кожнім стовпці цифрами від 0 до 7 (від 0 до N-1 у загальному випадку) зверху вниз. Вихідні значення ДПФ - комплексні числа, які можна зберігати у деякому масиві. Перехід відповідно до базових співвідношень від ДПФ попереднього стовпця до ДПФ із подвоєним числом точок наступного стовпця назвемо кроком перетворення. З огляду на послідовний характер кроків перетворення, вихідні значення ДПФ кожного стовпця можна зберігати в тому самому масиві (це повинен бути масив комплексних чисел). Остаточно масив буде містити результуючі значення ШПФ.
Графи базових співвідношень на кожнім кроці візуально групуються. Групи складаються з окремих графів чи декількох графів, що перетинаються між собою. У прикладі на першому кроці графа є 4 групи, у другому – 2 і на третьому - 1.
Введемо позначення:
- число відліків сигналу, що обробляється, чи число точок перетворення;
v - число кроків перетворення;
step - номер кроку, step = 0,..., v - 1;
group - номер групи графів, group = 0, ..., (group_max - 1), де group_max - число груп (залежить від номера кроку);
іnput - номер ДПФ, вихід якого з'єднаний з верхнім входом графа базового співвідношення;
іnput +add - номер ДПФ, вихід якого з'єднаний з нижнім входом цього графа;
add - різниця відповідних номерів;
pow - ступінь множника повороту (у тексті pow відповідає показнику k в ) .
Змінні зв'язані між собою в такий спосіб:
Ці співвідношення отримані при аналізі графа ШПФ.
Для приведеного графа побудована таблиця, у якій для кожного кроку перетворення занесені індекси і їхні межі зміни, що використані в циклах програми. Нижче приведений алгоритм програми. Його особливість - результат виходить у масиві вхідних даних. Алгоритм побудований для випадку комплексних вхідних даних, як більш загальний випадок.
2.2.3.Алгоритм ШПФ із проріджуванням за часом
Вхідні дані в алгоритмі:
Х(k), k = 0,..., N -1 - масив комплексних вхідних і вихідних даних.
Для k = 0,..., N -1 встановити:
Встановити:
1. Перестановка елементів масиву Х(k).
2. Якщо step = v, перейти до кроку 9.
3. group←0/
4. Якщо group = group_max: step←step+1, add←add * 2,
group_max←group_max/2, перейти до кроку 2.
5. іnput ← group 2step +1.
6. Якщо іnput = group 2step+1+2step: group←group+1, перейти до кроку 4.
7. [Базова операція]
pow ←group_max *<input>2step+1,
t← X(іnput + add) * W(pow +1),
X(іnput +add) ←X(іnput) - t,
X(іnput) ←X(іnput) + t.
8. input←input+1, перейти до кроку 6.
9. Завершення.
2.2.4.Алгоритм двійково-інверсної перестановки
1. k←0.
2. Якщо k = N, то виконання алгоритму припиняється.
3. n← m← 0 (n←0,m←0 ).
4. Якщо m = v : перейти до кроку 7.
5. Якщо біт з номером m (0-ої біт - наймолодший) числа k дорівнює 1,
то n←n + 2v - m - 1.
6. m←m+1, перейти до кроку 4.
7. Якщо k < n : X(k)↔ X(n).
8. k←k+1, перейти до кроку 2.
3. Аналіз (розробка) блок-схеми виконання заданої функції обробки сигналів та зображень на заданому типі процесора
Алгоритм базової операції ШПФ за основою 4 і проріджування за часом можна представити у вигляді
А'1 = А1 + A2W1 + A3W2 + A4W3 = (А1 + A3W2) + (A2W1 + A4W3),
A'2 = A1 jA2W1 – A3W2 ± jA4W3 = (A1 – A3W2 ) j(A2W1 - A4W3 ),
A'3 = A1 - A2W1 + A3W2 - A4W3 = (A1 + A3W2) - (A2W1 + A4W3),
A'4 = A1 ± jA2W1 – A3W2 jA4W3 = (A1 - A3W2) ± j(A2W1 - A4W3),
де A'1, A'2, А'з, А'4 - результати базової операції; А1, А2, А3, А4 - вхідні відліки; W1, W2, W3 - комплексні коефіцієнти; j - уявна одиниця, верхній знак перед j відповідає прямому, нижній - оберненому ШПФ.
ReA'1 = [ReА1 + (ReA2*ReW2 - ІmA3*ImW2)] + [(ReA2*ReW1 - ImA2*ImW1) + (ReA4 *ReW3 – ImA4*ImW3)],
ІмA'1 = [Іm A1 + (Re A3*Im W2 + ImA3*Re W2)] + [(ReA2*ImW1 + ІmA2*ReW1) + (ReA4*ImW3 + ImA4*Re W3)],
RеA'2 = [ReA1 - (ReA3*ReW2 - ImA3*ImW2)] ± [(ReA2*ImW1 + ImA2*ReW1) – (ReA4*ІmW3 + ImA4*ReW3)],
ImA'2 = [ImA1 - (ReA3*ІmW2 + ImA3*ReW2)] [(ReA2*ReW1 - ImA2*ImW1) - (ReA4*ReW3 - ImA4*ImW3)],
ReA'3 = [ReA1 + (ReA3*ReW2 - ImA3*ImW2)] - [(ReA2*ReW1 - ImA2*ImW1) + (ReA4*ReW3 - ImA4*ImW3)],
ІmA'3 = [ІmA1 + (ReA3*ImW2 + ImA3*ReW2)] - [(ReA2*ImW1 + ImA2*ReW1) +(ReA4*lmW3 + ImA4*ReW3)],
ReA'4 = [ReA1 - (ReA3*ReW2 – ImA3*ImW2)] [(ReA2*ImW1 + ImА2*ReW1) - (ReA4*ReW3 -ImA4*ReW3)],
ImA'4 = [ImA3 - (ReA3*ImW2 + ImA3*ReW2)] ± [(ReA2*ReW1 - ImA2*Im W1) - (ReA4*ReW3 – ImA4*ImW3)]
Для виконання базової операції вимагається виконати 12 операцій множення і 22 додавання.
Порядок перед відповідною ітерацією
Номер метелика в ярусі
1
2
3
0
0
0
1
16
4
1
32
8
2
48
12
3
1
1
4
2
17
5
5
33
9
6
49
13
7
2
2
8
3
18
6
9
34
10
10
50
14
11
…
15
51
60
16
31
55
61
47
59
62
63
63
63
Блок-схема перетворення виглядає так:
Рис.2.1 Блок-схема 64-точкового перетворення за основою 4
Рис.2.2 Граф 64-точкового ШПФ за основою 4 з прорідженням по часу
4.Розрахунковий розділ
Частота роботи процесора: , звідси цикл виконання команди: .
base – основа базової операції «метелик»;
N – кількість точок вхідного перетворення;
base=4;
N=64;
– кількість етапів перетворення;
– кількість базових операцій «метелик» на одному етапі;
– кількість базових операцій у всьому перетворенні;
Для виконання базової операції «метелик» необхідно:
12 операцій множення;
22 операцій додавання;
14 операцій читання з пам`яті:
- 4*2=8 (для читання дійсної та уявної частини вхідних відліків);
- 3*2=6 (для читання дійсної та уявної частини комплексних коефіцієнтів);
8 операцій запису:
- 4*2=8 (для запису дійсної та уявної частини вхідних відліків);
В результаті на одну базову операцію припадає 56 операцій: Nна 1 мет=56 (оп)
Тривалість виконання обчислення ШПФ:
Тривалість надходження даних у процесор для обробки:
Тнадх=50нс – такт надходження даних;
Тривалість надходження даних у процесор та тривалість обчислення ШПФ:
Ця величина менша за заданий час обробки: , тобто для виконання обчислення достатньо одного процесора.
Внутрішня ОЗП процесора ADSP-2181 є 80Кб, і може бути сконфігурована наступним чином: 16К слів по 16 розрядів.
Для розв`язання поставленої задачі необхідно 64 слова по 16 розрядів, або 128 слів по 8 розрядів, тому внутрішньої ОЗП цілком вистачає.
Для зберігання повертаючих множників необхідно пам’ять об`ємом: 384х8,
оскільки на 1 операцію метелик необхідно 3 повертаючих множника, кожен з яких містить дійсну та уявну частину, всього є 64 базових операцій, розрядність операндів є 8. Тому 3*2*64=384.
Окремо використовуємо завантажувальну пам`ять, де буде зберігатися лістінг програми.
Необхідна пам’ять розміром 64х16 для постійного приймання даних, що надходять.
Рис.3.2.Ділянка пам`яті внутрішньої ОЗП, що приймає участь в обробці
5. Розробка функціональної схеми
5.1. Вибір мікропроцесора
Згідно варіанту завдання було вибрано процесор ADSP-2181, тактова частота роботи якого – 40МГц і відповідна тривалість такту – 25нс. На вхід сигнального мікропроцесора надходять такі сигнали:
CLKIN – сигнал синхронізації, що надходить з внутрішнього тактового резонатора;
– глобальний сигнал аппаратного скиду;
– сигнал зовнішнього переривання;
Вихідними та двонапрямленими сигналами по відношенню до мікропроцесора є:
– сигнал вибору кристалу мікросхеми пам`яті даних;
– сигнал вибору кристалу мікросхеми пам`яті констант;
ADDR[13:0] – шина адреси;
DATA[23:0] – шина даних;
– строб читання даних з зовнішнього пристрою у мікропроцесор;
– строб запису даних у зовнішній пристрій з мікропроцесора;
Інші сигнали або не задіяні, або їх використання не розглядається у данному курсовому проекті.
5.2. Пам’ять завантаження
Для завантаження програми використовується зовнішня пам’ятьь програм в яку попередньо завантажується програма і перед початком роботи з цієї пам’яті програма завантажується у внутоішню пам’ять програм. Для цього використовуються сигнали BMODE=0 і MMAP=1 і =1.Також в цій пам’яті зберігаються повертаючі множники.
5.3. Зовнішнє ОЗП
В якості зовнішнього ОЗП використовуєтьтся зовнішній RAM має 8Кб слів і доступ до нього йде коли =1.
5.4. Розробка керуючого пристрою
Призначення даного вузла – арбітраж доступу до зовнішнього ОЗП між обчислювальним процесором та давачем сигналу.
Сигнали, що надходять на керуючий пристрій:
CLK – сигнал синхронізації з давача;
– глобальний сигнал апаратного скиду;
STRD – строб даних, що надходить з давача;
ADDR_IN – шина адреси з обчислювального процесору;
– строб читання зовнішньої пам`яті, надходить з процесору;
Сигнали, що виходять з керуючого пристрою:
– дозвіл видачі даних з мікросхеми зовнішньої пам`яті;
– дозвіл запису до мікросхеми зовнішньої пам`яті;
– вибір кристалу зовнішньої пам`яті;
ADDR_OUT – шина адреси, що скеровується на мікросхему зовнішньої пам`яті;
– сигнал маскованого переривання;
Регістр стану зберігає значення сигналів , , , Вихід можна використати як сигнали дозволу у буферній розв`язці для шини даних, оскільки сигнали та взаємовиключаючі і суперечать один одному. Сигнал можна також використати як сигнал ACK, що надходить на процесор і є підтвердженням доступу до зовнішньої пам`яті зовнішнім пристроєм.
Рис.4.2. Структура керуючого пристрою
Рис.4.3. Структура обчислювальної системи
Табл..4.1. Таблиця станів керуючого автомату
Сигнал
Стан очікування
Запис в ОЗП
Читання з ОЗП
RD#
x
X
0
DMS#
x
x
0
STRD
0
1
0
OE#
1
1
0
WE#
1
0
1
CS#
1
0
0
MEM_RD
0
0
1
MEWR
0
1
0
ADDR_OUT
Z
ADDR
ADDR_IN
Сигнал формується лічильником при досягненні межі лічби і направляється на процесор, де обробляється програмою обробки переривань. Лічильник сам скидає сигнал переривання на початку нового циклу лічби. Лічильник сам скидається у початковий стан при надходження стробу даних, а також при подачі апаратного скиду.
Дані надходять із сенсора 16-розрядними і по черзі записуються у зовнішнє ОЗП
6. Розробка програми виконання алгоритму ШПФ
Структуру програми, що виконує обчислення за алгоритмом ШПФ можна уявити наступним чином:
Рис.5.1. Узагальнена блок-схема алгоритму
Кожен з трьох циклів призначений для правильного визначення номеру відліку в конкретний момент обчислення. Перший цикл визначає номер ярусу, другий – номер базової операції у ярусі, третій – номер відліку у базовій операції.
Вводиться масив, що зберігає відліки, в програмі названий matrix, його номер відповідно – N (кількість точок перетворення). Кожен елемент массиву – комплексне число. Інший массив W призначений для зберігання повертаючи множників. На кожну базову операцію припадає 3 повертаючих множника (четвертий фактично дорівнює 1), тому його розмір: 3*16*3=144 (3 - кількість ярусів, 16 - кількість базових операцій в ярусі). Елемент цього масиву є так само комплексним числом.
Текст програми, написаної на мові С, поданий нижче
N=64;
struct complex
{ double re;
double im;
};
complex W[3*N];
complex matrix[N];
int i,imax,j,x1,x2,x3,x4;
int f; //номер етапу
int p; //номер операції в етапі
complex temp1,temp2,temp3,temp4;
f=0;
for(imax = N-1; imax <=0; imax = (imax+1)/4 – 1)
{ p=0;
for(j = 0; j < N; j = j + (imax+1)*4)
{
x1 = j + (imax+1)*0;
x2 = j + (imax+1)*1;
x3 = j + (imax+1)*2;
x4 = j + (imax+1)*3;
for (i = 0; i <= imax; i++)
{
temp1.re=
matrix[x1+i].re + matrix[x3+i].re*W[n*64+p*3+1].re –
- matrix[x3+i].im*W[n*64+p*3+1].im +
+ matrix[x2+i].re*W[n*64+p*3+0].re –
- matrix[x2+i].im*W[n*64+p*3+0].im +
+ matrix[x4+i].re*W[n*64+p*3+2].re –
- matrix[x4+i].im*W[n*64+p*3+2].im;
temp1.im=
matrix[x1+i].im + matrix[x3+i].re*W[n*64+p*3+1].im +
+ matrix[x3+i].im*W[n*64+p*3+1].re +
+ matrix[x2+i].re*W[n*64+p*3+0].im +
+ matrix[x2+i].im*W[n*64+p*3+0].re +
+ matrix[x4+i].re*W[n*64+p*3+2].im +
+ matrix[x4+i].im*W[n*64+p*3+2].re;
temp2.re=
matrix[x1+i].re - matrix[x3+i].re*W[n*64+p*3+1].re +
+ matrix[x3+i].im*W[n*64+p*3+1].im +
+ matrix[x2+i].re*W[n*64+p*3+0].im +
+ matrix[x2+i].im*W[n*64+p*3+0].re -
- matrix[x4+i].re*W[n*64+p*3+2].im –
- matrix[x4+i].im*W[n*64+p*3+2].re;
temp2.im=
matrix[x1+i].im - matrix[x3+i].re*W[n*64+p*3+1].im -
- matrix[x3+i].im*W[n*64+p*3+1].re -
- matrix[x2+i].re*W[n*64+p*3+0].re +
+ matrix[x2+i].im*W[n*64+p*3+0].im -
- matrix[x4+i].re*W[n*64+p*3+2].re +
+ matrix[x4+i].im*W[n*64+p*3+2].im;
temp3.re=
matrix[x1+i].re + matrix[x3+i].re*W[n*64+p*3+1].re –
- matrix[x3+i].im*W[n*64+p*3+1].im -
- matrix[x2+i].re*W[n*64+p*3+0].re +
+ matrix[x2+i].im