Курсовая работа: Метод вейвлет-перетворення

АНОТАЦІЯ

В даній роботі представлений один із перспективних методів математичного аналізу – вейвлет-перетворення, застосування якого дозволяє оброблювати сигнали будь-якого виду (в даному випадку медико-біологічний, а саме – фотоплетизмограма). Проводиться порівняння з Фур′є-аналізом і обґрунтовано доведено переваги вейвлет-перетворення. Розроблено програмний комплекс для обробки фотоплетизмограми.


ВСТУП

 

Сьогодні в медичну діагностику впроваджується все більша кількість методів, основаних на застосуванні лазерних та оптико-електронних приладів. До них відноситься і фотоплетизмографічний метод (ФПМ), що дозволяє вимірювати кровонаповнення та кровострум як в потужних венах і артеріях, так і в периферійних судинах і капілярах.

ФПМ у порівнянні з іншими методами діагностики біологічного об'єкту (БО) за оптичними показниками, наприклад з фотоакустичним методом, дозволяє підвищити достовірність реєстрації гемодинамічних показників кровонаповнення, а також те, що введенням в прилади, які реалізують даний метод, елементів світловолоконної техніки і джерел з різноманітними довжинами хвиль зондуючого випромінювання можна достатньо точно вирішувати задачі фотодинамічних досліджень, дистанційних вимірів тих або інших гемодинамічних показників БО.

Розробка нових більш ефективних лазерних та оптико-електронних комп'ютеризованих систем та комплексів та методів диференціальної діагностики стоматологічних захворювань залишається однією із актуальних задач сьогодення.

Оптичний метод діагностики мікроциркуляції судин характеризується достатньо широким діапазоном можливостей реєстрації найрізноманітніших фізіологічних функцій тканин, органів і систем організму. Також відмінною рисою параметрів є їх висока вибірність і точність. Оптичний метод також дозволяє використовувати поряд з лазерними та оптико-електронними датчиками гнучкі скловолоконні світловоди для дослідження мікроциркуляції.

Даний метод дозволяє проводити комплексну оцінку мікроциркуляторного русла по двох важливих показниках: морфологічним ознакам і функціональним характеристикам. Комплексний аналіз дозволяє одержати досить повну інформацію про стан мікроциркуляторного русла в нормі і патології.

За допомогою оптичного методу дослідження визначають ряд функціональних показників, що властиві усередині судин (рівень кровонаповнення, швидкість і характер кровопотоку, тромбоутворення).

Широке розповсюдження серцево-судинних захворювань підтверджує актуальність розробок сучасних приладів діагностики та моніторингу, спрямованих на підвищення ефективності методів та розвиток технічних засобів діагностики таких захворювань.

За останні роки на основі досягнень медичної фізики сформувався новий напрямок - біоінженерія, основною задачею якої є розробка технічних систем і нових високоефективних технологій для діагностики, профілактики, лікування патологічних станів, і реабілітації. Біотехнізація сучасної медицини вимагає нової взаємодії між фізико-технічними і медико-біологічними науками. В багатьох країнах чітко проглядається тенденція до формування біонженерних (медико-технічних) центрів [1].

У даній роботі приводиться огляд одного із сучасних напрямків розвитку вейвлет-аналіза. Насамперед, актуальність даної роботи варто розглядати в контексті бурхливого розвитку вейвлет-аналіза й найширшого кола сфер його застосування. Так, уже зараз вейвлет-аналіз зарекомендував себе як один з ефективних методів кодування сигналів, обробки зображень будь-якої природи, супутникові зображення, рентгенограми внутрішніх органів, архівації даних, аналізу складних особливостей сигналів, об'єднання й поділи сигналів, створення множинного доступу, прихованого зв'язку, спільного кодування джерела й каналу зв'язку, виділення сигналів на фоні шумів, а також інтерес викликає його застосування й у сфері контролю якості передачі інформації.[2]

Взагалі, реально працюючі у додатках математичні методи завжди (чомусь) опираються на чисту математику - це експериментальний факт. А от прикладна сторона вейвлетів проста на стільки, що далі нікуди. При цьому вейвлет-перетворення не тільки працює швидко, але і його програмна реалізація незрівнянно проста[3].


1. ОСОБЛИВОСТІ ВЗАЄМОДІЇ ОПТИЧНОГО ТА ЛАЗЕРНОГО ВИПРОМІНЮВАННЯ З БІОЛОГІЧНИМИ СИСТЕМАМИ

 

Використання лазерів у біології та медицині може здійснюватися в кількох напрямках, одним з яких можна вважати розробку на основі лазерної техніки приладів та методів для виявлення, ідентифікації, дослідження будови біологічних об’єктів, а також для вивчення природи процесів, що відбуваються в них [4].Застосування лазерів у біології і медицині засновано на використанні широкого кола явищ, пов'язаних із різноманітними проявами взаємодії світла з біологічними об'єктами. Лазерне випромінювання, так само як і звичайне світло, може відбиватися, поглинатися, розсіюватися, перевипромінюватися біологічним середовищем, і кожний із цих процесів несе інформацію про мікро- і макроструктуру цього середовища, рух і форму окремих його складових. Червоне, інфрачервоне (ІЧ) та ультрафіолетове (УФ) світло можуть надавати фотобіохімічну дію. Яскравими прикладами цього є фотосинтез рослин і бактерій, а також механізм зору. Високоінтенсивне світлове випромінювання ультрафіолетового (УФ), видимого червоного та інфрачервоного (ІЧ) діапазонів довжин хвиль робить руйнівну (деструктивну) дію на біологічні об'єкти. Необхідні інтенсивності можна створити і не тільки за допомогою лазерів [3,4]. Таким чином, процеси, що характеризують види взаємодій лазерного випромінювання з біооб'єктами, можна розділити на три групи. До першої відносять усі неспотворювальні взаємодії (принаймні, у межах похибок вимірів, що не здійснюють помітної дії на біооб'єкт), до другого - процеси, у яких виявляється фотохімічна дія, і до третього - процеси, що призводять до фотодеструкції. На рисунку 1 подана класифікація основних принципів застосування лазерів у біології і медицині, що враховує зазначені групи процесів.

Оскільки ми маємо справу з живими об'єктами, то крім фізико-хімічних проявів дії лазерного випромінювання необхідно враховувати його вплив і на функціонування живої матерії. Цей вплив визначається ступенем гомеостазу живого об'єкта [5].

Ступінь гомеостазу характеризує стани і процеси, що забезпечують стабільність організму до зовнішніх втручань, вона є функцією еволюційного розвитку і виявляється найнижчою у біологічних молекул і найвищою в хребетних тварин.

Світло малої інтенсивності не запускає адаптаційні механізми біосистеми, з ростом інтенсивності спочатку це стосується гомеостазу живої системи на локальному рівні, потім включаються загальні адаптаційні і регуляційні механізми системи, що повністю її відновлюють, далі вони вже не справляються з повним відновленням і частково відбуваються необоротні процеси, що наростають і призводять до руйнацій у системі. Проте об'єкт можна ще вважати «живим». При високих інтенсивностях руйнації виявляються настільки значними, що об'єкт уже не може вважатися «живим» [5,6].

У дослідах по порівнянню поглинання червоного випромінювання з різними фізичними властивостями було встановлено, що просторова когерентність не впливає на поглинання, а поляризоване випромінювання поглинається менш активно ніж неполяризоване. Встановлено також, що розсіювання видимого світла при проходженні його через біотканину значно перевищує поглинання. Це означає, що лазерне світло має досить високу здатність проникнення в тканини. Якщо врахувати можливість транспортування випромінювання вглиб тканини при допомозі волоконної оптики і можливе наступне його розсіювання то можна сподіватися на подальше розширення сфери клінічного використання лазерів [6].

1. 1 Принципи розповсюдження оптичного та лазерного випромінювання в багатошарових тканинах

Вплив лазерного випромінювання на біологічний матеріал або реакція живої тканини на це випромінювання обумовлені взаємодією фотонів і молекул, або з'єднань молекул тканини. Атомарні і молекулярні процеси і наступні біологічні реакції вияснені ще не цілком. Відомі процеси можуть бути підрозділені на фотохімічну взаємодію, термічну взаємодію і нелінійні процеси.

Ступінь того або іншого впливу залежить:

а) від властивостей лазерного випромінювання (довжина хвилі, густота енергії, тривалість опромінення і частота повторення);

б) від властивостей біологічного матеріалу (коефіцієнт поглинання, коефіцієнт розсіювання, густота і т.д.).

У залежності від довжини хвилі, густоти енергії і часу впливу лазерного випромінювання ефект визначається в основному двома внутрішніми параметрами тканини: з одного боку, оптичними властивостями тканини, що опромінюється і, з іншого боку, її термічними властивостями.

При попаданні лазерного променя на тканину можуть спостерігатися три процеси: відбиття, поглинання і/або пропускання - тільки незначний відсоток випромінювання відбивається безпосередньо від поверхні (рисунок 1.1).

Рисунок 1.1 - Оптичні властивості прошарку матерії. Падаючий променевий потік розділяється на три частини: відбита частина Rф, поглинена частина Аф і пропущена частина Тф: Рф+Аф+Тф=1

Промені, що проникають в тканину, частково поглинаються, частково розсіюються і частково пропускаються (рисунок 1.2).

Рисунок 1.2 - Оптичні властивості лазерного променя на шкірі

В залежності від довжини хвилі випромінювання, що падає, відбивається до 60% випромінювання. Розсіювання залежить від негомогенних структур тканини і визначається різними показниками заломлення в різних шарах і різницею між шарами і їх навколишнім середовищем. Хвилі з довжиною набагато більшою, ніж діаметр шару (³ 10 мкм), розсіюються клітинними структурами лише в незначній мірі. Але тому, що електромагнітний спектр широко використовуваних лазерів простягається від ІЧ (1 мм-0,78 мкм) до УФ (0,38-0,10 мкм) діапазону довжин хвиль, ми практично завжди маємо справу з розсіюванням. Глибину проникнення для довжини хвиль більше 1,0 мкм можна розрахувати на основі закону Ламберта-Бера в першому наближенні [7].

Інтенсивність I випромінювання, що пройшло через прошарок товщиною d визначається співвідношенням:

I=I0e-ad,       (1.1)

 

де I0 - інтенсивність при вході в речовину і a - коефіцієнт поглинання.

При застосуванні монохроматичного випромінювання довжиною хвилі l для коефіцієнта поглинання дійсне таке співвідношення:


a = 4pnk/l, (1.2)

 

причому показники переломлення n і поглинання k є константами для даного середовища. Співвідношення Ламберта-Бера справедливе в тому випадку, коли поглинання набагато перевищує розсіювання [8].

Частіше всього пропонується рішення опису взаємодії лазерного випрмінювання з біотканинами з позицій теорії радіаційного переносу [9], при цьому бiотканина аналізується як випадково-неоднорідне середовище, яке розсіює та поглинає, а випромінювання, що розповсюджується в ній, – як потік енергії, тобто всі ефекти, зв'язані з хвильовою природою світла (дифракція, iнтерференція, поляризація ), не приймаються до уваги.

Основне рівняння теорії радіаційного переносу може бути записане в вигляді

, (1.3)

де I (z, q) – потужність випромінювання, що розповсюджується на глибині z через одиничний майданчик і в одиничному тілесному куті в напрямку, який складає з нормаллю до цього майданчика кут, конус якого рівний q, Вт×м-2×стер-1; mа тa ms – коефіцієнти поглинання і розсіювання, м-1; Р ((q', q) – фазова функція розсіювання, що описує вірогідність того, що світло розповсюджується в напрямку q.

Найкращим чином співвідношення поглинання і розсіювання описане в теорії Кубелки-Мунка [8,9]. Рівняння, що описує поширення випромінювання в середовищах з врахуванням поглинання і розсіювання має вигляд:

dLc(r,z)/dz = -gLc(r,z), (1.4)

де Lc(r,z) - щільність потужності випромінювання [Вт/м2] колімованого променя в місці р (вектор місця) у напрямку z, g - коефіцієнт ослаблення (сума коефіцієнтів розсіювання [м-1] і поглинання [м-1]).

Розсіювання в біологічній тканині залежить від довжини хвилі лазерного променя. Випромінювання ексимерного лазера УФ діапазону (193, 248, 308 і 351 мкм), а також ІЧ-випромінювання 2,9 мкм ErYAG-лазера і 10,6 мкм СО2-лазеру мають глибину проникнення від 1 до 20 мкм [10,11]. Тут розсіювання грає другорядну роль. Для світла з довжиною хвилі 450-590 нм, що відповідає лініям аргону, глибина проникнення складає в середньому 0,5-2,5мм. Як поглинання так і розсіювання грають тут значну роль. Лазерний промінь цієї довжини хвилі хоча і залишається в тканині колімованим у центрі, але він оточений зоною з високим розсіюванням. Від 15% до 40% енергії падаючого пучка світла розсіюється. У області спектра між 590 і 1500 нм, у яку входять лінії Nd:YAG лазера 1,06 і 1,32 мкм, домінує розсіювання. Глибина проникнення складає від 2,0 до 8,0 мм.

 

1.2 Аналіз оптико-електронних ІІС для аналізу гемодинамічних показників

Реанімаційно-хірургічні монітори ЮМ-300 мають вбудовану систему автоматизованого кардіо- і реоаналізу. В основу аналізу покладений метод математичної обробки плетизмограми і кардіоінтервалів, зареєстрованих протягом визначеного часу.

У результаті аналізу ритму серця будуються такі графіки:

·          РИТМОГРАМА - це послідовність вертикальних ліній, висота яких відповідає тривалості відповідного RR-інтервалу (в секундах). На осі абсцис відкладаються порядкові номери RR-інтервалів;

·          ГІСТОГРАМА (варіаційна пульсометрія) - східчаста функція розподілу RR-інтервалів у досліджуваному ряді їхніх значень;

·          СКАТТЕРГРАМА - послідовне нанесення на графік у прямокутній системі координат двох сусідніх RR-інтервалів. Скаттерграма особливо ефективна, при діагностиці аритмії.

·          Крім того, проводиться спектральний аналіз ритму серця, будується відповідний графік (спектр ритму) і визначаються частотні складові спектра.

Система кардіоаналізу дозволяє визначати наступні параметри:

ДХ =RRmax - RRmin - варіаційний розмах; Мо - мода (значення RR-інтервалу, що найчастіше зустрічається); АМо - амплітуда моди (число реалізації (у відсотках) даної моди стосовно загального числа аналізованих кардіоінтервалів); RR - середня тривалість (мс) синусових RR-інтервалів за 5 хв: SDNN - середньоквадратичне відхилення від середньої тривалості всіх синусових RR-інтервалів (за 5 хв.); Cv - коефіцієнт варіації (%) - відношення SDNN до RR; SDANN - середньоквадратичне відхилення від середніх тривалостей синусових інтервалів, розрахованих на всіх 5-хвилинних інтервалах запису; PNN50 - частка сусідніх синусових інтервалів, що розрізняються більш ніж на 50 мс; RMSSD - середньоквадратичне відхилення між тривалістю сусідніх синусових інтервалів; ULF - потужність (енергія) складової спектра ультранизької частоти (0-0,04 Гц); LF - потужність низькочастотної складової спектра (0,04-0,15 Гц); HF - потужність високочастотної складової спектра (0,15-0,4 Гц); LF/HF - відношення потужностей низькочастотної і високочастотної складових спектра.

З аналізу плетизмографічних даних визначаються наступні параметри: t1 - тривалість RR-інтервалу; t2 - анакротична фаза; t3 - катакротична фаза; t4 - період швидкого кровонаповнення; t5 - період повільного кровонаповнення; t6 - період венозного відтоку; Vsf - сфігмографічна швидкість; 11 - дикротичний індекс; 12 - діастолічний індекс; 13 - індекс периферійного опору; 14 - інтегральний гідравлічний індекс; 15 - інтегральний артеріальний індекс; 16 - інтегральний венозний індекс;

Всі накопичені дані можна передати в комп'ютер і роздрукувати за допомогою спеціального програмного забезпечення [10].


2 . ВИКОРИСТАННЯ ПЕРЕТВОРЕННЯ ФУР’Є ДЛЯ АНАЛІЗУ ПУЛЬСОВОЇ ХВИЛІ

Програма обчислення миттєвого значення частоти ударів пульсу входить у склад спеціалізованого програмного забезпечення (СПЗ) та проводить обробку аналогових сигналів, що надходять з датчиків пульсової хвилі різного типу (оптоелектронних, ємнісних, тензометричних і т.д.), з метою обчислення періоду пульсової хвилі та перерахунка миттєвого значення частоти ударів пульсу за хвилину. Задача ускладнюється тим, що пульсовій хвилі, як і іншим біомедичним сигналам, що повторюються, притаманний квазіперіодичний характер [5]. Це означає, що кожний наступний період сигналу лише приблизно відповідає попередньому, особливо за амплітудою відповідних ділянок (наприклад, екстремальних значень). Крім того, форма сигналу може різко змінюватись від періоду до періоду у зв'язку із загальним хвилюванням дослідженого хворого, що був поміщений у незвичну для нього обстановку. Тому неможливо використовувати прості методи обчислення періоду сигналу, що полягають у пошуку екстремальних точок з однаковою амплітудою. Хороші результати отримують при використанні перетворення Фур'є та аналізі періоду першої гармоніки розкладеної у спектр сигналу, однак дані методи потребують значних обчислювальних витрат за часом та об'ємом оперативної пам'яті.

В результаті моделювання запропонований достатньо простий швидкодіючий засіб обчислення миттєвого значення частоти пульсу, що використовує прості операції складання, віднімання та порівняння.Алгоритм обчислення миттєвого значення частоти ударів пульсу

Запропонований алгоритм обчислення миттєвого значення частоти ударів пульсу складається з таких етапів:

- виконується настройка таймера ТО ОМЕВМ на заданий дискрет часу t у нс роботи АЦП;

- виконується установка коефіцієнта перерахунка лічильника для коректування правильної роботи системних годинників;

- виконується установка максимального розміру місця у ОЗП, що відводиться під буфер відліків з АЦП N буфера;

- виконується установка програмного флага роботи АЦП, тобто дозволяється робота АЦП, що виконується програмою обробки переривання від таймера;

- відбувається циклічний аналіз стану програмного флага переповнення таймера ТО. Якщо підпрограма обробки переривання від таймера, то ще не заповнений весь буфер відліків у ОЗП, програмний флаг роботи зберігає одиничне значення. Таким чином, відбувається зупинка роботи основної програми до повного заповнення буферів відліку, коли підпрограма обробки переривання спрацьовує програмний флаг роботи та зупиняє запис відліків у буфер;

- відбувається перегляд всього буфера підрахунків та пошук максимального за величиною відліку Аmax, наповнення суми відліків та обчислення середнього за буфером відліку

A = sum [Aі, де i=0...Nбуф-1] / Nбуф

де і - номер відліку,

Аі - і-й відлік буфера,

Nбуф - розмір буфера відліків.

Якщо розмір буфера обліку Nбуф обраний кратним степеню двійки, тобто Nбуф = 2^j, де j=1,2,3..., то операція поділу накопичення суми замінюється на просте відкидання j молодших розрядів накопиченої суми відліків за буфером;

- відбувається обчислення величини порога Ф, який дозволяє виділити характерні фрагменти буфера відліку


Ф = A + (Amax - A)/2 = (Amax+A)/2 ;

- відбувається порогова обробка вихідного масиву відліку та занулення тих відліків, у яких величина менша порога, тобто порогових відліків Аі:

 

 Аі, якщо Аі>P

Аі' =

 0, якщо Аі<=P

Дана обробка дозволяє впустити лише ті фрагменти, які містять систолічний/діастолічний піки пульсової хвилі, а виключає невірні та випадкові викидає;

- відбувається пошук найближчого від початку буфера фрагмента A1 = {A1', де i=m...n} для якого А1'<>0 та у даному фрагменту визначається координата Хмах за величиною відліку А'(x) = max [A']. Якщо існує декілька розташованих поряд та рівних за величиною екстремальних відліків, тобто є у наявності горизонтальна "поличка" хвилі, відбувається обчислення ширини полички d та корекції координати Х на половину ширини полички Х1 = Х + d/2;

- координата Х1 записується як координата першої хвилі пульсу;

- відбувається обнуління знайденого фрагмента А1, тобто А1'=0 для всіх i=m...n;

- відбувається повторний пошук найближчого від початку буфера нульового фрагмента А11 = {Ai', i=k...1}, причому використовується одна й та ж підпрограма пошуку, що й у попередньому випадку;

-у фрагменту А11, як і раніше, визначається координата Х11максимального за величиною відліку з урахуванням можливої наявності горизонтальної полички. Координата Х11 запам'ятовується як координата наступної пульсової хвилі;

- вираховується період Т пульсової хвилі шляхом визначення різниці Х11-Х1 та множення його на тривалість одного такту роботи АЦП.

T = (X11 - X1) Xt;

- оскільки тривалість такту t вимірюється в одиницях мілісекунд, то миттєва частота ударів пульсу за хвилину F визначається як

F = 60 000 / T

Оскільки результат обчислень F представлений у загальному двійковому форматі, відбувається його перетворення у двійково-десяткову форму, зручну для людського сприйняття, та вивід результату обчислень на дисплей, тобто запис кодів Семи-сегментних індикаторів у буферний ЗП контролера клавіатури та дисплея.

Але під час виконання роботи був знайдений більш ефективний метод для аналізу пульсової хвилі – вейвлет-аналіз, якому і присвячений наступний розділ.


3. СУТНІСТЬ ВЕЙВЛЕТ-АНАЛІЗУ

 

Вейвлет-перетвореня сигналів є узагальненням спектрального аналізу, типовий представник якого - класичне перетворення Фур'є. Застосовувані для цієї мети базиси названі вейвлетами. Термін “вейвлет” пішов від англійського wavelet, що на українську мову переводиться як “коротка хвиля''. У математичній літературі поняття “вейвлет” позначають іноді словом “сплеск”, що звужує саме поняття, тим більше, що вейвлети й призначені для аналізу сплесків - сигналів нестаціонарного характеру.

Введені порівняно недавно, в 80-х роках, вони в наступні роки одержали швидкий теоретичний розвиток і широке застосування в різних областях обробки сигналів і зображень. На відміну від традиційного перетворення Фур'є, вейвлет-перетворення забезпечує двовимірне подання досліджуваного сигналу в частотній області в площині частота-положення. Аналогом частоти при цьому є масштаб аргументу базисної функції (найчастіше часу), а положення характеризується її зрушенням. Це дозволяє розділити великі й дрібні деталі сигналів, одночасно локалізуючи їх на тимчасовій шкалі. Іншими словами вейвлет-аналіз можна охарактеризувати як локалізований спектральний аналіз або - спектральний аналіз локальних збурювань. Апаратурним аналогом одного з видів вейвлет-аналіза є багато канальна смугова фільтрація сигналу при постійному відношенні ширини смуги фільтра до центральної частоти.

Вейвлет-аналіз розроблений для рішення завдань, які виявилися занадто складними для традиційного аналізу Фур'є. Перетворення Фур'є представляє сигнал, заданий у тимчасовій області, у вигляді розкладання по ортогональних базисних функціях (синусам і косинусам) з виділенням частотних компонентів. Недолік перетворення Фур'є полягає в тому, що частотні компоненти не можуть бути локалізовані в часі, його застосовують тільки в аналізі стаціонарних сигналів, у той час як багато сигналів мають складні частотно-часові характеристики. Як правило, такі сигнали складаються із близьких за часом, коротких високочастотних компонентів і довгих, близьких по частоті низькочастотних компонентів. Для аналізу таких сигналів необхідний метод, здатний забезпечити одночасний дозвіл як по частоті, так і за часом. Перше необхідно для локалізації низькочастотних складових, друге - для виділення компонентів високої частоти. Існує два підходи до аналізу нестаціонарних сигналів такого типу. Перший заснований на локальному перетворенні Фур'є. Прямуючи цим шляхом, нестаціонарний сигнал зводиться до стаціонарного шляхом його попереднього розбиття на сегменти (фрейми), статистика яких не змінюється з часом. Другий підхід полягає у використанні вейвлет-перетворення.

Всім відомо, що будь-який сигнал можна розкласти в суму гармонік (синусоїд) різної частоти. Але синусоїдальні хвилі нескінченні, і не дуже добре відслідковують зміни сигналу в часі. Щоб вловити ці зміни, замість нескінченних хвиль можна взяти зовсім однакові, але розподілені за часом короткі "сплески". Однак, як виявилося, цього недостатньо, треба додати ще їхні стислі копії. От тепер сигнал можна розкласти на суму таких сплесків різного розміру й місця розташування. Коефіцієнти розкладу, які несуть інформацію про еволюції сигналу, залежать від вибору початкового сплеску. Для кожного прикладного завдання можна підібрати найбільш пристосований (саме для неї) сплеск, що і називається вейвлетом. Математична сторона вейвлет-аналіза – річ досить тонка, хоча й достатньо наочна[11].


4. АНАЛІЗ ВЕЙВЛЕТ-ПЕРТВОРЕННЯ. ПОРІВНЯННЯ З ФУР’Є-АНАЛІЗОМ

Протягом багатьох десятиліть і по теперішній час основним засобом аналізу реальних фізичних процесів був гармонійний аналіз. Математичною основою аналізу є перетворення Фур'є. Перетворення Фур'є розкладає довільний процес на елементарні гармонійні коливання з різними частотами, а всі необхідні властивості й формули виражаються за допомогою однієї базисної функції exp(jwt) або двох дійсних функцій sin(wt) і cos(wt). Гармонійні коливання мають широке розповсюдження в природі, і тому зміст перетворення Фур'є інтуїтивно зрозумілий незалежно від математичної аналітики.

Перетворення Фур'є володіє рядом чудових властивостей. Оператор зворотного перетворення Фур'є збігається з вираженням для комплексно - сполученого оператора. Областю визначення перетворення є простір L2 інтегрувальних із квадратом функцій, і багато реальних фізичних процесів, спостережувані в природі, можна вважати функціями часу, що належать цьому простору. Для застосування перетворення розроблені ефективні обчислювальні процедури типу швидкого перетворення Фур'є (ШПФ). Ці процедури входять до складу всіх пакетів прикладних математичних програм і реалізовані апаратно в різних процесорах обробки сигналів.

Вейвлетне перетворення має багато спільного з перетворенням Фур'є. У той же час є ряд досить істотних відмінностей. Як приклад розглянемо застосування вейвлет-аналіза до синусоїд f(t)=sin(2πt/T1)+α sin(2πt/T2) , що дозволяє легко порівняти з результатами звичайного перетворення Фур'є.

На рисунку 4.1 показаний сигнал у вигляді суми синусоїд, що відрізняються частотами: (y=sin(30*x)+sin(100*x)).


Рисунок 4.1 - Сума синусоїд , що відрізняються частотами

Вейвлет-перетворення такого сигналу виявляє періодичну структуру не гірше й не краще перетворення Фур'є. На рисунку 4.2 видні дві широких смуги, що відповідають двом різним частотам.

Рисунок 4.2 - Вейвлет перетворення суми синусоїд з різними частотами

Однак відмінність цих двох спектральних аналізів проявляється, коли сигнал являє собою дві послідовні синусоїди з різними частотами ( рисунок 4.3).


Рисунок 4.3 - Дві послідовні в часі синусоїди з різними частотами

Як видно з рисунку 4.3 вейвлет-перетворення в цьому випадку дозволяє простежити еволюцію частоти сигналу в часі, тоді як Фур'є-спектр (рисунок 4.5) в обох випадках дасть нам тільки два піки й ніяк не відіб'є сам момент зміни частоти сигналу[12].

Рисунок 4.4 - Вейвлет-перетворення двох послідовних у часі синусоїд з різними частотами

Рисунок 4.5 - Спектр Фур'є двох послідовних у часі синусоїд з різними частотами

4.1 Перетворення Фур'є (ПФ)

В основі спектрального аналізу сигналів лежить інтегральне перетворення й ряди Фур'є. Нагадаємо деякі математичні визначення.

У просторі функцій, заданих на кінцевому інтервалі (0,T), норма, як найбільш загальна числова характеристика довільної функції s(t), по визначенню обчислюється як корінь квадратний зі скалярного добутку функції. У загальному випадку, для комплексних функцій, квадрат норми (енергія сигналу) відповідає виразу:

||s(t)||2 = ás(t), s(t)ñ = s(t)·s*(t) dt, (4.1.1)

де s*(t) – функція, комплексно сполучена з s(t).

Якщо норма функції має кінцеве значення (інтеграл сходиться), то говорять, що функція належить простору функцій L2[R], R=[0,T], інтегрувальних із квадратом (простір Гильберта), і, відповідно, має кінцеву енергію. У просторі Гильберта на основі сукупності ортогональних функцій з нульовим скалярним добутком

áv(t), w(t)ñ = v(t)·w*(t) dt = 0 (4.1.2)

завжди може бути, створена система ортонормованих "осей" (базис простору), при цьому будь-який сигнал, що належить цьому простору, може бути представлений у вигляді вагової суми простих складових, проекцій сигналу на ці "осі" - базисних векторів. Значення проекцій визначаються скалярними добутками сигналу з відповідними функціями базисних "осей".

Базис простору може бути утворений будь-якою ортогональною системою функцій. Найбільше застосування в спектральному аналізі одержала система комплексних експонентних функцій. Проекції сигналу на даний базис визначаються виразом:

Sn = (1/T) s(t) exp(-jn··t) dt, n Î (-∞, ∞), (4.1.3)

де =2/T – частотний аргумент векторів. При відомих виразах базисних функцій сигнал s(t) однозначно визначається сукупністю коефіцієнтів Sn і може бути абсолютно точно відновлений (реконструйований) по цих коефіцієнтах:

s(t) =Sn exp(jn·Dw·t). (4.1.4)

Рівняння (4.1.3) і (4.1.4) називають прямим і зворотним перетворенням Фур'є сигналу s(t). Таким чином, будь-яка функція гильбертова простору може бути представлена у вигляді комплексного ряду Фур'є (4.1.4), що називають спектральним представленням сигналу або його Фур'є-образом.

На практиці ряд Фур'є обмежується певною кількістю членів N. Обмеження числа членів ряду значенням N означає апроксимацію нескінченного сигналу N - мірною системою базисних функцій спектра сигналу з певною погрішністю залежно від фактичного спектра сигналу. Ряд Фур'є рівномірно сходиться до s(t) по нормі (4.1.1):

||s(t) -Sn exp(jnDwt)|| = 0. (4.1.5)

Таким чином, ряд Фур'є - це розкладання сигналу s(t) по базисі простору L2(0,T) ортонормированных гармонійних функцій exp(jnDwt) зі зміною частоти, кратним частоті першої гармоніки w1=Dw.. Звідси, ортонормований базис простору L2(0,T) побудований з однієї функції v(t) = exp(jDwt) = cos(Dwt)+j·sin(Dwt) за допомогою масштабного перетворення незалежної змінної так, що vn(t) = v(nt).

Для коефіцієнтів ряду Фур'є справедлива рівність Парсеваля збереження енергії сигналу в різних представленнях:

 (1/T) |s(t)|2 dt = |Sn|2. (4.1.6)

Розклад в ряд Фур'є довільної функції y(t) коректно, якщо функція y(t) належить цьому ж простору L2(0,T), тобто квадратично інтегрувальна з кінцевою енергією:

|y(t)|2 dt < ¥ , t Î (0,T), (4.1. 7)

при цьому вона може бути періодично розширена й визначена на всій тимчасовій осі простору R(-¥, ¥) так, що

y(t) = y(t-T), t Î R,

за умови збереження кінцівки енергії в просторі R(-¥, ¥).

З позицій аналізу довільних сигналів і функцій у частотній області й точному відновленні після перетворень можна відзначити ряд недоліків розкладання сигналів у ряди Фур'є, які привели до появи віконного перетворення Фур'є й стимулювали розвиток вейвлетного перетворення. Відзначимо основні з них:

·          Обмежена інформативність аналізу нестаціонарних сигналів і практично повна відсутність можливостей аналізу їхніх особливостей (сингулярностей), тому що в частотній області відбувається «розмазування» особливостей сигналів (розривів, сходів, піків і т.п.) по всьому частотному діапазоні спектра.

·          Гармонійні базисні функції розкладу не здатні в принципі відображати перепади сигналів з нескінченною крутістю типу прямокутних імпульсів, тому що для цього потрібно нескінченно велика кількість членів ряду. При апроксимації стрибків нелокалізованими в часі базисними функціями необхідно, щоб суперпозиція цих функцій не тільки відновила стрибок, але й знищила один одного за межами стрибка, що робить рівнозначними всі компоненти його спектра. При обмеженні числа членів ряду Фур'є на околицях стрибків і розривів відновленого сигналу виникають осцилляции (явище Гіббса).

·          Перетворенням Фур'є відображаються глобальні відомості про частоти досліджуваного сигналу, оскільки базисні функції перетворення визначені на нескінченному тимчасовому інтервалі. ПФ не дає представлення про локальні властивості сигналу при швидких тимчасових змінах його спектрального складу. Так, наприклад, перетворення Фур'є не розрізняє сигнал із сумою двох синусоїд. Перетворення Фур'є в принципі не має можливості аналізувати частотні характеристики сигналу в довільні моменти часу.

4.2 Віконне перетворення Фур'є

Частковим виходом із цієї ситуації є так зване віконне перетворення Фур'є з віконною функцією, що рухається по сигналі, що має компактний носій. Повний часовий інтервал сигналу, особливо при великій його тривалості, розділяється на підінтервали, і перетворення Фур'є виконується послідовно для кожного вікна окремо. Тим самим здійснюється перехід до частотно-тимчасового (частотно-координатному) поданню сигналів і результатом віконного перетворення є сімейство спектрів, яким відображається зміна спектра сигналу по інтервалах зрушення вікна перетворення. Це якоюсь мірою дозволяє виділяти на координатній осі й аналізувати особливості нестаціонарних сигналів. Розмір носія віконної функції w(t) звичайно встановлюється порівнянним з інтервалом стаціонарності сигналу. Власне кажучи, таким перетворенням один нелокалізований базис розбивається на певну кількість базисів, локалізованих у межах функції w(t), що дозволяє представляти результат перетворення у вигляді функції двох змінних - частоти й тимчасового положення вікна. Віконне перетворення виконується відповідно до виразу:

S(w,bk) = s(t) w(t-bk) exp(-jwt) dt. (4.2.8)

Функція w(t-b) являє собою функцію вікна зрушення перетворення по координаті t, де параметром b задаються фіксовані значення зрушення. При зрушенні вікон з рівномірним кроком bk = kDb. В якості вікна перетворення може використовуватися як найпростіше прямокутне вікно ( w(t)=1 у межах вікна й 0 за його границями), так і спеціальні вагові вікна (Бартлетта, Гаусса, Кайзера та ін.), що забезпечують малі перекручування спектра за рахунок граничних умов вирізки віконних відрізків сигналів і нейтралізуюче явище Гіббса.

4.3 Приклад віконного перетворення

 

Приклад віконного перетворення для нестаціонарних сигналів на великому рівні шуму наведений на рисунку , наведеного у додатку А. По спектрі сигналу в цілому можна судити про наявність у його складі гармонійних коливань на трьох частотах. Віконне перетворення не тільки підтверджує даний висновок, але й показує конкретну локальність коливань по інтервалі сигналу й співвідношення між амплітудами цих коливань.

Координатна розв'язна здатність віконних перетворень визначається шириною віконної функції й, у силу принципу невизначеності Гейзенберга, обернено пропорційна частотній розв'язній здатності. При ширині віконної функції, рівної b, частотна розв'язна здатність визначається значенням Dw = 2p/b. При необхідній величині частотного дозволу Dw відповідно ширина віконної функції повинна, бути дорівнює b = 2p/Dw. Для віконних перетворень Фур'є ці обмеження є принциповими. При розмірі масиву даних N = 300 і ширині віконної функції Db = 100 частотна розв'язна здатність результатів перетворення зменшується в N/Db = 3 рази в порівнянні з вихідними даними, і графіки Sw(nDwSw) по координаті n для наочного зіставлення із графіком S(nDwS ) побудовано із кроком по частоті DwSw = 3DwS, тобто по точках n = 0, 3, 6, … , N...

4.4 Частотно-часові віконні перетворення

Функція віконних перетворень (4.2.8) може бути, переведена в тривимірний варіант із незалежними змінними й за часом, і по частоті:

S(t,w) = s(t-t) w(t) exp(-jwt) dt. (4.4.9)

На рисунку, наведеного у додатку Б, наведений приклад обчислення й представлення (модуль правої частини головного діапазону спектра) результатів тривимірної спектрограми при дискретному задані вхідного сигналу sq(n). Сигнал являє собою суму трьох послідовних радіоімпульсів з різними частотами без пауз, з відношенням сигнал/шум, близьким до 1. Віконна функція wi задана в однобічному варіанті з ефективною шириною вікна b @ 34 і повним розміром М =50. Установлений для результатів крок по частоті Dw = 0.1 трохи вище фактичної розв'язної здатності 2p/M = 0.126.

Для забезпечення роботи віконної функції по всьому інтервалі сигналу задавалися початкові й кінцеві умови обчислень (продовження на M крапок обох кінців сигналу нульовими значеннями).

Як видно за результатами обчислень, віконне перетворення дозволяє досить точно локалізувати інформативні особливості сигналу за часом і по частоті[13].

Використання дискретного вейвлет-перетворення дозволяє провести доведення багатьох положень теорії вейвлетів, пов'язаних з повнотою й ортогональністю базису, збіжністю рядів і т.д. Доказовість цих положень необхідна, наприклад, при стиску інформації або в завданнях чисельного моделювання, тобто у випадках, коли важливо провести розклад з мінімальним числом незалежних коефіцієнтів вейвлет-перетворення й мати точну формулу зворотного перетворення. Використання безперервного вейвлет-перетворення для аналізу сигналів більш зручно, а його деяка надмірність, пов'язана з безперервною зміною масштабного коефіцієнта а й параметра зрушення b, стає тут позитивною якістю, тому що дозволяє більш повно й чітко представити й проаналізувати інформацію, що міститься у вихідних даних. Зокрема, стає можливим проведення локалізації й класифікації особливих крапок і обчислення різних фрактальних характеристик сигналу, а також виконання частотно-часового аналізу нестаціонарних сигналів. Наприклад, у таких сигналів, як мовний сигнал, спектр радикально міняється в часі, а характер цих змін являє собою дуже важливу інформацію при розпізнаванні мови.

На основі вейвлетів створюються й такі елементи, як високочастотний і низькочастотний вейвлет-фільтри, за допомогою яких відбувається фільтрація сигналу по алгоритму Малла (рисунок 4.6). При цьому для збільшення дозволу вейвлет-фільтрів по частоті використається простий і досить ефективний прийом. Опишемо його для ортогонального випадку[2].


Рисунок 4.6 – Розклад по вейвлет-пакетам.

Сімейства вейвлетів у тимчасовій або частотній області використаються для представлення сигналів і функцій у вигляді суперпозицій вейвлетів на різних масштабних рівнях декомпозиції (розкладання) сигналів. Перші теоретичні роботи з основ вейвлетних перетворень були виконані в 90-х роках минулого століття Мейером (Mayer Y.), Добеши (Daubechies I.) і Маллатом (Mallat S.A.). Математичний апарат вейвлет-перетворення перебуває в стадії активної розробки, однак спеціальні пакети розширень по вейвлетам уже існують в основних системах комп'ютерної математики (Matlab, Mathematica, Mathcad, і ін.).

У цей час вейвлет-перетворення й вейвлетний аналіз використовуються в багатьох галузях науки й техніки для всяких завдань: для розпізнавання образів, для чисельного моделювання динаміки складних нелінійних процесів, для аналізу апаратної інформації й зображень у медицині, космічній техніці, астрономії, геофізиці, для ефективного стиску сигналів і передачі інформації з каналів з обмеженою пропускною здатністю й т.д.

4.5 Розклад по піддіапазонам

 

Іноді буває корисно розкласти сигнал на компоненти, енергія яких зосереджена в різних частотних піддіпазонах (тобто істотно відмінна від нуля на різних під відрізках відрізка ), і кодувати їх з різним ступенем детальності (наприклад, залежно від чутливості людського вуха до звуків різної частоти). Розподіл «енергії» сигналу по частотах характеризує , Задовго до створення вейвлет-аналіза для цього використалася схема, що ми зараз опишемо.

Ми хочемо знайти два фільтри,  (придушуючий високі частоти) і ( придушуючий низькі частоти), які дозволяли б розкласти сигнал на два компоненти,  і , удвічі їх прорідити (половина значень стає зайвою – адже частотний діапазон скоротився вдвічі!), а потім, за допомогою транспонованих фільтрів, точно відновити за цими даними вихідний сигнал (цю операцію можна застосовувати рекурсивно). Умови на шукані фільтри зручно записати в термінах z-перетворення.

Нехай  – z-перетворення однієї з компонентів. Перед кодуванням вона проріджується вдвічі, а перед відновленням вихідного сигналу доводить до вихідної довжини вставкою нулів між сусідніми значеннями. При цьому z-перетворення з  перетворюється в . Підставивши дане рівняння для кожного з фільтрів, одержимо z-перетворення компонентів перед відновленням

 (4.5.10)

z-перетворення транспонованих фільтрів мають вигляд  і . Сигнал відновиться з їхньою допомогою точно, якщо:

.

Одержуємо умови точного відновлення :

 (4.5.11)


У матричній формі вони записуються так:

,

де

 (4.5.12)

Підставивши , одержимо умови на ДПФ шуканих фільтрів:

    (4.5.13)

 

Допустимо, що ми знайшли  такий, що

 (4.5.14)

Тоді, підставивши

 (4.5.15)

ми бачимо, що умова виконується. Завдання звелося до знаходження тригонометричного багаточлена , що задовольняє умові. На методах побудови таких багаточленів ми зупинимося в наступній лекції. Фільтри  і , що задовольняють умові, називаються квадратурними дзеркальними фільтрами. На рисунку 4.7 (a) і (б), показані ДПФ такої пари фільтрів  і , а також вихідний сигнал до й після фільтрації (без проріджування)[12].



Рисунок 4.7(а) – Сигнал до фільтрації


Рисунок 4.7 (б) – Сигнал після фільтрації


5. ЗАСТОСУВАННЯ ВЕЙВЛЕТ-АНАЛІЗА ДЛЯ ОБРОБКИ СИГНАЛІВ

 

5.1 Огляд існуючих методів

 

5.1.1 Пірамідне представлення сигналів

На рисунку 5.1 схематично зображене пірамідне представлення одномірного сигналу. Сигналові ставляться у відповідність дві піраміди: піраміда гауссіанів (ПГ) і піраміда лапласіанів (ПЛ). Ці назви відбивають аналогію з популярними в графіку операціями згладжування (згортки з колоколообразним фільтром) і виділення перепадів (обчислення “дискретного оператора Лапласа”). Можна вважати цю конструкцію спрощеним варіантом попередньої.

В основі ПГ знаходиться вихідний сигнал. Наступний поверх ПГ – вихідний сигнал, профільтрований низькочастотним фільтром  і проріджений після цього вдвічі – передбачається, що фільтр h «убиває» верхню половину частотного діапазону, тому густоту вибірки можна відповідно зменшити. До цього поверху застосовується та ж операція, і так далі. У випадку кінцевих сигналів кожний наступний поверх удвічі коротше попереднього.

Рисунок 5.1 – Пірамідне представлення сигналів

Поверхи ПЛ – різниці між послідовними поверхами ПГ. Вони обчислюються так. Нехай, наприклад,  і  – перший і другий поверхи ПГ,  – перший поверх ПЛ, що ми хочемо обчислити. Для цього спочатку вирівнюються довжини поверхів:

а потім виконується фільтрація сполученим фільтром  (його коефіцієнти – переставлені у зворотному порядку коефіцієнти , у Фур'є-області це рівнозначно переходу к) . У результаті виникає вектор  . По визначенню, .Тепер замість вихідного сигналу ( ) досить запам'ятати пари ( ). Вихідний сигнал можна точно відновити по формулі:

Сигнал  удвічі коротше вихідного, а сигнал , як правило, майже цілком складається з дуже малих величин. Багато хто із цих величин можна без помітного збитку для точності відновлення замінити нулями, а інші закодувати більш короткими словами, чим компоненти вихідного сигналу. За рахунок цього загальна довжина запису ( ) буде істотно меншої довжини запису вихідного сигналу. Це скорочення стане ще більшим, якщо обчислити кілька поверхів ПЛ, і запам'ятовувати замість вихідного сигналу кілька поверхів ПЛ і останній поверх ПГ.

Ступінь стиску інформації цим методом залежить від вибору фільтра . При експериментах з пірамідними представленням було зроблене спостереження: «якість» фільтра зручно виражати в термінах еквівалентної вагової функції. Ця функція виникає так. Неважко обчислити коефіцієнти фільтрів, згортка сигналу з якими дає відразу другий поверх ПГ, третій поверх, і т.д. Виявляється, що при відповідній нормуванню вектори цих коефіцієнтів сходяться до якоїсь граничної «форми» – графікові функції , що повинні задовольняти функціональному рівнянню[12]:

 (3.1)


Процес одержання  зображений на рисунку 5.2.

Рисунок 5.2 – Процес одержання графікові функції  

 

5.1.2 Напівортогональний багато масштабний аналіз

Вейвлет-базис називається напівортогональним, якщо для будь-якого рівня дозволу  простір вейвлетів  ортогональний простору  (і, отже, всім просторам , , ... )[2]. Очевидно, що під класом напівортогональних вейвлетів є клас ортогональних вейвлетів, для якого додатково потрібна ортогональність базисних функцій . Відсутність такого обмеження дозволяє будувати, наприклад, гладкі симетричні вейвлети з компактним носієм (помітимо, що єдиними ортогональними симетричними вейвлетами з компактним носієм є вейвлети Хаара, які не володіють навіть безперервністю). У матричній формі умова напівортогональності можна записати в такий спосіб:


 

Якщо замість індексу j записати, то маючи
й , умова напівортогональності буде виглядати так:

.

Якщо  й  задані, то  є рішенням однорідної системи рівнянь , де  — відома матриця. Якщо однорідна система має нетривіальні рішення, то їх нескінченно багато, тобто  визначається неоднозначно. Тому для визначеності на накладається ряд додаткових умов. Наприклад, ми хочемо, щоб побудовані нами вейвлети мали компактний носій і були симетричні. Це значить, що стовпці матриці  повинні мати найменш можливе число підряд ідучих ненульових елементів, причому самі ланцюжки ненульових елементів повинні бути симетричними.

Прикладом напівортогональних вейвлетів є сплайнові вейвлети Сплайнові вейвлети будуються на основі B-сплайнів [3]. Існують різні види сплайнових вейвлетів. Ми розглянемо вейвлети, побудовані на основі нерівномірних B-сплайнів, що інтерполюють кінцеві крапки. Далі для стислості такі сплайни будемо називати просто B-сплайнами, а відповідні вейвлети — B-сплайновими вейвлетами. Будемо будувати B-сплайнові вейвлети на одиничному відрізку. Нехай m — ступінь сплайна, j — рівень дозволу. Простір  породжується  B-сплайнами, побудованими на послідовності вузлів


Неважко показати, що побудовані в такий спосіб простори ,  вкладені один в одного й задовольняють всім вимогам багато масштабного аналізу. На рис. 1 показані набори кубічних ( ) B-сплайнових скейлинг-функцій просторів  і . Матриця  має  стовпців і  рядків, всі стовпці, за винятком m перших і m останніх є зсуненими копіями стовпця , причому ненульові елементи цих стовпців є біноміальними коефіцієнтами, помноженими на . Нижче приводяться матриці , , і  для кубічного випадку. [11, 12].

  .

 

Рисунок 5.3 – B-сплайнові скейлинг-функції просторів  і .

Рисунок 5.4 – B-сплайнові вейвлети просторів  і .

Скейлинг-функції  й матриці  задані. Взято стандартний скалярний добуток в Тепер можна шукати матрицю . Помітимо, що ця матриця повинна мати  стовпців (розмірність простору ) і  рядків. Як було відзначено вище, матриця  (і, отже, вейвлет-базис) визначається неоднозначно. Матриця  побудована таким чином, щоб вона була розрідженою й містила мінімальне число підряд ідучих ненульових елементів у стовпцях. Структура такої матриці схожа на структуру матриці : вона розріджена і її стовпці крім m перших і m останніх є зсунененими копіями один відносно одного. Нижче приводяться матриці ,  і  для кубічного випадку, на рисунку 5.4 показані вейвлети просторів  і .

При  ми одержимо квадратичні B-сплайнові вейвлети, при  — лінійні, а при  — уже добре відомі ортогональні вейвлети Хаара[14].

Редагування кривої: Редагування здійснюється в такий спосіб: виконується декомпозиція вихідної кривої, отримані в результаті цього коефіцієнти деяким чином змінюються, після чого виробляється відновлення, але вже по модифікованому наборі коефіцієнтів. Можливі два принципово різних підходи: змінювати низькочастотну частину перетворення або змінювати високочастотну частину. У першому випадку можна міняти форму кривій “у цілому”, зберігаючи її дрібні особливості (рисунок 5.3), у другому - навпаки - зберігаючи форму, міняти деталі (рисунок 5.4). Очевидно, що при редагуванні кривих активно використається властивість локалізації вейвлетов у просторі, що дає можливість робити маніпуляції з окремими частинами кривої.


Рисунок 5.5 – Згладжування кривої

Рисунок 5.6 – Редагування кривої:


Рисунок 5.7 – Редагування кривої


6. ТЕХНІЧНІ ДАНІ

·          Діапазон вхідних напруг каналу реєстрації ЕКГ - 0.03.. 5 мВ

·          Чутливість каналу ЕКГ встановлюється з ряду - 2.5, 5, 10, 20, 40, 80 мм/мВ

·          Відхилення встановленої чутливості від номінальної - не більш ±10 %

·          Вхідний імпеданс каналу ЕКГ - не менше 15 МОм

·          Коефіцієнт ослаблення синфазних сигналів каналу ЕКГ - не менше 28000

·          Напруга внутрішніх шумів каналу ЕКГ - не більше 20 мкВ

·          Відхилення величини відображуваного на екрані каліброваного імпульсу (1 мВ) від номінальної - не більш ±10%

·          Нерівномірність АЧХ каналу ЕКГ в діапазоні 0.5.. 60Гц - 85..105%

·          Коефіцієнт придушення фільтру каналу ЕКГ на частоті 25 Гц - не більше 3 дВ не менш 2,5 с.

·          Постійна часу каналу фотоплетизмограми - не менше 1 с

·          Тривалість фронту фотоплетизмограми - не більше 0,1 с

·          Швидкість розгортки - 25, 50 мм/с

·          Постійна часу каналу ЕКГ Відхилення встановленої швидкості розгортки від номінальної - не більше ±10 %

·          Діапазон визначення SaО2 - 0 ..100 %

·          Відхилення значення SaО2 в діапазоні 80 . . 99 - не більше ± 2%, в діапазоні 50 . . 79 - не більше ± 4%, у діапазоні 0 . . 49 - не нормується

·          Діапазон установки значень порогу сигналізації по SaО2 - 50 .. 95 (із кроком 5)

·          Діапазон визначення ЧСС - 30 . . 250 уд/хв

·          Відхилення значення ЧСС від фактичного значення в діапазоні 30 . . 99 - не більш ± 2 уд/хв, в діапазоні 100 . . 250 - не більше ± 3 уд/хв, -

·          Діапазон установки значень порогів сигналізації по ЧСС (із кроком 10)

·          Живлення приладу від мережі змінного струму, напруга живлення - 220 ± 22 В, частота - 50 ± 0.5 Гц, споживана потужність - не більш 12 ВА

·          Габаритні розміри приладу - 290x213x119 мм

·          Габаритні розміри первинного перетворювача каналу SaО2 - 71x25x25 мм

·          Довжина кабелю первинного перетворювача каналу SaО2 - не менше 2 м

·          Довжина кабелю відведень - не менше 2.5 м

·          Маса приладу - не більше 2.5 кг

Прилад дозволяє автоматично здійснювати побудова гістограми розподілу КІ з кроком 8 мс, обсягом вибірки від 20 до 150 КІ [15].


7. ПАРАМЕТРИ Й ОБРОБЛЮВАНА ІНФОРМАЦІЯ

 

Оброблювальною інформацією являються: електрокардіограма, артеріальний тиск, фотоплетизмограма, реоенцефалограма, реовазограма.

Особливості системи:

·          можливість обробки й аналізу розширеного обсягу інформації від одного або декількох пацієнтів; обсяг додатково необхідної інформації може уточнюватися на етапі укладання договору на постачання;

·          наявність розширеної бази даних, що включає, крім інформації про пацієнтів і лікувальні сеанси, довідково-інформаційні і методичні матеріали по ГБО-терапії обсягом близько 200 найменувань;

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

Система Б-001.5 дозволяє включати в свій склад пристрої і програмне забезпечення "Мультимедіа" і здійснювати під час проведення лікувального сеансу голосовий супровід оброблюваних параметрів пацієнта і середовища.

Комплекс забезпечує наступні функціональні можливості:

Багатоканальну реєстрацію і графічне представлення фізіологічних показників у реальному масштабі часу з автоматичною установкою посилення і калібруванням сигналів, що вводяться, у діапазоні можливих змін [16].

Представлення в додатковому графічному вікні динаміки повільно протікаючих змін параметрів, що реєструються в реальному масштабі часу у вигляді середніх складових грудного і діафрагмального дихання, величини периферійного тиску, провідності шкірних покривів.

Запис у пам'ять комп'ютера даних, що вводяться, для наступного перегляду, редагування, аналізу і формування бази даних.

Введення з клавіатури оператором спеціальних маркерних оцінок у ході обстеження.

Аналіз даних тестування і подання результатів у вигляді гістограми виразності змін фізіологічних показників, гістограми імовірності і значимості відхилень фізіологічних функцій, представлення їхніх кількісних значень у табличному й образно-графічному вигляді (паттерни реакцій в полярних координатах).


8. РЕАЛІЗАЦІЯ ПРОГРАМНОГО ЗАБЕЗПЕЧЕННЯ ДЛЯ ОБРОБКИ ФОТОПЛЕТИЗМОГРАФІЧНИХ ДАНИХ

8.1 робота з базою даних

 

Програма роботи з базою даних представляє функції для створення картки пацієнта, пошуку, збереження ПГ й огляду інформації в базі.

База даних представляє собою список файлів. Використання довгих імен файлів в ОС Windows дозволяє називати файли прізвищами пацієнтів, що створює одразу зручну и гнучку систему збереження даних.

Реєстрація пацієнта:

Для реєстрації пацієнта необхідно вибрати пункт меню "Файл-Новий файл" або натиснути інструментальну кнопку (IК) .

Після цих дій на екрані з’явиться діалогове вікно зображене на рисунку 8.1.

Рисунок 8.1 – Вибір опції НОВИЙ ФАЙЛ

Кожний пункт реєстраційної карти при встановленому на ньому маркері виділяється [32].

Заповніть пункт "Прізвище, ім‘я", виберіть мишкою або за допомогою клавіатури стать, дату народження шаблон.

Поле шаблона вміщує список доступних шаблонів структур. Так в даній версії програми є такі версії як osteo.shb (вміщує структуру аналогічну структурі хребта) й dant.shb (стоматологія).

Встановіть маркер на полі "Затвердити". Після заповнення всіх пунктів процес реєстрації й створення нового файла закінчено. Корекція введених даних можлива в ході роботи програми в пункті меню "Вікна-Персона" або натискання IК

Рисунок 8.2 – Карта розміщення хребців

Після реєстрації з’явиться вікно "Навігатор" в якому активні елементи розміщенні по групам. Так навігатор для вертебродіагностики представлений на рисунку 8.2. Можливим є створення структур для навігатора довільної форми (різне число секцій й підсекцій) [17, 18].

Переміщення по навігатору відбувається з допомогою клавіатури або маніпулятора "миша". Після вибору певного елемента 3-го рівня навігатора ("До сеансу" або "Після сеансу") можна приступати до зняття ФПГ даних.

Реєстрація ФПГ:

Після вибору певного елемента навігатора (як, наприклад зображено на рисунку 8.3 ви можете здійснити зняття даних, про що свідчить активний пункт меню "Процес-Моніторинг" та IК .

При знятті ПХ за допомогою вбудованих функцій програми (про що свідчить вимкнений прапорець "Використовувати зовнішній пакет" у вікні настроювань ) необхідно вказати необхідні активні і видимі канали у вікні каналів . Якщо використовується зовнішній модуль по зняттю даних (як на дійсний момент), то вибір каналів не потрібний. Число каналів для цього випадку визначається зміною в конфігураційному файлі зовнішнього модуля.

Рисунок 8.3 – Карта розташування хребців з вибраним елементом

Рисунок 8.4 – Вікно каналів

Регулюючи силу притиснення датчиків у міжхребетних западинах і вибираючи коефіцієнт підсилення датчиків (для зовнішнього знімного пристрою це регулятор на передній панелі), досягнути появи стійких ПХ із максимальним рівнем амплітуди сигналу в двох каналах [19].

При використанні убудованих функцій зняття даних натисніть IК . При цьому відбудеться перехід від простого моніторингу до запам'ятовування даних, що знімається.

З цього моменту дані почнуть записуватися в пам'ять (у випадку зовнішнього модуля дані запам'ятовуються із самого початку). При цьому колір кривої зміниться. Запис виробляється до моменту натискання IК . Знята ФПГ автоматично записується в пам'ять для обраного елемента навігатора.

Для реєстрації наступної ФПГ виберіть новий елемент у навігаторі і повторіть описані дії.

Якщо Ви хочете перезаписати вже зняту ФПГ, виберіть той же елемент.

Подальші дії аналогічні тим, що виконувалися для першої ФПГ.

По закінченні зняття даних необхідно занести зняті ФПГ у базу даних. Для цього в меню виберіть процедури "Файл - Зберегти файл". Відкриється стандартне вікно, як показано на рисунку 8.5.

Рисунок 8.5 – Запис знятих ФПГ у базу даних


Перегляд ФПГ:

Зняті ФПГ зберігаються у виді, представленому на рисунку 8.5. Для перегляду значень ФПГ використовується інформаційна панель, що знаходиться під графіками і містить інформацію про координату, значення і похідні ФПГ у місці, де знаходиться курсор маніпулятора "миша".

Передбачено функцію масштабування що переглядаються ФПГ:  - збільшити,  - зменшити. Відображення похідних контролюється IК  (вихідний графік, 1-я і 2-я похідні).

Пересування по ФПГ здійснюється за допомогою смуги прокручування, розташованої над графіками.

Розміщення маркерів:

Для розрахунку основних параметрів пульсової хвилі, необхідно розставити маркери в характерних точках ФПГ як показано на рисунок 8.6. Для цього виберіть пункт "Обробка - Розставити маркери".

На кожній ПХ установлюється 5 маркерів:

1. Початок анакорти

2. Рівень швидкого кровонаповнення.

3. Рівень максимального кровонаповнення.

4. Рівень інцизури.

5. Венозний відтік.

Рисунок 8.6 – Розміщення маркерів

Після автоматичного розміщення маркерів їхнє розташування може бути відкоректоване.

Для точної установки маркерів необхідно, натиснувши ліву кнопку миші, пересунути курсор, стежачи за чисельними значеннями похідних у нижній інформаційній панелі.

Обробка фотоплетизмограм:

Даний режим призначений для виконання наступних функцій: програмної фільтрації перешкод, обчислення параметрів ПХ.

Для фільтрації низьких частот використовується метод параболічного наближення по n точкам (число точок для інтерполяції вказується у вікні настроювань  "Опції - Загальні - Кількість точок при фільтрації").

Діалогове вікно "Опції":

Всі опції розбиті на 5 класів (рисунок 8.7):

"Загальні" - опції, використовувані при обчисленнях і при збереженні даних:

Рисунок 8.7 – Діалогове вікно "Опції" (загальні)


"Нечутливість до перепадів" – параметр, що впливає на якість розміщення маркерів, "Архівувати дані при збереженні" - якщо параметр включений, те дані записуються у файл даних у стиснутому LZH-компресією вигляді, що приводить до збереження вільного місця на носії інформації, "Кількість точок при фільтрації" - фільтрація полягає в параболічній апроксимації (нерекурсивний фільтр низьких частот), тому це число (5,7 чи 9 точок) указує скільки точок буде впливати на значення в даній точці. 9 точок дає максимальне згладжування, "Схема розрахунку" - визначає вигляд звіту - як при діагностиці остеохондроза (дані до і після сеансу) чи щодо однієї еталонної хвилі, "Моніторинг" - визначає параметри інтерфейсу з апаратною частиною комплексу (рисунок 8.8).

Рисунок 8.8 – Діалогове вікно "Опції" (моніторинг)

Перший список, що випадає, відповідає за вибір джерела даних: плата АЦП – дані надходять із плати АЦП, використаної в комплексі для діагностики остеохондрозу, зовнішній пакет – для знімання даних запускається зовнішній пакет, результати роботи якого зчитуються потім з його зовнішнього тимчасового файлу в програму, COM-порт – взаємодія з останнім розробленим модулем, що підключається до COM2, калібрувальний сигнал – дані не надходять ні з якого пристрою. Синусоїда генерується програмно. Використовується для оцінки вірогідності даних знімаються при даних установках на конкретної ЕОМ, "Коеф. масштабування" і "Зміщення" – відповідають за перетворення даних, що надходять у випадку зняття даних зовнішнім пакетом. Зовнішній пакет видає дані в діапазоні від –2048 до 2048. У даній же програмі дані представляються 9-ю бітами, тобто діапазоном чисел від 0 до 511, "Точне число даних", "Моніторинг" та "Зняття" – у випадку малопродуктивного комп'ютера всій програмі необхідно дати найвищий пріоритет tpTimeCritical і скасувати команду Sleep, однак тоді програма не буде реагувати на дії користувача, а значить він не зможе зупинити процес зняття даних [20]. Для цього включається режим "Точне число даних", у ході якого зчитується визначене число дискретів в обох режимах, "Застосовувати відразу" – зміни внесені користувачем в опції будуть відразу ж враховані програмою, "Використовувати Sleep" – програма після кожної мс віддає комп'ютерний час операційному середовищу, що навіть при найвищому пріоритеті роботи програми (tpTimeCritical ) дозволяє користувачу взаємодіяти з програмою, "Затримка після зняття зі шкіри " – визначає частоту дискретизації даних. Частота обернено пропорційна даному параметру, "Усереднювати при знятті" – виконувати фільтрацію на етапі зняття даних. На якість впливає параметр "Кількість точок при фільтрації", описаний вище, "Дозволити доповнювати каналами" – у виключеному режимі, дані сусідніх неактивних каналів будуть знищені, "Процесор Pentium" – у включеному режимі на машині тільки з процесором intel Penium виконує більш точну прив'язку вчасно при зчитуванні даних, "Пріоритет потоку зняття даних" – установлює пріоритет потоку зняття даних (Вправо – найвищий пріоритет), "Кольори" – визначають набір кольорів використовуваних у програмі, "Графік" – визначає параметри відображення графіка: "Горизонтальні шкали", "Вертикальні шкали" – число шкал на графіку, "Чутливість "миші"" – визначає число точок, на які покажчик може відхилитися від маркера, не "загубивши" цей маркер, "Помічати дані при великому збільшені" – при великому масштабі дані позначаються хрестиком, "Маркери лише до значення" – при виключеному режимі маркери малюються не тільки до значення, а у всю висоту даного каналу, "Навігатор" – опції щодо навігатора: "Помічати в навігаторі кольором елементи з даними" – об'єкти (стрілочки), що містять дані будуть виділені червоним якщо вони містять дані.


ВИСНОВКИ

 

В даній роботі розглянуті принципи взаємодії оптичного випромінювання з біотканиною, що дозволяє стверджувати про доцільність використання фотоплетизмограми (ФПМ) для оцінки мікроциркуляцій крові. Зняття двох ФПГ з двох розташованих рядом ділянок шкіри дозволяє відкинути такі фактори як колір, температуру і деколи геометрію поверхні шкіри. Абсолютний та порівняльний аналіз двох ФПГ дозволяє стверджувати про порушення кровоносної системи, присутніх в одному з каналів.

Перегляд існуючих методів дозволяє стверджувати складність діагностування вертеброзахворювань особливо в вітчизняних умовах стану медицини. Проведена критеріальна шкала показує доцільність застосування нового неінвазивного методу в медичній практиці. Розробка та вдосконалення неінвазивних дешевих методів діагностики та моніторингу може покращити стан діагностики в регіональних лікарняних установах. У роботі фільтрація вихідного сигналу відбувалась при допомозі Фур′є-аналізу, але було запропоновано більш гнучкий метод обробки сигналу вейвлетний аналіз. Це один із перспективних методів математичного аналізу . При допомозі порівняння були обґрунтовано доведені його переваги в порівнянні з Фур’є-аналізом. Різновидність вейвлет-аналіза у випадку багато масштабного представлення свідчить про широке застосування даного аналізу в усіх сферах обробки сигналів.

Описана в останньому розділі апаратно-програмна реалізація комплексу реалізує запропонований метод та відповідає висунутим вимогам. Даний комплекс є зручним та недорогим інструментальним засобом оцінки порушення мікроциркуляції в судинних системах.


ЛІТЕРАТУРА

 

1.      Малая Л.Т., Микляев И.Ю., Кравчук П.Г. Микроциркуляция в кардиологии. - Харьков: Вища школа, 1977. – 142c.

2.      Добеши И. Десять лекций по вейвлетам перевод – Ижевск НИЦ Регулярная и хаотическая динамика М: 2001 , 464 с.

3.      Засецкий А В, Иванов А Б, Постников С Д, Соколов И В. Контроль качества в телекоммуникациях и связи М: 2001, 335 с.

4.      Павлов С.В. Разработка оптоэлектронного комплекса для преобразования биомедицинской информации. Автореферат диссертации кандидата наук, 1995. - 16c.

5.         Чандпасекар С. Перенос лучистой энергии – М., 1953.

6.         Кузьміч В.В. Основные принципы и особенности транскутанной "отражательной" оксиметрии // Мед. Техника. – 1993. - №3. – С. 36-42.

7.         Van Gemert M.I.C. Star W.M. // IEEE Trans. Biol. Med. Eng., 1989 – Vol. 36. №2– P.1146-1154.

8.         Samsung electronic. S3C44B0X RISC MICROPROCESSOR datasheet.

9.         Samsung electronic. K6R4016V1C dataseet.

10.      ГОСТ 23751-79. Платы печатные. Общие технические условия.

11.    Новиков Л.В. Основы вейвлет-анализа сигналов. Учебное пособие. Санкт-Петербург: 1999, 152 с.

12.    Л. Л. Маслюк, А.М. Переберен Введение в вейвлет-анализ. Учебный курс М: 1998, 370с.

13.    http://prodav.narod.ru/

14.    Шикин Е В, Боресков А В. Компьютерная графика. Динамика, реалистические изображения. Диалог – Мифи, М: 1996, 289 с.

15.      Разевиг В.Д. Применение программ PCAD и PСSpiсe для схемотехнического моделирования на ПЭВМ: В 4 выпусках.- М: Радио и связь, 1992. – 56c.

16.      Яншин А.А. Теоретические основы конструирования и надежности ЭВА, - М.: Радио и связь. - 1983. - 276 с.

17.      Павлов С.В., Превар А.П., Матохнюк М.В., Чернуха А. Застосування оптико-електронних та лазерних технологій при аналізі мікроциркуляторних змін у вогнищі гострого гнійного запалення в ділянці нижніх кінцівок / OPTOELECTRONIC INFORMATIJN-POWER TECHNOLOGIES / Вінниця, 2002.- №2 - С.148.

18.      Павлов С.В., Мисловський І.Д., Ахмед Авад. Оптико-електронні інформаційні технології дослідження рівня перифепійного кровообігу/ Матеріали I Всеукраїнської науково-практичної конференції “Медичні технології і вища освіта”.- Луцк, 2004. – С. 137-144.

19.      Петрук В.Г., Черноволик Г.О., Шевчук С.В., Безсмертний Ю.О. Математична модель перетворення проміння при поверхневим шаром біотканини з системною патологією / OPTOELECTRONIC INFORMATIJN-POWER TECHNOLOGIES / Вінниця, 2002. -№2 - С.154.

20.      Н.І. Заболотна, С.В. Павлов, В.В. Шолота Комп’ютерне моделювання задач лазерної та оптоелектронної техніки /Вінниця-ВНТУ, 2003. – 149 с.


Додаток А (Віконне перетворення сигналу)


Додаток Б (Трьохвимірна спектограма)


Додаток В (Лістинг програми)

 

procedure TData. AutoSetMarkers;

const

SC=513;

var

X,

x1,x2, // кінцеві координати для циклів

x11,x22, // поправочні для попередніх

xmp1,xmp2 // координати для максимуму производной

: Longint;

XI : Array[1..5] of integer; // буде містити координати // перепадів 3-х к ряду

sign : shortint; // важливий лише знак

max,maxp, // знайдене максимальне значення

HMW // скільки вже було

: integer;

was : Boolean;

A : Array of integer;

begin

x1:=0;

x2:=fcount-1;

Was:=false;

SetLength(A,fcount);

FillChar(XI,SizeOf(XI),0);

for x:=0 to fcount-1 do begin

if Markers[x] then

if not was then begin x1:=x; was:=true end

else begin x2:=x; end;

A[x]:=0;

end;

ClearAllMarkers; // очищаємо всі маркери

while Der[1,x1]=0 do inc(x1);

sign:=Der[1,x1];

x11:=0;

for x:=x1+1 to x2-1 do

if (Sign*Der[1,x]<0) then begin

// знайшли перепад між - и +

if x11=0 then x11:=x;

x22:=x;

A[x]:=SC; // помічаємо перепад

sign:=-sign

end;

x1:=x11;// скорочуємо діапазон

x11:=0;

x2:=x22;

XI[1]:=x1;

for x:=x1+1 to x2 do

if A[x]=SC then

if XI[2]=0 then XI[2]:=x

else begin

XI[3]:=x;

if x11=0 then x11:=xi[2];

if (items[XI[2]]-items[XI[1]]=0) or

(abs((items[XI[3]]-items[XI[2]])/(items[XI[2]]-items[XI[1]]))<=GO.MarkSens/100)

then continue;

A[XI[2]]:=items[XI[2]]-items[XI[1]];

A[XI[2]+1]:=items[XI[3]]-items[XI[2]];

XI[1]:=XI[2];

XI[2]:=XI[3];

end;

max:=0;HMW:=0; // знаходимо максимум по которому відловлювати

for x:=x1 to x2 do

if A[x]=SC then A[x]:=0;

for x:=x1 to x2 do

if (A[x]>0)and(A[x]<>SC) then begin

inc(HMW);

if A[x]>max then max:=A[x];

if HMW>=30 then break;

end;

x:=x1;

while x<=x2 do

if (abs(max-A[x])<=max*0.3) then begin

FillChar(XI,SizeOf(XI),0);

XI[1]:=x;

for x11:=x+1 to x2 do // знаходимо місце для 3-го маркеру

if A[x11]*A[x]>0 then begin

XI[3]:=x11;

break

end;

if XI[3]=0 then begin inc(x);continue;end;

maxp:=0;xmp1:=0;xmp2:=0;

for x11:=XI[1]+1 to XI[3] do // знаходимо місце для 2-го маркеру

if Der[1,x11]>maxp then begin

maxp:=Der[1,x11];

xmp1:=x11;

xmp2:=x11;

end else

if Der[1,x11]=maxp then xmp2:=x11;

XI[2]:=(xmp1+xmp2) div 2;

for x11:=XI[3]+2 to x2 do // знаходимо місце для 4-го маркеру

if A[x11]*A[XI[3]+1]>0 then begin

XI[4]:=x11;

break

end;

if XI[4]=0 then begin inc(x);continue;end;

for x11:=XI[4]+2 to x2-2 do // знаходимо місце для 5-го маркеру

if (Der[1,x11]<0)and(Der[1,x11+1]<0) then begin

XI[5]:=x11-1;

break

end;

if XI[5]=0 then begin inc(x);continue;end;

// корекція

XI[1]:=XI[1]-2;

For x11:=xi[2]-2 downto xi[1]+1 do

if (Der[1,x11]=0)or(Der[1,x11-1]*Der[1,x11]<0){or

(Der[2,x11]=0)or(Der[2,x11-1]*Der[2,x11]<0)} then begin

XI[1]:=x11-1;

break;

end;

XI[3]:=XI[3]-1;

XI[4]:=XI[4]-1;

for x11:=1 to 5 do

Markers[XI[x11]]:=True;

max:=(max+A[x])div 2;

x:=XI[5]+1;

end

else inc(x);

end;