CC BY
32
Обчислювальні методи в радіоелектроніці
УДК 621.372.061
УМОВНА ДЕКОНВОЛЮЦІЯ В ОБЛАСТІ ТРАНСФОРМАНТ ФУР’Є. ПОБУДОВА ОБЕРНЕНОЇ МАТРИЦІ ДЕГРАДАЦІЇ ОБРАЗУ
Рибін О.І., д.т.н. професор, ІванюкН.О., аспірантка
Національний технічний університет України «Київський політехнічний інститут», м. Київ, Україна
Вступ
Реставрація образів має велике значення при розв’язанні багатьох задач сучасної техніки і медицини. Серед методів реставрації найбільш простим і достатньо ефективним є метод умовної деконволюції, при реалізації якого однією з умом є те, що енергія шуму в реставрованому образі не повинна перевищувати енергії шуму образу вихідного (спотвореного) [1—3]. Реставрацію за методом умовної деконволюції звичайно виконують в області перетворення Фур’є [1,3]. В роботах [4,5] запропоновано алгоритм реставрації образу за методом умовної деконволюції в базисі перетворення Ада-мара. Показано [5] можливість використання символьних методів для формування матриць, математичної моделі фільтра, який виконує процедуру реставрації. Так, згідно з [4,5] оцінку fHad спектру Адамара реставрованого образу можна знайти за виразом
/Н а d
G2 Н а d + 1 C
2 Had
xG1HHad x^Had,
(1)
де LHad — дискретний спектр Адамара спотвореного сигналу з адитивним шумом розміру N*2n,n — ціле позитивне число; G[Had — спектральний матричний оператор перетворення Адамара, обчислюваний відповідно формулі G[Had = HadxGxHadТ; G2Had та C2Had спектральні матриці Адамара порядку N імпульсної характеристики деградації образу та функції «гладкості» відповідно, обчислювані за формулами G2Had = HadxGT X GxHadT; C2Had=HadxCT XCxHadT;G, C—матричні оператори дискретної згортки та «гладкості» відповідно, порядку N; Had — матричний оператор дискретного перетворення Адамара; X — коефіцієнт варіації Лагранжа; Т — знак транспонування.
Формування матриць G2Had, C2Had, G1Had можливо виконати за допомогою символьного методу, запропонованого в [5], що дозволяє підвищити точність обчислень при операціях з обмеженою розрядність операндів. Г о-ловну задачу становить обернення матриці в (1).
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№50
21
Обчислювальні методи в радіоелектроніці
Обернення матриці для реалізації алгоритму реставрації
Матриці G, С є циркулянтними, що й забезпечує блочно-діагональний вигляд матриць G2Had, C2Had, G±Had [5] та діагональний вигляд матриць перетворення Фур’є G2F, C2Had від G, С.
Так, нормований (діленням кожного рядка матричного оператора F на VN) матричний оператор Фур’є FH можна представити у вигляді добутку двох матричних операторів
F =Ртх Had; F* = HadхР*, (2)
де Had — нормований матричний оператор дискретного перетворення Адамара; Рт — нормований матричний оператор кратного [6 — 8] ортогонального перетворення; * — знак комплексного спряження. З (2) отримуємо
Рт = FH xHadT; Р*= HadxF*; Рт хр*= р*хрт = Е,
(3)
де Е — одинична матриця.
Тоді, якщо знайти прості символьні залежності для формування елементів матриці Р, то для достатньо прозорих і легко програмованих для безпосереднього формування матриць G2Had, C2Had знайти символьні залежності добутків
В = Рт X G2Had хР‘, (4)
які мають вигляд діагональних матриць (зручних і для обернення в символьному вигляді), не буде складно.
Результат обчислення матриці Рт має просту блочно-діагональну структуру, що дозволяє у подальшому в символьному вигляді реалізувати обчислення за виразом (1). Покажемо це на прикладі підматриць-блоків діагоналі невисоких порядків. Так, перший блок діагоналі (з номерами рядків та стовпців 0,0) має елемент, що дорівнює 1. Другий елемент діагоналі (з номерами рядків та стовпців 1,1) теж дорівнює 1. Відзначене вірно, якщо порядок матриці Ртдорівнює двом. Якщо порядок матриці Рт дорівнює чотирьом, другий елемент діагоналі буде в першому рядку, але вже в другому стовпці (з номером 2). Для матриці восьмого порядку елемент першого рядка знаходиться вже в (2x2=4) в стовпці з номером 4 і т.д. Третій елемент діагоналі — це матриця другого порядку, яка має вигляд (якщо порядок матриці Рт дорівнює чотирьом)
Рт
(2,3),(1,3)
COSU/Ax
e~jn/4
ejn/4
ejn/4 '
e-jn/4_ ,
тобто в цьому випадку ненульові елементи 2, 3 рядків знаходяться в стовпцях 1, 3. Якщо порядок матриці збільшити вдвічі, то номери стовпців теж
22
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№50
Обчислювальні методи в радіоелектроніці
збільшаться вдвічі (2, 6). Тому у подальшому будемо наводити мінімальні значення номерів стовпців, що відповідає випадку, коли відповідні рядки мають номери, максимальні для даного формату матриці РТ .
Так підматриця четвертого порядку
^Т(4 . . . 7),(1,3,5,7)
С О S7T/4 X
cos п/8 х e~J3n/8 cos Зп/8 x e-n78 cos 3n/8 x en78 cosn/8 x ej3n/8
j sirni/8 х е~]3п/8 /s/пЗп/8 x e~jn/8 —jsin Зп/8 x e^/8 —jsinn/8x e]3n/8 (5)
jcos п/8 x e~]3n/8 jcos 3n/8x e~]n/8 —jcos Зп/8 x e^/8 —jcosn/8 x e3n!8
—sinп/8 x e~]3n/8 —sin 3n/8x e~]n/8 —sin Зт/8 x e^/8 —s іпт/ 8 x е і37г/8 .
Аналогічно, підматриця восьмого порядку має елементи першого рядка (множник с о s т/4 винесено за знак матриці)
р8,1= соs т/ 1 6хсоs 2 т/ 1 6хе~і7п/1 6; р83= соs 3 т/ 1 6 хсо s 6т/ 1 6хе_і 5я/ 1 6; р8,5= — со s 5 т/ 1 6 хсо s 1 От/ 1 6хе_і 3я/ 1 6;р8,7= — со s 7т/ 1 6 хс о s 1 3 т/
1 6хе^я/ 1 6;р8,9= р8,7; р8,1 1 = Р8,5; Р 8,1 3 = Р8 3 Р8,1 5= Ps, 1-Елементи першого стовпця (починаючи з р 9 , l)
Р9 , і. = / s in т / 1 6хс о s 2 т / 1 6хе_ 7я/ 1 6; р1 о , -і = / с о sт / 1 6 хs in 2 т/ 1 6хе_ 7я/ 1 6; р1 1 і. = -s і пт / 1 6хs in 2 т/ 1 6хе_;” 7я/ 1 6; р1 2 , l = /со s^ 1 6хсо s2 т/ 1 6хе_;'77Г/ 1 6; р1 3 , l= -s іпт/ 16хс о s 2 т/ 1 6хе_J'7я/ 16; р1 4 , l = — с о sт / 1 6 хs і п2 т / 1 6хе_ 7я/ 1 6; р1 5 , l = — s in т / 1 6хs і п 2 т / 1 6хе_ 7я/ 1 6.
Усі елементи кожного з рядків (8...15) мають ті самі функції, що й перший елемент рядка, але їх аргументи відповідно т/ 1 6 ,2 т/ 16; 3т/16, 6т/16; 5т/16, 10т/16; 7т/16, 14т/16; 7т/16, 14т/16; 5т/16, 10т/16; 3т/16, 6т/16; т/ 1 6 ,2 т/ 16. Знаки перед функціями отримуємо множенням знака (1, -1, /, -/) перед елеме н том першого стовпця н а 1 , 1 , -1, -1, -1, -1, 1, 1, а аргументи при експоненті змінюються від / (з кроком до
+/ 7 т/ 16. Знаки перед функціями першого стовпця для перших двох рядків 1 та відповідно, а для другої пари рядків утворюються множенням першої пари на . Тригонометричні функції вздовж по стовпцю змінюються наступним чином. Множник з найменшим аргументом (т / 1 6 ) послідовно вздовж стовпця, починаючи з міняється на , знову на
і т.д. Множник з подвоєним аргументом є функцією, яка повторюється двічі вздовж стовпця, тобто , ,
, , і т.д.
Вісник Національного технічного університету України "КПІ" 23
Серія — Радіотехніка. Радіоапаратобудування.-2012.-№50
Обчислювальні методи в радіоелектроніці
Аналогічно, для матриці РТ 32-го порядку матриця Рт (1 6 _ _ _ з 1 ) ( 1 ( з ( 5 ( 7> 9> 1 1 ( 1 з ( 1 5 _ _ _ _ 3 1 ) формується безпосередньо за знайденими вище закономірностями (номери стовпців ненульових елементів, обчислені для матриці вдвічі меншого формату подвоюються). Для цього побудуємо перший елемент першого стовпця додаткової (для подвоєного формату) підматриці (множник с о s л/4 виносимо за знак матриці)
Рі6і= с о s л/ 3 2хс о s 2 л/ 3 2хс о s 4л/3 2 е~і 1 5я/3 2.
Далі вздовж кожного з рядків аргументи перших множників добутку тригонометричних функцій мають в чисельниках множники перед , які дорівнюють 3, 5, 7, ... NI2-1, NI2-1,...7, 5, 3, 1. Друга тригонометрична функція в добутку має подвоєний аргумент, третя — помножений на 4 і т.д. Експоненти мають аргументи, починаючи від -(N/2-1), які зростають з кроком 2 л / N до значення +(N/2-1), де N — порядок матриці Рт (в даному випадку N=32). Знаки перед добутком тригонометричних функцій вздовж рядка до стовпця з номером NI2-1 змінюються наступним чином: +1, +1, -1, -1, +1, +1, -1, -1, ... і т.д. Починаючи зі стовпця з номером NI2, усі множники слід замінити на протилежні. Ці знаки (вздовж рядка) слід для різних рядків (вздовж стовпця) помножити на 1, (перша пара рядків) наступну пару рядків слід помножити на , -1, (що відповідає множенню першої пари знаків на У). Отриману четвірку знаків слід помножити на j і отримати множники , -1, -1, - . Отриману вісімку знаків слід знову помножити на j і т.д. Тригонометричні функції вздовж кожного стовпця чере-дуються наступним чином. Функції з мінімальним аргументом чередують-ся за кожним кроком, тобто якщо в першому рядку підматриці першим множником є функція с о s( 2 п — 1 ) /N, то наступною є s і п( 2 п — 1 ) / N, а потім знову с о s( 2 п — 1 ) / N. Функції, починаючи з с о s2 ( 2 п — 1 ) / N череду-ються що два кроки, а функції, починаючи з с о s 4( 2 п — 1 ) / N — що чотири кроки і т.д.
Отже, запропонований символьний алгоритм формування матриці Рт є достатньо простим. Переставленням стовпців матриця Рт зводиться до блочно-діагональної. Тоді кожен блок діагоналі цієї матриці та спряженої транспонованої можна множити праворуч та ліворуч на відповідний до
нього блок діагоналі матриць(22Яай, С2Яad, ^іяad. Результатом такого множення буде діагональна підматриця. Розглянемо тепер можливості безпосереднього формування елементів діагональної матриці D.
Обчислення діагональної матриці
Оскільки добуток (4) дає діагональну матрицю, то у подальшому будемо обчислювати лише діагональні елементи. Обчислення діагональних елементів dj jматриці D для кожного блока діагоналі в (4) можна організувати за виразом
24
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№50
Обчислювальні методи в радіоелектроніці
dU= [If *{(А*Х Pt) О Сгнагіблока) X[1] _ (б)
де Pt*, Pj — і-ті стовпець та рядо к відповідн их бло ків матриць Р* та Р; Т — зн а к тра нспо нува н ня; О — зн а к м н оже н н я, та кого, при я кому коже н елеме н т матриці Р/хpf множиться на відповідний йому елемент матриці G 2 я a d б ло ка; * — знак комплексного спряження; [ 1 ] — стовпець одиниць.
Так, для підматриць другого четвертого порядку (5) матимемо для першого діагонального елемента
Я* х Pf = [ СOST/ 8 X eJ3т/8 — _/S ІП7Г/ 8 X eJ3т/8 —/ СOST/ 8 X eJ3 Т/8 — S ІП7Г/
8 X e]'37T/8]f X [ С OST/ 8 X e—3 71/8 JS І П Т/8 X e—3?т/8у С О ST/ 8 X e—3 71/8 — SІП7Г/
COS 7Г/8 X COS 7Г/8 —jsinn/Q х cosn/Q —jcosn/Q х cosn/Q —sinn/Q x cosn/Q
j COS n/Q X sin 7T/8 sinn/Q x sinn/Q cos n/Q x sinn/Q —jsinn/Q x sinn/Q
j cos n/Q x cos 7t/8 sinn/Q X COS7T/8 COS 7T/8 X COS 7T/8
—jsinn/Q x cosn/Q
8 X e—3т/8]=
— COS 7T/8 x sin n/Q
-jsinn/Q x Sinn/Q jcos n/Q x sinn/Q sinn/Q x sinn/Q
(7) '
Оскільки для блока G 2 я a d б ло ка (7) елементи G ! , 2, G ! , з, G 2,1, G з , !, G 2,4, G3,4, G4,2, G4, з дорівнюють нулю, матрицю (7) можна записати у вигляді
п п
cos — X cos — 8 8
О
О
п п
—sin — х cos — 8 8
О
п п
sin— х sin — 8 8
п п
cos — х sin — 8 8
О
О
п п sin — х cos — 8 8 п п cos — X cos — 8 8
О
п П~і
— cos — х sin — 8 8
О
О
п п
sin— х sin —
8 8
(8)
Нехай для визначеності імпульсну характеристику деградації образу описано вісімкою дискретних відліків , тоді, блок четвертого порядку матриці G 2 я має вигляд
де Gi, і= G2,2 Gl,4
1 2 Я a d бло ка
G4 4 0 0 - Gl,4 ’
0 G2,2 Gl,4 0
0 <7і,4 Gl,l 0
— Gl,4 0 0 G2,2 -
-/5 + 0+Р 7;
+ ^+0-/7;
0+/? і + 0—/?з+0—/?5+ 0+/? ^
(9)
причому /о = S=09^ /і = Sk=o9kX9k+і; /2 = Sk=o9kX9k+2;...
З=о 9k X 9s - 1 = 9k X 9s - і; S<N [9].
Тоді результат обчислень за виразом (6) має вигляд dі, і= /о+V2х/1+0 -V2х/3-/4 -V2х/5 + 0 +V2х/7;
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№50
25
Обчислювальні методи в радіоелектроніці
d44 dlt±; d33- /0+0 +0 +0 - /4+0 +0 +0; d2і2— d33 і т . д.
Аналогічно, для матриці восьмого порядку.
Для ілюстрації запропонованого алгоритму розглянемо сигнал рис. 1, а, спектр потужності Фур’є якого зображено на рис. 1,6.
1 2 3 4 5 6 7 8
1 2 3 4 5 6 7 8
Рис. 1, а
Рис. 1,б
Для сигналу рис.1,а матриця D в (4) має вигляд
1,8557 0 0 0 0 0 0 0
0 0,7981 0 0 0 0 0 0
0 0 1,1161 0 0 0 0 0
0 0 0 1,1161 0 0 0 0
0 0 0 0 1,454 0 0 -0,2416
0 0 0 0 0 0,9708 0,2416 0
0 0 0 0 0 0,2416 1,454 0
0 0 0 0 0,2416 0 0 0,9708
Значення елементів цієї діагональної матриці (при відповідному впорядкуванні трансформант рис. 2) співпадають зі складовими спектру потужності Фур’є сигналу рис.1,а.
Усі обчислення проводилися з використанням Mathlab.
1 5 8 4 2 3 6 7
Рис. 2
26
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№50
Обчислювальні методи в радіоелектроніці
Як видно з вище приведених рисунків, спектри сигналу, отримані різними шляхами, співпадають з розумною похибкою (табл.1).
Таблиця 1
Спектр потужності сигналу по Фур’є Спектр потужності сигналу по Адамару (матриця D )
1,0 1,0
0,8375 0, 7835
0,6015 0,6015
0,4692 0,5232
0,4301 0,4301
0,4692 0,5232
0,6015 0,6015
0,8375 0, 7835
Для обернення матриці D достатньо обчислити зворотні значення для елементів діагоналі, а для обчислення зворотної матриці (G2Had) — пом-
ножити матрицю (D) на РТ праворуч та на Р* ліворуч. Певні відміни між деякими трансформантами спектрів потужності, обчислених різними способами, пояснюється різними способами обчислень та відносно невеликим форматом їх представлення (з використанням перетворення Адамара або Фур’є). Зі збільшенням формату ці відміни стають знакомо малими.
Висновки
1. Використання символьних залежностей при формуванні матриці G2Had потребує мінімально можливої кількості операцій множення та додавання.
2. Використання інформації про жорстку структуру одержуваних матриць довільного порядку дозволяє обчислювати лише їх ненульові елементи, що значно зменшує кількість арифметичних операцій, усуває дублювання , які виникають при безпосередньому перемноженню матриць перемноженні матриць.
3. Розроблений алгоритм, легко програмується та дозволяє з достатньо високою точністю обчислити обернену матрицю будь-якого порядку.
Література
1. Jan Jin Cislicovafiltrace, Analyzaarestauracesignlu / VUTvBRNE,1997,438s.
2. Рибін О.І.Реставрація образів методом умовної деконволюції в області просторових частот / О.І. Рибін, В.Ю. Корольов // Вісник Технічного університету Поділля.— 2000.— C. 145 — 147.
3. Рыбин А.И. Реставрация образов в частотной области методом взвешенной фильтрации / А.И. Рыбин, В.Ю. Королев // Радиоэлектроника — 2001.— №4.— С. 51 — 56 (Изв. вузов).
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№50
27
Обчислювальні методи в радіоелектроніці
4. Рибін О.І. Реставрація образів за методом умовної деконволюції в області трансформант перетворення Адамара/ О.І. Рибін, Н.О. Іванюк // Вісник НТУУ «КПІ». Сер. Радіотехніка. Радіоапаратобудування. — 2011. — №46.— С. 51 — 58.
5. Рибін О.І.Умовна деконволюція в області трансформант Фур’є. Побудова матриці деградації образу / О.І. Рибін, Н.О. Іванюк // Вісник НТУУ «КПІ». Сер. Радіотехніка. Радіоапаратобудування. — 2011. — №47.— С. 30 — 41.
6. Рибін О.І. Аналіз лінійних систем в області трансформант кратного перетворення EIWAL / О.І. Рибін, А.П. Ткачук // Вісник НТУУ «КПІ». Сер. Радіотехніка. Радіоапаратобудування.— 2006. — №33. — С.31 — 38.
7. Рибін О.І. Аналіз лінійних систем з використанням кратних перетворень / О.І. Рибін, Ю.Х. Ніжебецька, І.О. Рибіна // Вісник НТУУ«КПІ». Сер. Радіотехніка. Радіоа-паратобудування.— 2010. — №40.— С. 5 — 11.
8. Рибіна І.О. Аналіз лінійних систем з використанням кратних перетворень шляхом розкладу реакції системи в ряд Тейлора / І.О. Рибіна, О.Я. Вівчарик, О.А. Якубенко // Вісник НТУУ «КПІ». Сер. Радіотехніка. Радіоапаратобудування. — 2011.— №44.— С.37—48.
9. Рибін О.І. Алгоритми побудови матриці деградації образу з використанням дискретного перетворення Адамара/ О.І. Рибін, Н.О. Іванюк // Вісник НТУУ «КПІ». Сер. Радіотехніка. Радіоапаратобудування. — 2012. — №49.— С. 18— 32.
Рибін О.І., Іванюк Н.О. Умовна деконволюція в області трансформант Фур’є. Побудова оберненої матриці деградації образу. В статті запропоновано алгоритм символьного формування матриці кратного перетворення для приведення матриці параметрів фільтра, який корегує спотворений образ (деградований за рахунок неточковості імпульсної характеристики вхідного пристрою), до діагонального вигляду. Це, в свою чергу, дозволяє легко знайти зворотну матрицю параметрів корегуючого фільтра, що забезпечує можливість реалізувати реставрацію образів за методом умовної деконволюції в базисі перетворення Адамара. Використання символьних залежностей при формуванні матриць з урахуванням їх жорсткої структури дозволяє значно скоротити кількість арифметичних операцій при програмній реалізації алгоритмів методу та підвищити точність обчислень. Внаслідок використання запропонованого методу обернення матриці корекції, метод стає конкурентоспроможним (за швидкодією) з методом умовної деконволюції в області перетворення Фур’є. На відміну від базису перетворення Фур’є розв’язання задачі умовної деконволюції в базисі перетворення Адамара для сигналів, які мають розриви першого роду (скоки), є більш прийнятним, оскільки самі функції Адамара є скоковими.
Ключові слова: реставрація образу, перетворення Адамара, умовна деконволюція, точність, символьні методи, обернення матриці.
Рыбин А.И., Иванюк Н.А. Условная деконволюция в области трансформант Фурье. Построение обратной матрицы деградации образа. В статье предложен алгоритм символьного формирования матрицы кратного преобразования для приведения матрицы параметров фильтра, который корректирует искаженный образ (деградировавший за счет неточечности импульсной характеристики входного устройства), к диагональному виду. Это, в свою очередь, позволяет легко найти обратную матрицу параметров корректирующего фильтра, что обеспечивает возможность реализовать реставрацию образов методом условной деконволюции в базисе преобразования Ада-мара. Использование символьных зависимостей при формировании матриц с учетом их жесткой структуры позволяет значительно сократить количество арифметических
28
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№50
Обчислювальні методи в радіоелектроніці
операций при программной реализации алгоритмов метода и повысить точность вычислений. В следствии использования предложенного метода обращения матрицы коррекции, метод становится конкурентоспособным (по быстродействию) с методом условной деконволюции в области преобразования Фурье. В отличие от базиса преобразования Фурье решения задачи условной деконволюции в базисе преобразования Ада-мара для сигналов, которые имеют разрывы первого рода (скачки), является более приемлемым, поскольку сами функции Адамара есть скачкообразными.
Ключевые слова: реставрация образа, преобразования Адамара, условная деконволюции, точность, символьные методы обращения матрицы.
Rybin A.I., Ivaniuk N.A. Conventional deconvolution in the Fourier transforms. The image degradation inverse matrix construction. Algorithm for symbolic formation of the multiple transformation matrix for the filter parameters matrix reduction to diagonal form is proposed in this paper. The distorted image is adjusted by this filter (degraded thru input device non-point impulse response). This leads it easy to find the correction filter parameters inverse matrix, which provides an opportunity to realize the image restoration by conditional deconvolution in the basis of Hadamard transform. Using the symbolic dependences in the matrices formation regarding to their tough structure allows significantly reducing the arithmetic operation number in the software implementation of method algorithms and to improve accuracy calculation. From the use of proposed correction matrix inversion method it becomes competitive (for speed) with the conventional deconvolution method in the Fourier transform. In contrast to the Fourier transform basis the solution of the conventional deconvolution problem on Hadamard transform basis for signals with first-order discontinuities (jumps) is more acceptable, since the Hadamard functions are jumplike themselves.
Keywords: image restoration, Hadamard transform, conventional deconvolution, accuracy, symbolic method of matrix inversion.
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№50
29
