Подякувати Студентському архіву довільною сумою
Admin
26.02.2023 12:38
Дякуємо, що користуєтесь нашим архівом!
Інформація про навчальний заклад
ВУЗ:
Національний університет Львівська політехніка
Інститут:
Не вказано
Факультет:
Не вказано
Кафедра:
Кафедра ЕОМ
Інформація про роботу
Рік:
2008
Тип роботи:
Навчальний посібник
Предмет:
Проектування комп’ютерних засобів обробки сигналів та зображень
Частина тексту файла (без зображень, графіків і формул):
Міністерство Освіти і Науки України
Національний Університет “Львівська Політехніка”
Кафедра ЕОМ
ЦИФРОВЕ ОПРАЦЮВАННЯ СИГНАЛІВ ТА ЗОБРАЖЕНЬ:
АЛГОРИТМИ ТА РЕАЛІЗАЦІЯ
Навчальний посібник
до лекційного курсу
з дисципліни “Проектування комп’ютерних засобів обробки сигналів та зображень” для студентів спеціальностей
7.091501 і 8.091501 "Комп'ютерні системи та мережі"
7.091503 і 8.091503 “Спеціалізовані комп'ютерні системи“
Львів – 2008
Цифрове опрацювання сигналів та зображень: Алгоритми та реалізація. Навчальний посібник до лекційного курсу з дисципліни “Проектування комп’ютерних засобів ообробки сигналів та зображень” для студентів спеціальностей 7.091501 і 8.091501 "Комп'ютерні системи та мережі", 7.091503 і 8.091503 “Спеціалізовані комп'ютерні системи“ / Укладачі: Р. Попович, Є. Ваврук – Львів: Національний університет “Львівська політехніка”, 2008, 147 с.
Укладачі: Є. Ваврук, к.т.н., доцент
Р. Попович, к.т.н., доцент
Відповідальний за випуск: Мельник А. О., професор, завідувач кафедри
Рецензенти: Ткаченко Р.О., д. т. н, професор
Голембо В.А. В.В., к. т. н, доцент
Анотація
Даний навчальний посібник укладений відповідно до навчальної програми з дисципліни "Проектування комп’ютерних засобів обробки сигналів та зображень". В ньому розглянуті основні питання проектування алгоритмічних, програмних та апаратних засобів опрацювання сигналів та зображень, вивчення яких допоможе студентам спеціальностей 7.091501 і 8.091501 /"Комп'ютерні системи та мережі"/ та 7.091503 і 8.091503 /"Спеціалізовані комп'ютерні системи"/ здобути початкові практичні навички у оптимальному виборі алгоритмів, їх програмній та апаратній реалізації.
ЗМІСТ
Вступ
5
1. Особливості задач і алгоритмів цифрового опрацювання сигналів та зображень
7
2. Алгоритми швидкого перетворення Фур’є та їх програмна реалізація
13
3. Організація DSP- процесорів для задач опрацювання сигналів та зображень
26
4: Інтерфейси DSP-процесорів
30
5. Проектування процесора ШПФ на ПОС
38
6: Проектування засобів опрацювання сигналів та зображень на ПЛІС
53
7. Реалізація алгоритмів опрацювання сигналів та зображень на нейропроцесорах
62
8. Стиск нерухомих зображень з використанням перетворень різного типу дискретних косинусних перетворень
71
9. Опрацювання мовних сигналів
95
10. Використання вікон для опрацювання сигналів
122
11. Діагностика і контроль процесорів і систем опрацювання сигналів та зображень
135
Висновки
145
Література
146
ВСТУП
В останні роки цифрове опрацювання сигналів та зображень (ЦОСЗ) суттєво впливає на такі ключові технологічні галузі, як телекомунікації, цифрове телебачення і засоби інформації, біомедицина і цифровий звукозапис. Сьогодні ЦОСЗ є основою багатьох новітніх видів цифрових розробок і різних застосувань в інформаційному середовищі (наприклад. цифровий мобільний зв'язок, цифрові відеокамери, телебачення і системи звукозапису). Поряд з тим ЦОСЗ ширше впроваджується і в класичні галузі застосування (радіолокація, геофізика, опрацювання мовних сигналів, сейсмологія, системи зв'язку і передачі даних, медична і технічна діагностика, системи керування). Такому широкому застосуванню сприяли успіхи в розробці швидких алгоритмів і методів ЦОСЗ (наприклад, хвилькові перетворення), досягнення в технології конструювання інтегральних схем (архітектура нових процесорів ЦОСЗ, програмовані логічні інтегральні схеми – ПЛІС, системи на кристалі), використання програмних пакетів (MATLAB, база програм на мові С).
Суть ЦОСЗ як області науки пролягає у розв'язку на обчислювальній машині чотирьох основних задач:
- представлення сигналів та зображень в зручній для сприйняття формі;
- виділення із сигналів та зображень корисної інформації;
- внесення в сигнали та зображення корисної інформації;
- формування сигналів та зображень із заданими параметрами.
Основні переваги ЦОСЗ полягають в: можливості реалізації складних методів та алгоритмів опрацювання сигналів та зображень, недоступних для аналогових пристроїв; забезпеченні високої точності опрацювання; гнучкості і універсальності засобів, розвиненому користувацькому інтерфейсові тощо.
Головна проблема ЦОСЗ полягає у підвищенні швидкодії при реалізації певного набору математичних операцій над сигналами та зображеннями. Цей набір, як правило, визначається проблемною орієнтацією засобів ЦОСЗ. Для розв'язання цієї задачі комплексно застосовуються два напрямки - підвищення ефективності обчислювальних алгоритмів і вдосконалення архітектури комп'ютерних засобів. При цьому виходять з детального аналізу особливостей проблемної області. Власне, незалежно від вибраного напрямку і теми, студенти повинні виконувати курсові роботи з даного курсу виходячи з головної проблеми ЦОСЗ.
Метою вивчення дисципліни є засвоєння основних методів та алгоритмів опрацювання сигналів та зображень, принципів та шляхів проектування апаратних і програмних комп’ютерних засобів опрацювання сигналів та зображень, набуття початкових практичних навиків проектування таких засобів.
В результаті вивчення дисципліни студенти повинні:
знати: типи і особливості алгоритмів опрацювання сигналів та зображень, математичні основи їх реалізації, характеристики елементної бази та основи реалізації алгоритмів на її основі, особливості і етапи розробки цифрових пристроїв на процесорах і вузлах, орієнтованих на задачі опрацювання сигналів та зображень; мови програмування та мови опису апаратних засобів; пакети моделювання, опрацювання і фільтрації сигналів та зображень
вміти: вибирати найефективніший метод і алгоритм, здійснювати математичне формування алгоритму та методу його розв’язання, розробляти блок-схеми алгоритмів, розробляти програмну реалізацію виконання алгоритмів, розраховувати технічні параметри апаратних і програмних засобів, проектувати структурні та функціональні схеми процесорів та вузлів опрацювання сигналів; моделювати і опрацьовувати сигнали різної форми.
Для тих, хто хоче краще орієнтуватися в перспективних напрямках розвитку комп’ютерних засобів опрацювання сигналів та зображень, краще розуміти специфіку і особливості їх проектування радимо, насамперед, скористатись матеріалами Internet-сайтів
1. Особливості задач і алгоритмів цифрового опрацювання сигналів та зображень
Аналіз задач і алгоритмів
До основних галузей, де використовується опрацювання сигналів та зображень відносяться:
1. Радіолокація (РЛ) — виявлення, фільтрація сигналу з режекцією завад та накопичення сигналу.
2. Гідроакустика (ГА) — контроль водяного простору, підводна сейсмологія, гідронавігація.
3. Зв’язок (ЗВ) — підвищення надійності, пропускної здатності, скритності, виділення символів, згортка символьних послідовностей, скорочення надлишковості, підвищення завадостійкості, керування транспортними засобами.
4. Геофізика (ГФ) — пошук нафтоносних (водоносних) шарів.
Біомедицина (БМ) — візуалізація органів, діагностика.
Системи керування виробничими процесами — керування, знімання даних, передача інформації.
7. Промислова діагностика — неруйнівний контроль, візуалізація стану вузлів, діагностика віддалених об’єктів.
Кількість основних алгоритмів і операцій, що використовуються для їх виконання, як це видно з табл.1.1 є невеликою.
Таблиця 1.1. Перелік основних алгоритмів і застосованих операцій
Алгоритми обробки і
Області застосування
застосовувані операції
РЛ
ГА
ЗВ
ГФ
БМ
СК
Згортка
+
+
+
ШПФ
+
+
+
+
+
Перемноження елементів векторів
+
+
+
+
Додавання елементів векторів
+
+
+
Перетворення координат
+
+
+
Обчислення тригонометричних функцій
+
+
+
+
+
Пошук максимальних значень
+
+
Операції над матрицями
+
+
+
+
Табличне перетворення відліків зображення
+
Сортування
+
+
Отримання квадратного кореня
+
+
+
+
+
Піднесення до степеня
+
+
+
+
+
Статистична обробка
+
Нормування
+
+
Розрахунок нормованих коефіцієнтів кореляції
+
Операції повертання координат
+
Логарифмування
+
Потенціювання
+
Номінальні розкиди координат
+
Інтерполяція даних
+
+
Розглянемо основні типи алгоритмів.
Лінійна фільтрація. Лінійна фільтрація використовується для подавлення шумів, спектр яких не перетинається зі спектром сигналу і дозволяє виділити в чистому вигляді функцію, що описує явище, яке досліджується. На практиці розрізняють два основних типи лінійних фільтрів: з скінченою (СІХ) і нескінченою (НІХ) імпульсними характеристиками. НІХ фільтри мають рекурсивну структуру, що описується різницевими рівняннями згідно з формулою (1.1).
, n ( 0; (1.1)
де x та y - вхідна та вихідна реалізації процесу; ai і bi - постійні коефіцієнти, що характеризують властивості фільтра, причому aм(0; M – порядок фільтру.
СІХ - фільтри реалізуються на основі нерекурсивної структури, яка описується формулою (1.2), і може бути отримана з формули (1.1) якщо ai прирівняти до нуля,.
, n ( 0; (1.2)
причому, bi можна назвати дискретною імпульсною характеристикою системи (фільтра).
Медіанна фільтрація. Медіанна фільтрація є нелінійним способом опрацювання одномірних та двомірних послідовностей (виборок) і використовується для зменшення рівня імпульсних шумів. В порівнянні з лінійною, медіанна фільтрація має такі переваги: зберігає різкі перепади сигналу і добре згладжує імпульсний шум. Процес виконання фільтрації цього типу виконується в три етапи: сортування виборок реалізації в ковзаючому «вікні»; вибір середнього значення у «вікні» (медіана); заміна відліку, розташованого в середині вікна значенням медіани.
Дискретне перетворення Фур’є. Дискретне перетворення Фур’є (ДПФ) скінченої періодичної послідовності відліків сигналу {х(n)}, 0 n N-1, визначається формулою (1.3):
,k=0,1,…,N-1 (1.3)
де - комплексний дискретний ортонормований базис.
Таким чином, ДПФ представляє собою множину спектральних коефіцієнтів Фур’є, що відповідають N рівномірно рознесеним частотам, починаючи від нульової і до (але не включно) частоти 2/T .Обернене дискретне перетворення Фур’є (ОДПФ) визначається формулою (1.4):
,n=0,1,…,N-1; (1.4)
ОДПФ дозволяє однозначно перейти від Фур’є-образу X(k) до функції х(n).
Швидке перетворення Фур’є (ШПФ) є узагальненою назвою множини алгоритмів, що призначені для виконання ДПФ і ОДПФ. Основна ідея алгоритмів ШПФ полягає в розбитті процедури виконання ДПФ на m=log2N етапів, на кожному з яких вхідний дискретний набір розділяється на дві частини. До кожної такої частини вхідного дискретного набору, застосовується операція ДПФ в 2 рази меншої розмірності: в результаті застосувань розбиття m разів на останньому етапі, матимемо тривіальне 2-х точкове ДПФ. Застосування ШПФ можливе лише тоді, коли N є складним числом. Найбільшого розповсюдження набули алгоритми ШПФ для послідовностей довжини N, що є степенем числа 2 (N=2m). Загальне число операцій в алгоритмі ШПФ складає приблизно Nlog2N додавань і 0,5Nlog2N множень комплексних чисел. Алгоритми ШПФ з проріджуванням в часі і з проріджуванням по частоті принципово мало відрізняються між собою. Тому розглянемо тільки алгоритм ШПФ з проріджуванням в часі. Формула розкладу алгоритму, в результаті застосування якої отримується шуканий дискретний набір X(k), k=0,1,...,N-1, що є зображенням вихідного дискретного набору перетворюваної функції x(n), n=1,2,...,N-1 в Фур’є просторі визначається згідно з формулою (1.5).
, k=0,1,...,N-1; (1.5)
де X0(k) та X1(k) відповідно N/2 - точкові ДПФ послідовностей x(2n) та x(2n+1), n=0, 1, ..., (N/2-1). При проріджуванні, N-точкове ДПФ зводиться до обчислення двох N/2-точкових ДПФ. Рекурсивно застосовуючи вказаний спосіб розділення до перетворень меншої розмірності, отримуємо алгоритм ШПФ за основою два з проріджуванням в часі. Аналізуючи формулу (1.5), і беручи до уваги, що послідовності X0(k) та X1(k), періодичні з періодом N/2, та WNN-k = - WNk , можна встановити, що при кожному розділенні реалізується однотипна базова операція:
A*=A+WNkB; B*=A – WNk B; (1.6)
причому k{0, 1, ..., (N/2-1)}, а величини A, B, A*, B*, і WNk - комплексні числа. Для виконання базової операції (1.6) необхідно виконати одне множення, одне додавання і одне віднімання комплексних величин.
Кореляція. Для визначення подібності між сигналами в різні моменти часу або виділення сигналу на фоні шуму проводять кореляційну обробку, важливе місце в якій займає обчислення функцій взаємної кореляції двох числових послідовностей x і y, кожна з яких містить N відліків, записується у вигляді:
, r = 0,1,…,N-1; (1.7)
Тобто, взаємна кореляційна функція двох сигналів обчислюється як сума їх поелементних добутків з відносною затримкою r.
Автокореляційна функція записується у вигляді:
, r = 0,1,…,N-1; (1.8)
Обчислення автокореляційної і взаємокореляційної функцій двох сигналів складається з трьох основних операцій: часової затримки, множення і сумування.
Згортка сигналів. При згортці сигналів x(n) та h(m), де n = 0,1,…,N1-1, m = 0,1,…,N2-1, виконується обробка згідно з формулою (1.9).
; (1.9)
де h(m) і x(n) рівні нулю поза відповідними інтервалами, h(m) - дискретна імпульсна характеристика системи (фільтра).
Нормування сигналів. Одним з прикладів операції нормування сигналів є ділення комплексних чисел, що виконується згідно з формулою (1.10).
(1.10)
Сортування. В основі алгоритмів сортування лежать дві основні операції: порівняння і пересилання даних.
Більшості з названих алгоритмів властиві регулярність, рекурсивність і локальність, що робить їх придатними для реалізації на НВІС. Зокрема, до локально рекурсивних алгоритмів відносяться алгоритми премноження матриць, згортки, НІХ-фільтрація, сортування вибором; до глобально-рекурсивних - алгоритми ШПФ, бітове сортування.
На основі аналізу алгоритмів в табл. 1.2 виділено набір базових операцій для опрацювання сигналів та зображень.
Таблиця 1.2. Базові операції для задач керування та опрацювання інформації
Алгоритм
Базові операції
Реалізація фільтрів, обчислення кореляційної і взаємокореляційної функцій
Додавання, віднімання, множення, обчислення суми парних добутків
Алгоритми перетворення
Множення, додавання комплексних чисел; обчислення операцій ШПФ, ШПХ, множення послідовностей комплексних чисел
Алгоритми обчислення: коефіцієнтів взаємної кореляції, координат, відстаней, виконання вагового множення, тощо
Обчислення тригонометричних функцій, добування квадратного кореня, піднесення до степені, ділення, сортування
Детальніше конкретні алгоритми описані і реалізовані в подальших розділах посібника.
Основні параметри типових операцій. Основні параметри типових операцій наведені в табл. 1.3, де N - розмірність масиву інформації.
Таблиця 1.3. Параметри типових операцій
Тип операції
Складність виконання
Тип обчислюва-льних засобів
Необхідна швидкодія (MOPS)
Лінійні операції: просторова фільтрація, згортка, виявлення контурів, скалярний добуток, КІХ-фільтрація.
N
Скалярний
102 - 105
Операції другого порядку: лінійні перетворення, перетворення Фур'є, кореляція, перемноження матриці на вектор, сортування, медіанна фільтрація.
N2
Векторний
103 - 107
Операції вищих порядків: матричні операції, спектральні обчислення, адаптивні операції, розв’язання задач лінійних алгебраїчних рівнянь
N3
Матричний
104 - 108
Особливості задач і алгоритмів.
Аналіз наведених задач і алгоритмів їх розв’язання показує, що вони мають такі особливості:
- широкий динамічний і частотний діапазон сигналів, що обробляються;
- велика інтенсивність і постійність потоків даних;
- великий об’єм обчислень з перевагою обчислювальних операцій над логічними;
- регулярність і рекурсивність алгоритмів;
- можливість розпаралелення алгоритмів в часі та в просторі;
- в основі більшості алгоритмів лежить математична операція множення/сумування;
- постійне ускладнення нових алгоритмів і підвищення вимог до точності даних;
- структура даних дозволяє застосувати їх векторне опрацювання з використанням просторового та часового паралелізму.
Особливості організації обчислювальних засобів
1.2.1. Методи аналізу обчислювальних засобів архітектур. Технічно системи керування та опрацювання інформації реалізуються як комплекс спеціалізованих і універсальних засобів обчислювальної техніки. Головними вимогами при їх проектуванні є адаптація структури на виконання задач в РРЧ і забезпечення вагогабаритних показників, споживаної потужності та оптимального співвідношення між обсягами апаратних засобів та програмного забезпечення. Такий підхід забезпечує найкращі швидкісні та вагогабаритні параметри вузлів, однак робить їх частково вузькоспеціалізованими. Тому при проектуванні та аналізі систем вибирається метод, який забезпечує виконання вищенаведених вимог, і одночасно дозволяє реалізувати апаратні та програмні засоби на базі універсальних вузлів, або можливості оперативного переналаштовування спеціалізованих вузлів на виконання інших задач. При синтезі бортових систем користуються змістовними (евристичними, інтуїтивними) та формальними (алгоритмічними) методами, більшість з яких орієнтована на не взаємозв’язані моделі систем. Такі моделі відображають різні боки системної організації, але не враховують цілісність систем в функціональному і структурному аспектах та концепцію їх розвитку. В цьому напрямку розглядаються два підходи до проектування систем: структурно-функціональний (СФП) і функціонально-структурний (ФСП). Перший напрямок зводиться до розбиття складних систем на окремі структурні рівні (підсистеми) для визначення їх функціонального призначення. Він є обмеженим, як концепція синтезу систем, що не дозволяє ефективно використовувати його при створенні систем цифрового опрацювання. ФСП розглядає питання розвитку та взаємодії програмних і апаратних засобів, аналізу функцій систем і передбачає: цілісний підхід до аналізу і синтезу багаторівневих систем; облік взаємозв'язку функцій і структури об'єктів (при визначаючій ролі функції по відношенню до структури); облік енергетичних і інформаційних зв'язків як між елементами системи, так і системи з зовнішнім середовищем; - розгляд системи в розвитку.
1.2.2. Процесори в задачах опрацювання сигналів та зображень. Структура системи визначається в основному специфікою задач і алгоритмів, що вирішуються комплексом. Загальний варіант структурної схеми СКОІ наведений на рис.1.1, де Д - давачі, П - приймач, ВАО - вузол аналогової обробки, ВД - вузол дискретизації, ОВ - обчислювальний вузол, ВК - вузол керування, ВР - вузол реєстрації, відображення, додаткового опрацювання, ВП - виконавчий пристрій, - дискретна вибірка. Блок давачів виконує функції перетворення вхідних даних i (t) в вихідні сигнали S1(i,,t), де - вектор неінформаційних параметрів сигналів. В приймачі виконується попереднє опрацювання суміші сигналів і завад (t)=S1(i,,t)+n1(t), різні типи фільтрації зовнішніх n1(t) і внутрішніх n2(t) завад, нормування вихідних процесів і аналогово-цифрове перетворення x(t). Функції інших вузлів залежать від області застосування системи. Причому, частина обчислювальних функцій і алгоритмів виконується вузлами ОВ і ВК. В даній структурі вузли ВД, ОВ, ВК, ВР можуть бути реалізовані на базі універсальних процесорів, процесорів опрацювання сигналів (DSP-процесори), ПЛІС, систем на кристалі.
Шляхи та засоби підвищення продуктивності засобів цифрового опрацювання наведені в табл..1.4.
Напрямки реалізації апаратних засобів на базі яких проектуються засоби опрацювання є такими:
- програмовані процесори;
- процесори на базі комплектів спеціалізованих НВІС функціональних вузлів;
- спеціалізовані НВІС процесорів, які виконують конкретні алгоритми цифрового опрацювання сигналів (ЦОС);
- процесори на базі НВІС програмованих середовищ;
- процесори ЦОС (ПОС), які конструктивно вмонтовані в ЕОМ;
- процесори на базі ПЛІС.
Таблиця 1.4. Шляхи та засоби підвищення продуктивності
Шляхи підвищення ефективності
Засоби підвищення ефективності
Інші
Підвищення регулярності структури
Розширення функціональних можливостей
Підвищення продуктивності
—
—
Розпаралелення процесу обчислення, апаратно-мікропрограмна реалізація укрупнених операторів
Мінімізація втрат при обміні інформацією
—
Режим розподілу ре-сурсів, укрупнення операторів
—
Зменшення складності апаратних засобів
Зменшення складності зв’язків за рахунок викорис-тання НВІС
—
—
Зменшення склад-ності апаратних засобів за рахунок збільшення об’ємів виробництва
Використання модульного конструювання
Інваріантність до типу ЕОМ, збільшення числа операторів, що ефективно реалізуються, використання багатофунк- ціональних модулів
—
Рис.1.1.
2. Алгоритми швидкого перетворення Фур’є та їх програмна реалізація
Оскільки алгоритми ШПФ є одними з найуживаніших при опрацюванні сигналів їх розгляд є доцільним в даному посібнику.
Алгоритм швидкого перетворення Фур'є – є оптимізованим за швидкодією способом обчислення дискретного перетворення Фур'є (ДПФ), що має складність O(Nlog2N) на відміну від складності ДПФ порядку O(N2).
2.1. Основні положення алгоритму ШПФ
Визначення 1. Дано кінцеву послідовність x0, x1, x2,..., xN-1 (у загальному випадку комплексних чисел). ДПФ полягає в пошуку послідовності X0, X1, X2,..., XN-1, елементи якої обчислюються за формулою:
(2.1)
Визначення 2. Зворотне ДПФ полягає в пошуку послідовності x0, x1, x2,..., xN-1, елементи якої обчислюються за формулою:
(2.2)
Основною властивістю перетворень (2.1) і (2.2) є те, що з послідовності {x} отримується (при прямому перетворенні) послідовність {X}, а якщо застосувати до {X} зворотне перетворення, то знову отримується вихідна послідовність {x}.
Визначення 3. Величина
називається повертаючим множником.
Властивості повертаючих множників
При k = 1
Пряме перетворення Фур'є можна виразити через повертаючі множники. Тоді формула (2.1) матиме вигляд:
(2.3)
Геометричне тлумачення повертаючих множників наведене на рис.1. Комплексне число (rejφ) представлене у вигляді вектора, що виходить із початку координат (r - модуль числа, а φ – аргумент). Модуль відповідає довжині вектора, а аргумент - куту повороту. Модуль повертаючого множника . дорівнює одиниці, а фаза - 2π/N. Оскільки при множенні комплексних чисел, представлених у показниковій формі, їхні модулі перемножуються, а аргументи підсумовуються, множення вихідного числа на повертаючий множник не змінить довжину вектора, але змінить його кут. Тобто, відбудеться повертання вектора на кут 2π/N.
З формули (2.3) можна визначити геометричний зміст перетворення Фур'є таким чином: представити N комплексних чисел-векторів з набору {x}, кожне у вигляді суми векторів з набору {X}, повернених на кути, кратні 2π/N.
Рис.2.1.
Основні формули
Теореми, що пояснюють суть перетворення Фур’є (наведені без доведення).
Теорема 1. Якщо комплексне число представлене у вигляді e j2πN, де N - ціле, то e j2πN = 1.
Теорема 2. Величина періодична по k і по n з періодом N. Тобто, для будь-яких цілих l і m виконується рівність:
(2.4)
Теорема 3.
Для величини справедлива формула:
(2.5)
З наведених теорем визначається основна ідея алгоритму ШПФ:
Необхідно розділити суму (1) з N доданків на дві суми по N/2 доданків, і обчислити їх окремо. Для обчислення кожної з підсум, треба їх теж розділити на дві і т.д.
Необхідно повторно використовувати вже обчислені доданки.
При обчисленні алгоритму ШПФ застосовують або "проріджування за часом" (в першу суму попадають доданки з парними номерами, а в другу - з непарними), або "проріджування за частотою" (в першу суму попадають перші N/2 доданків, а в другу - інші). Обидва варіанти рівноцінні. В силу специфіки алгоритму доводиться застосовувати тільки N, що є ступенями 2. Розглянемо випадок проріджування за часом.
Теорема 4. Визначимо ще дві послідовності: {x[even]} і {x[odd]} через послідовність {x} таким чином:
x[even]n = x2n, x[odd]n = x2n+1, (2.6)
n = 0, 1,..., N/ 2-1
Нехай до цих послідовностей застосовані ДПФ і отримані результати у вигляді двох нових послідовностей {X[even]} і {X[odd]} по N/2 елементів у кожній.
Елементи послідовності {X} можна виразити через елементи послідовностей {X[even]} і {X[odd]} за формулою:
(2.7).
Формула (2.7) дозволяє скоротити число множень удвічі (не враховуючи множень при обчисленні X[even]k і X [odd]k), якщо обчислювати Xk не послідовно від 0 до N - 1, а попарно: X0 і XN/2, X1 і XN/2+1,..., XN/ 2-1 і XN. Пари утворяться за принципом: Xk і XN/2+k.
Теорема 5. ДПФ можна обчислити і за формулою:
(2.8)
З теореми випливає, що можна не зберігати обчислені значення X[even]k і X[odd]k після обчислення чергової пари, і одне обчислення можна використовувати для обчислення двох елементів послідовності {X}.
На цьому кроці буде виконане N/2 множень комплексних чисел. Якщо застосувати подібні схеми для обчислення послідовностей {X[even]} і {X[odd]}, то кожна з них вимагатиме N/4 множень, разом ще N/2. Продовжуючи аналогічно далі log2N разів, можна дійти до сум, що містять лише один доданок, так що загальна кількість множень рівна (N/2)log2N.
Розглянемо ШПФ для різних N. Додамо ще один нижній індекс, який буде вказувати на загальну кількість елементів послідовності, до якої цей елемент належить. Тобто X{R}k - це k-ий елемент послідовності {X{R}} з R елементів. X{R}[even]k - це k-ий елемент послідовності {X{R}[even]} з R елементів, обчислений по парних елементах послідовності {X{2R}}. X{R}[odd]k - це k-ий елемент послідовності {X{R}[odd]}, обчислений по непарних елементах послідовності {X{2R}}.
У випадку, коли доданок лише один (N = 1) формула (1) спрощується до:
,
Оскільки в цьому випадку k може бути рівне тільки 0, то X{1}0 = x{1}0, тобто, ДПФ над одним числом дає це ж саме число.
Для N = 2 за теоремою (2.5) отримаємо:
Для N = 4 за теоремою (2.5) отримаємо:
Звідси виходить, що якщо елементи вихідної послідовності були дійсними, то при збільшенні N елементи ДПФ стають комплексними.
Для N = 8 за теоремою (2.5) отримаємо:
Необхідно звернути увагу, що на попередньому кроці використовувалися ступені W4, а на цьому - ступені W8. Зайвих обчислень можна уникнути, якщо врахувати, що
Тоді формули для N=4 будуть використовувати ті ж W-множники, що і формули для N=8:
Висновок:
В основі алгоритму ШПФ лежать такі формули:
(2.9)
2.2. Програмна реалізація основних елементів ШПФ
Алгоритм попередньої перестановки
Розглянемо конкретну реалізацію ШПФ. Нехай є N=2T елементів послідовності x{N} і треба одержати послідовність X{N}. Розділимо x{N} на дві послідовності: парні і непарні елементи і т.д з кожною послідовністю. Ітераційний процес закінчиться, коли залишаться послідовності довжиною по 2 елементи. Приклад процесу для N=16 показаний нижче:
Разом виконується (log2N)-1 ітерацій.
Розглянемо двійкове представлення номерів елементів і займаних ними місць. Елемент із номером 0 (0000) після всіх перестановок займає позицію 0 (0000), елемент 8 (1000) - позицію 1 (0001), елемент 4 (0100) - позицію 2 (0010), елемент 12 (1100) - позицію 3 (0011). І так далі. Позиції до і після перестановок дзеркально симетричні. Двійкове представлення кінцевої позиції виходить із двійкового представлення початкової позиції перестановкою бітів у зворотному порядку. І навпаки.
Це є закономірним для всіх N. На першому кроці парні елементи з номером n перемістилися в позицію n/2, а непарні з позиції в позицію N/2+(n-1)/2, де n=0,1,..., N-1. Таким чином, нова позиція обчислюється зі старої позиції як:
ror(n,N) = [n/2] + N{n/2},
де [x] - ціла частина числа, {x} - дробова.
В асемблері ця операція - циклічний зсув вправо (ror), якщо N - ступінь двійки. Спочатку береться двійкове представлення числа n, потім всі біти, крім молодшого зміщуються на 1 позицію вправо, а молодший біт переміщається на місце, що звільнилося, найстаршого (лівого) біта (див. рис.2.2).
Рис. 2.2
Подальші розбиття виконуються аналогічно. На кожному наступному кроці кількість послідовностей подвоюється, а число елементів у кожній з них зменшується вдвічі. Операції ror піддаються вже не всі біти, а тільки кілька молодших. Старші j-1 бітів залишаються недоторканими (зафіксованими), де j - номер кроку (див. рис.2.3):
Рис. 2.3
Простежимо за довільним бітом номера позиції. Нехай біт перебував в j-му двійковому розряді. Біт буде послідовно зсуватися вправо на кожному кроці доти, поки не виявиться в крайній правій позиції (після j-го кроку). На j+ 1-му кроці буде зафіксовано j старших бітів, а той біт, за яким стежимо, переміститься в розряд з номером T-j-1, що і необхідно для дзеркальної перестановки біт.
Алгоритм перестановки бітів наведений нижче:
for(I = 1; I < N-1; I++)
{
J = reverse(I,T); //reverse переставляє біти в I у зворотному порядку
if (I >= J) // пропустити вже переставлені
conitnue;
S = x[I]; x[I] = x[J]; x[J] = S; // перестановка елементів xI і xJ
}
Операція зворотної перестановки біт може бути реалізована через інші бітові операції. Нижче наведений алгоритм функції перестановки T молодших бітів у числі I:
unsigned int reverse(unsigned int I, int T)
{
int Shift = T - 1;
unsigned int LowMask = 1;
unsigned int HighMask = 1 << Shift;
unsigned int R;
for(R = 0; Shift >= 0; LowMask <<= 1, HighMask >>= 1, Shift -= 2)
R |= ((I & LowMask) << Shift) | ((I & HighMask) >> Shift);
return R;
}
У змінних LowMask і HighMask зберігаються маски, що виділяють два біти, які переставляються. Перша маска у двійковому поданні виглядає як 0000...001 і в циклі змінюється, зсуваючи одиницю щораз на 1 розряд вліво:
0000...001
0000...010
0000...100
...
Маска HighMask за кожну ітерацію зсуває одиничний біт на 1 розряд вправо.
1000...000
0100...000
0010...000
...
Ці зсуви виконуються інструкціями LowMask <<= 1 і HighMask >>= 1.
Змінна Shift показує відстань (у розрядах) між бітами, що переставляються. Спочатку вона дорівнює T-1 і кожну ітерацію зменшується на 2. Цикл припиняється, коли відстань стає менше або дорівнює нулю.
Операція I & LowMask виділяє перший біт, потім він зсувається на місце другого (<<Shift). Операція I & HighMask виділяє другий біт, потім він зсувається на місце першого (>>Shift). Після чого обидва біти записуються в змінну R операцією "|".
Замість того щоб переставляти біти позицій місцями, можна застосувати і інший метод. Для цього треба вести відлік 0,1,2,...,N/ 2-1 уже зі зворотним проходженням бітів. Алгоритм збільшення на одиницю можна реалізувати програмно, приклад для T=4 наведений нижче.
I = 0;
J = 0;
for(J1 = 0; J1 < 2; J4++, J ^=1)
for(J2 = 0; J2 < 2; J3++, J ^=2)
for(J4 = 0; J4 < 2; J4++, J ^=4)
for(J8 = 0; J8 < 2; J8++, J ^=8)
{
if (I < J)
{
S = x[I]; x[I] = x[J]; x[J] = S; // перестановка елементів xI і xJ
}
I++;
}
У алгоритмі враховано, що при збільшенні числа від 0 до нескінченності (зі збільшенням на одиницю) кожний біт міняється з 0 на 1 і назад з певною періодичністю: молодший біт - щораз, наступний - кожний другий раз, наступний - кожний четвертий і так далі. Ця періодичність реалізована у вигляді T вкладених циклів, у кожному з яких один з бітів позиції J перемикається туди і назад за допомогою операції XOR (В C/C++ записується як ^=). Позиція I використовує звичайний інкремент I++.
Даний алгоритм залежить від кількості вкладених циклів залежно від T. Оскільки T звичайно обмежено деякою розумною межею (16..20), то можна написати декілька варіантів алгоритму. Нижче наведений варіант цього алгоритму, що емулює вкладені цикли через стеки Index і Mask.
int Index[MAX_T];
int Mask[MAX_T];
int R;
for(I = 0; I < T; I++)
{
Index[I] = 0;
Mask[I] = 1 << (T - I - 1);
}
J = 0;
for(I = 0; I < N; I++)
{
if (I < J)
{
S = x[I]; x[I] = x[J]; x[J] = S; // перестановка елементів xI і xJ
}
for(R = 0; R < T; R++)
{
J ^= Mask[R];
if (Index[R] ^= 1) // еквівалентно Index[R]^=1; if (Index[R] != 0)
break;
}
}
Величина MAX_T визначає максимальне значення для T і в найгіршому випадку дорівнює розрядності цілочисельних змінних. Даний алгоритм повільнішим від попереднього, але економніший за обсягом коду.
І ще один алгоритм, який використовує класичний підхід до багаторозрядних бітових операцій: треба розділити 32-біта на 4 байти, виконати перестановку в кожному з них, після чого переставити самі байти.
Перестановку бітів в одному байті можна робити по таблиці, для якої потрібно мати масив reverse256 з 256 елементів (записані числа від 0 до 255 з переставленими в кожному порядок бітів).
Даний масив можна застосувати для реалізації функції reverse:
unsigned int reverse(unsigned int I, int T}
{
unsigned int R;
unsigned char *Ic = (unsigned char*) &I;
unsigned char *Rc = (unsigned char*) &R;
Rc[0] = reverse256[Ic[3]];
Rc[1] = reverse256[Ic[2]];
Rc[2] = reverse256[Ic[1]];
Rc[3] = reverse256[Ic[0]];
R >>= (32 - T);
Return R;
}
Звертання до масиву reverse256 переставляють у зворотному порядку біти в кожному байті. Покажчики Ic і Ir дозволяють звернутися до окремих байтів 32-бітних чисел I і R і переставити у зворотному порядку байти. Зсув числа R вправо наприкінці алгоритму усуває розходження між перестановкою 32 біт і перестановкою T біт. Геометрична ілюстрація алгоритму, де стрілками показані перестановки бітів, байтів і зсув наведена на рис.2.4.
Рис. 2.4
У всіх випадках є N перестановок з порівнянням I і J, що запобігає повторній перестановці деяких елементів.
Перевагу необхідно надати останньому алгоритму - він має усього N переходів.
Оптимізований останній алгоритм наведено нижче.
static unsigned char reverse256[]=
{
0x00, 0x80, 0x40, 0xС0, 0x20, 0xA0, 0x60, 0xE0,
0x10, 0x90, 0x50, 0xD0, 0x30, 0xB0, 0x70, 0xF0,
0x08, 0x88, 0x48, 0xC8, 0x28, 0xA8, 0x68, 0xE8,
0x18, 0x98, 0x58, 0xD8, 0x38, 0xB8, 0x78, 0xF8,
0x04, 0x84, 0x44, 0xC4, 0x24, 0xA4, 0x64, 0xE4,
0x14, 0x94, 0x54, 0xD4, 0x34, 0xB4, 0x74, 0xF4,
0x0C, 0x8C, 0x4C, 0xCC, 0x2C, 0xAC, 0x6C, 0xEC,
0x1C, 0x9C, 0x5C, 0xDC, 0x3C, 0xBC, 0x7C, 0xFC,
0x02, 0x82, 0x42, 0xC2, 0x22, 0xA2, 0x62, 0xE2,
0x12, 0x92, 0x52, 0xD2, 0x32, 0xB2, 0x72, 0xF2,
0x0A, 0x8A, 0x4A, 0xCA, 0x2A, 0xAA, 0x6A, 0xEA,
0x1A, 0x9A, 0x5A, 0xDA, 0x3A, 0xBA, 0x7A, 0xFA,
0x06, 0x86, 0x46, 0xC6, 0x26, 0xA6, 0x66, 0xE6,
0x16, 0x96, 0x56, 0xD6, 0x36, 0xB6, 0x76, 0xF6,
0x0E, 0x8E, 0x4E, 0xCE, 0x2E, 0xAE, 0x6E, 0xEE,
0x1E, 0x9E, 0x5E, 0xDE, 0x3E, 0xBE, 0x7E, 0xFE,
0x01, 0x81, 0x41, 0xC1, 0x21, 0xA1, 0x61, 0xE1,
0x11, 0x91, 0x51, 0xD1, 0x31, 0xB1, 0x71, 0xF1,
0x09, 0x89, 0x49, 0xC9, 0x29, 0xA9, 0x69, 0xE9,
0x19, 0x99, 0x59, 0xD9, 0x39, 0xB9, 0x79, 0xF9,
0x05, 0x85, 0x45, 0xC5, 0x25, 0xA5, 0x65, 0xE5,
0x15, 0x95, 0x55, 0xD5, 0x35, 0xB5, 0x75, 0xF5,
0x0D, 0x8D, 0x4D, 0xCD, 0x2D, 0xAD, 0x6D, 0Xed,
0x1D, 0x9D, 0x5D, 0xDD, 0x3D, 0xBD, 0x7D, 0xFD,
0x03, 0x83, 0x43, 0xC3, 0x23, 0xA3, 0x63, 0xE3,
0x13, 0x93, 0x53, 0xD3, 0x33, 0xB3, 0x73, 0xF3,
0x0B, 0x8B, 0x4B, 0xCB, 0x2B, 0xAB, 0x6B, 0xEB,
0x1B, 0x9B, 0x5B, 0xDB, 0x3B, 0xBB, 0x7B, 0xFB,
0x07, 0x87, 0x47, 0xC7, 0x27, 0xA7, 0x67, 0xE7,
0x17, 0x97, 0x57, 0xD7, 0x37, 0xB7, 0x77, 0xF7,
0x0F, 0x8F, 0x4F, 0xCF, 0x2F, 0xAF, 0x6F, 0xEF,
0x1F, 0x9F, 0x5F, 0xDF, 0x3F, 0xBF, 0x7F, 0xFF,
}
unsigned int I, J;
unsigned char *Ic = (unsigned char*) &I;
unsigned char *Jc = (unsigned char*) &J;
for(I = 1; I < N - 1; I++)
{
Jc[0] = reverse256[Ic[3]];
Jc[1] = reverse256[Ic[2]];
Jc[2] = reverse256[Ic[1]];
Jc[3] = reverse256[Ic[0]];
J >>= (32 - T);
if (I < J)
{
S = x[I];
x[I] = x[J];
x[J] = S;
}
}
Основний цикл алгоритму
На рис.2.5 наведений другий етап обчислення ДПФ. Лініями зверху вниз показано використання елементів для обчислення значень нових елементів. Зручно, що два елементи на певних позиціях в масиві дають два елементи на тих самих місцях. Тільки належати вони будуть вже іншим, більш довшим масивам, що розміщені на місці попередніх (коротших). Це дозволяє обійтись одним масивом даних для вихідних даних, результату і зберігання проміжних результатів для всіх T ітерацій.
Рис. 2.5
Нижче наведені дії, які необхідно виконати після первинної перестановки елементів.
#define T 4
#define Nmax (2 << T)
complex Q, Temp, j(0,1);
static complex x[Nmax];
static double pi2 = 2 * 3.1415926535897932384626433832795;
int N, Nd2, k, mpNd2, m;
for(N = 2, Nd2 = 1; N <= Nmax; Nd2 = N, N += N)
{
for(k = 0; k < Nd2; k++)
{
W = exp(-j*pi2*k/N);
for(m = k; m < Nmax; m += N)
{
mpNd2 = m + Nd2;
Temp = W * x[mpNd2];
x[mpNd2] = x[m] - Temp;
x[m] = x[m] + Temp;
}
}
}
Змінна Nmax містить повну довжину масиву даних і кінцеву кількість елементів ДПФ. В таблиці наведена відповідність між виразами в формулі (2.9) і змінними в програмі.
В алгоритмі
В формулі
N
N
Nd2
N/2
k
K
m
k з поправкою на номер послідовності
mpNd2
k+N/2 з поправкою на номер послідовності
pi2
2π
j
j (уявна одиниця)
x[k] (на початку циклу)
X{N/2}[even]k
x[mpNd2] (на початку циклу)
X{N/2}[odd]k
x[k] (в кінці циклу)
X{N}k
x[mpNd2] (в кінці циклу)
X{N}N/2+k
W
Зовнішній цикл - це основні ітерації. На кожній из них 2Nmax/N ДПФ (довжиною по N/2 елементів кожне) перетворюється в Nmax/N ДПФ (довжиною по N елементів кожне).
Наступний цикл по k є циклом синхронного обчислення елементів з індексами k і k + N/2 у всіх результуючих ДПФ.
Внутрішній цикл перебирає Nmax/N штук ДПФ одне за другим: спочатку обчислюються елементы з номерами 0 и N/2, у всіх ДПФ в кількості Nmax/N штук, потім з номерами 1 і N/2 + 1 і т.д. На рис.2.5 показана послідовність обчислень для N = 8, в якій забезпечується однократне обчислення .
x[mpNd2] обчислюється раніше від x[k], щоб x[k] не було модифіковано перед використанням.
Алгоритм оберненого ДПФ відрізняється тим, що замість -j треба використовувати +j і після всіх обчислень поділити кожне x[k] на Nmax.
Оптимізація повертаючих множників
Для обчислення повертаючих множників можна використати формулу (10), що дозволить
. (2.10)
уникнути операцій піднесення до ступеня і обійтися одними множеннями, якщо заздалегідь обчислити WN для N = 2, 4, 8,…,Nmax. У міру збільшення k повертаючий множник буде змінюватися, але разом з тим буде рости похибка обчислень. І після N/2 підряд множень у повертаючому множнику, може нагромадитися велике відхилення від точного значення, оскільки при множенні величин їхні відносні похибки сумуються.
Цього можна було б уникнути, будь число WN цілим, але воно не ціле при N > 2, тому що обчислюється через значення синуса і косинуса:
Уникнути багатьох звертань до повільних функцій синуса і косинуса можна за допомогою алгоритму обчислення ступеня через багаторазові піднесення до квадрату і множення. Наприклад:
У цьому випадку необхідно всього 5 множень ( не потрібно обчислювати двічі) замість 13. У найгіршому випадку для піднесення до ступеня від 1 до N/ 2-1 потрібно log2N множень замість N/2, що дає цілком прийнятну похибку для більшості практичних задач.
Можна вдвічі зменшити кількість множень на кожному кроці, якщо використовувати результати минулих обчислень
,
для зберігання яких потрібно додатково log2(Nmax) комплексних комірок пам'яті:
Якщо чергове не було обчислено попередньо, то береться двійкове представлення k і аналізується. Кожному одиничному біту відповідає рівно один множник. У загальному випадку одиниці в біті з номером b відповідає множник , що зберігається в b-ій комірці масиву.
Для зменшення кількості множень при обчисленні до одного на два цикли необхідно відвести N/2 комплексних комірок для зберігання всіх попередніх . Непарні елементи обчислюються за формулою (2.10). Парні обчислюються за формулою , тобто, нічого не обчислюється, а береться одне зі значень, обчислених на попередньому кроці, коли N було вдвічі менше. Щоб не потрібно було копіювати величини на нові позиції досить їх відразу розташовувати в тій позиції, що вони займуть при N = Nmax і вводити просте виправлення Skew.
Пояснення до лістингу, що наведений нижче.
1. Реалізовані найпростіші операції над комплексними числами (класи Complex і ShortComp...
Ділись своїми роботами та отримуй миттєві бонуси!
0%
Оголошення від адміністратора