Модули статистической обработки анализатора 'Тензотрем'
Содержание
Введение
1. Тензометрический метод оценки состояния двигательных отделов
центральной нервной системы
1.2 Структурная организация тензометрического треморографа
1.3 Измерительный блок
1.4 Устройство управления
1.4 Структурная организация аппаратного обеспечения
тензометрического треморографа
1.5 Варианты реализации
2. Задачи статистической обработки изометрических данных
2.1 Алгоритм преобразования Фурье
2.2 Метод корреляции
2.3 Метод главных компонент
2.4 Построение графика функции
2.5 Чтение данных из файла
3. Технико-экономическое обоснование
Заключение
Список использованных источников
Приложения
Реферат
Цель работы - исследование и разработка программных модулей
статистической обработки измерительной информации тензометрического треморографа.
Объект исследования - тензометрический треморограф.
Тензометрический треморограф предназначен для оценки
активности моторной системы человека по параметрам тремора (непроизвольного
сокращения мышц) рук пациента на основе регистрации изометрического усилия.
Конструкция и схемотехника тензометрического треморографа
должны обеспечивать:
- работу в условиях стационарной
лаборатории;
- устойчивость к транспортной тряске и
вибрации;
- технологичность конструкции;
- удобство эксплуатации.
Были рассмотрены несколько вариантов реализации модулей
статистической обработки: MS Excel, язык программирования С, MathLab.
Выбран язык С, ввиду его близости к машинному уровню, наличия
бесплатного компилятора, возможности написания программы без использования
больших библиотек.
В конечном итоге целью разработки анализатора является
автономный прибор на базе микроконтроллера, что налагает ограничения на
используемые алгоритмы и методы их реализации.
Социально-экономический эффект применения тензометрического
треморографа - раннее обнаружение и более точная диагностика двигательных
нарушений, возможность выбора оптимальных методов лечения и своевременная их
коррекция, связанная с подбором действенной лекарственной терапии.
Определения,
обозначения, сокращения
БПФ - быстрое преобразование Фурье.
Реверс битов (Обратная перестановка) - преобразование
двоичного кода в противоположный: все нули в коде заменяются на единицы,
единицы, соответственно заменяются нулями.
ПК - персональный компьютер.
МГК - метод главных компонент.
ВКР - выпускная квалификационная работа.
БОС - биологическая обратная связь.
корреляция анализатор тензометрический треморограф
Введение
Тензометрический треморограф предназначен для оценки
активности моторной системы человека по параметрам тремора (непроизвольного
сокращения мышц) рук пациента на основе регистрации изометрического усилия.
Основные характеристики прибора:
- диапазон произвольных усилий - ± (50÷10000) г;
- время регистрации - от 10 с до 5 мин.
(типовое значение 30 с);
- частота съема информации - (100 ÷ 400) Гц;
- время задержки съема информации - от 1 до
10 с;
- чувствительность - 1 г;
- питание прибора - 220 В, 50 Гц.
Область применения:
- медицина (клиническая диагностика
заболеваний центральной нервной системы, оценка и регулярный контроль профессиональной
пригодности, спортивная медицина, научно-исследовательская деятельность в
области нейронных механизмов организации движений, экспресс-мониторинг
различных групп населения для раннего выявления патологических состояний
центральной нервной системы);
- экология (исследования влияния на человека
различных факторов среды).
Экономическая эффективность заключается в создании прибора
базирующегося на принципиально новом методе регистрации тремора подвижных
звеньев тела человека, позволяющего объективно оценивать состояние управляющих
структур мозга, удобного в использовании и доступного для медицинских
лабораторий любого уровня.
Социально-экономический эффект применения тензометрического
треморографа - раннее обнаружение и более точная диагностика двигательных
нарушений, возможность выбора оптимальных методов лечения и своевременная их
коррекция, связанная с подбором действенной лекарственной терапии.
Объектом исследования является прибор, предназначенный для
оценки активности моторной системы человека по параметрам тремора
(непроизвольного сокращения мышц) рук пациента на основе регистрации
изометрического усилия - тензометрический треморограф.
Цель работы - исследование и разработка программных модулей
статистической обработки измерительной информации тензометрического
треморографа. Необходимо реализовать методы математической обработки данных
анализатора: преобразование Фурье, метод главных компонент, нахождение
корреляции и автокорреляции.
Преобразование Фурье помогает выделить собственное и непроизвольное
усилие, а соответственно и наличие или отсутствие заболевания.
Корреляционная функция позволяет наблюдать взаимосвязь работы
рук пациента, а автокорреляционная функция показывает взаимосвязь во времени
для одной выбранной руки пациента с ней же.
Метод главных компонент позволяет уменьшить размерность
данных, потеряв при этом минимум информации.
Входными данными является информация с тензодатчиков о
прикладываемом усилии, снимаемая с заданной частотой заданное время.
1.
Тензометрический метод оценки состояния двигательных отделов центральной
нервной системы
Рассмотрим новый метод [1], [2] оценки активности моторной
системы человека.
Индикатором патологического состояния центральной нервной
системы является повышенный видимый невооруженным глазом тремор.
Тремор возникает в результате непроизвольного сокращения
мышц. Повышение амплитуды и изменение частоты и формы колебаний относительно
нормы (частота и амплитуда физиологического тремора) служат признаками
нарушения центральных и периферических нейронных механизмов регуляции движений.
Анализ этих параметров важен как для понимания роли дисфункции отдельных
структур мозга в процессах управления движениями, так и для клиники в аспектах
раннего обнаружения, более точной диагностики двигательных нарушений и выбора и
коррекции оптимальных методов лечения, связанных с подбором действенной
лекарственной терапии.
Таким образом, анализ форм тремора, как индикатора состояния
двигательных отделов центральной нервной системы, имеет большое значение в
неврологической практике.
Известные методы исследования активности моторной
(двигательной) системы человека не позволяют объективно оценивать состояние
управляющих структур мозга. При этом обеспечивается только качественная оценка
патологических состояний и затруднена дифференциация форм тремора, возникающего
в результате поражения различных структур мозга.
Остается много неразрешенных вопросов, касающихся клинической
вариабельности тремора, течения и исходов ряда заболеваний экстрапирамидной
системы и мозжечка, вызывающими нарушение выполнения движений и часто
сопровождаемыми патологическим тремором. Терапия таких нарушений представляет
значительные трудности из-за отсутствия методов объективного контроля
параметров треморного синдрома.
Разрабатываемый прибор будет базироваться на новом методе
оценки активности моторной системы человека по параметрам колебаний
изометрического усилия и позволит выявить характерные признаки в параметрах
тремора, связанных с функционированием пирамидной и экстрапирамидной систем,
что обеспечит возможность дифференциальной диагностики патологических состояний
и контроль эффективности различных видов терапии.
Указанная методика, позволяет, благодаря регистрации
параметров произвольного изометрического усилия методом тензометрии, оценивать
активность структур центральной нервной системы, управляющих моторной функцией.
Непроизвольный компонент произвольно удерживаемого усилия
На рисунке 1.1 представлены фрагменты записей непроизвольного
компонента произвольно удерживаемого усилия, средняя величина которого за
интервал регистрации 30 с указана цифрой над каждой осциллограммой (кг), для
левой (слева) и правой (справа) рук. А - здоровый испытуемый; Б, В и Г -
пациенты с диагнозом, соответственно, синдром паркинсонизма, болезнь Паркинсона
и симптоматический тремор. Ось абсцисс: время, с. Ось ординат: усилие, кг.
1.2
Структурная организация тензометрического треморографа
Конструкция и схемотехника [3], [4] тензометрического
треморографа должны обеспечивать:
- работу в условиях стационарной
лаборатории;
- устойчивость к транспортной тряске и
вибрации;
- технологичность конструкции;
- удобство эксплуатации.
В составе тензометрического треморографа для выполнения
требуемых функций, необходимы следующие узлы и устройства.
Измерительная часть, включающая в себя:
- преобразователь усилия в электрический
сигнал;
- устройство нормирования сигнала
тензодатчика;
- аналого-цифровой преобразователь.
Устройство управления, осуществляющее сбор информации от АЦП
и передачу данных в персональный компьютер (ПК).
Также в состав прибора входит блок питания, представляющий
собой АС/DC преобразователь.
Анализ устройств на предмет их конструктивной реализации
показывает, что они могут быть выполнены в виде функционально законченных
модулей. Такое конструктивное решение позволит значительно повысить
технологичность изготовления и настройки изделия, упростить процедуру и
сократить сроки восстановительного ремонта.
1.3
Измерительный блок
Тензометрический треморограф, предназначенный для оценки
состояния двигательных отделов центральной нервной системы, путем исследования
параметров колебаний изометрического усилия, состоит из измерительной части и
персонального компьютера с соответствующим программным обеспечением.
Измерительная часть реализует следующие функции:
регистрацию текущего изометрического усилия;
первичное преобразование измерительной информации в форму,
доступную для передачи в ПК;
формирование блоков передаваемой информации в соответствии с
выбранным протоколом передачи;
передачи измерительной информации в компьютер в соответствии
с определенным алгоритмом.
Регистратор текущего изометрического усилия может быть
выполнен в различных вариантах: как в виде отдельных аппаратных модулей
тензорезисторов и мостовой схемы, так и с использованием промышленного
тензодатчика, имеющего в своем составе оба элемента.
Критерием выбора конкретной реализации служат такие параметры
как: режим работы, максимальная нагрузка, чувствительность, сопротивление
изоляции, габариты, пыле- и влагозащищённость и особенности конструкции.
Аппаратный регистратор текущего изометрического усилия должен
обеспечивать преобразование изометрического усилия пальцев рук пациента в
электрический сигнал.
В качестве сенсорного элемента для определения
изометрического усилия необходимо использовать тензометрический датчик [1].
В основе работы тензодатчиков лежит явление тензоэффекта
(тензорезистивный эффект), заключающееся в изменении активного сопротивления
проводников при их механической деформации.
1.4
Устройство управления
Устройство управления реализовано на базе микроконтроллераSilabs 8051F320 обеспечивает
выполнение следующих функций.
. Начальная инициализация аппаратных средств контроллера,
задействованных в процессе работы:
- инициализация системного таймера и
настройка системного тактового сигнала;
- настройка портов ввода/вывода, приоритетного
декодера матрицы;
- распределение функций встроенных таймеров
счетчиков;
- инициализация USB0 для связи с ПК;
- инициализация SPI для управления АЦП.
2. Разработка процедур, вызываемых по прерыванию:
- получение данных от АЦП по SPI,
- передача данных в ПК по USB.
Дополнительно разработан управляющий автомат для анализа
состояния шины USB и выполнения действий в соответствии с общепринятым
протоколом обмена, а также дескрипторов, с помощью которых этот обмен будет
производиться.
Разрабатываемое устройство обладает следующими свойствами:
- работает в режиме Full speed;
- относится к классу и работает по
протоколу, определяемому производителем;
- максимальный размер пакета для контрольной
точки в фазе данных - 64 байта;
- идентификаторы имеют некоторые значения:
VID, PID, DID;
- имеет строковые описания производителя и
продукта;
- питается от своего источника, источник
шины USB не задействован.
Программное обеспечение
При анализе результатов регистрации произвольно удерживаемого
изометрического усилия особое внимание уделяется программному обеспечению для
ПК, на которое возлагаются следующие задачи:
- обработка массивов информации с получением
аналитических и диагностических результатов;
- статистическая обработка результатов;
- участие в организации обратной связи с
пациентом;
- представление результатов исследований в
удобной форме и пр.
Структура
программного обеспечения
Структура программного обеспечения, так же, как и аппаратных
средств, определяется возлагаемой функциональной нагрузкой и условиями эксплуатации,
а также, в значительной мере, и условиями самого проектирования.
В нашем случае программное обеспечение предназначено для
прибора, который будет внедрен в практику медицинских учреждений для
диагностики состояния центральной нервной системы относительно новым методом
оценки активности моторной системы человека по параметрам колебаний
изометрического усилия. Отсюда следуют два очевидных важных обстоятельства:
- принимая во внимание высокую
диагностическую значимость, необходимо обеспечить высокую достоверность
измерительных и аналитических процессов;
- учитывая новизну метода, необходимо дать
возможность модернизировать программное обеспечение в процессе эксплуатации.
Одним из основных признаков структуры программного
обеспечения, соответствующей этим требованиям, является модульность.
Разрабатываемая структура программного обеспечения должна
быть такой, чтобы любую её часть можно было изменить, не вызывая изменений в
остальных частях программы. Принцип построения таких структур известен, и носит
название "принцип IPO (Input-Process-Output) - вход-процесс-выход".
Задача состоит в выявлении таких программных модулей, в таком распределении
функциональной нагрузки, чтобы каждый такой блок-модуль имел один вход и один
выход - на входе программный модуль получал бы определенный набор исходных
данных, выполнял содержательную обработку и возвращал бы один набор данных
результата.
Каждый такой блок-модуль должен представлять собой
ограниченный по сложности и объему программный элемент, не зависящий логически
от работы других программных блоков. Межмодульные информационные связи должны
быть предельно ослаблены, взаимодействие между ними в рамках основного процесса
минимизируется. (Однако в случае организации резервирующих процессов или для
автотестирования эти связи могут быть активизированы). При этом важнейшим
свойством модуля является его локальная функциональная полнота и завершенность,
необходимая и достаточная для выполнения отдельного элемента функциональной
спецификации.
Качество разрабатываемой модульной структуры можно оценить по
степени независимости:
- от алгоритма программы, который будет
разработан с учетом выбранной модульной структуры;
- от аргументов (параметров) модуля;
- от структуры таблиц, протоколов и т.п.;
- от структуры и формата базы данных.
Программа, выполняемая на основе модульной структуры,
приобретает те свойства, которые необходимы для наших условий проектирования и
применения, а именно:
- модульную программу легко составлять и
отлаживать: функциональные компоненты такой программы могут быть написаны и
отлажены порознь;
- модульную программу легче адаптировать,
сопровождать и модернизировать;
- модульную программу проще разрабатывать в
связи с возможностью перераспределения независимых модулей ограниченного объема
среди программистов с учетом сложности модуля и квалификации исполнителя.
Применим известный принцип абстракции и будем рассматривать
структуру программного обеспечения, подлежащую разработке, как многоуровневую.
Поэтому естественным следует считать применение технологии проектирования
"сверху вниз", то есть путем постепенной, пошаговой детализации
выполняемой программы, когда на каждом этапе детализации должны быть
рассмотрены различные варианты решения и выбраны близкие к оптимальным. При
этом критериями оптимизации являются, кроме указанного выше принципа
модульности, также возможность сводить сложные задачи к простым, а также
критерий открытости для оперативных изменений, характеризующий прозрачную
программу, снабженную подробными комментариями.
В то же время критерий формализации, требующий строгого
математического описания каждого модуля структуры, представляется нам хотя и
желательным, но не приоритетным.
Исследование
вариантов программного обеспечения
При исследовании вариантов структурной организации
программного обеспечения тензотреморографа выделены следующие программные
модули, ответственные за реализацию основных функций прибора:
- блок ввода данных оператором (исходные
данные пациента, параметры проводимого теста, исходные данные о лице или
организации, проводящем обследование пациента, дата проведения теста,
комментарий);
- блок ввода результатов измерения
изометрического усилия непосредственно от измерительной системы через USB-порт
компьютера;
- блок визуализации сигнала обратной связи;
- блок фиксации на энергонезависимом
носителе (например, на жестком диске компьютера или флэш-накопителе) серии
отсчетов в течение заданного времени тестирования с заданным интервалом
квантования;
- блок выделения произвольного и
непроизвольного компонентов изометрического усилия и определения амплитуды,
формы и частотного спектра непроизвольных колебаний;
- блок специальных методов анализа с целью
сопоставления параметров зарегистрированных спектров колебаний вариантам нормы
и различным формам патологии обследуемых пациентов;
- блок выбора режима работы прибора и
управления процессом тестирования пациентов;
- блок встроенной функциональной диагностики
аппаратных компонентов прибора.
При этом просматриваются три основные подсистемы
тензометрического треморографа с различным сочетанием аппаратных и программных
средств для реализации соответствующих функций:
- подсистема ввода-вывода, ориентированная
преимущественно на использование терминала компьютера как основного аппаратного
средства;
- подсистема управления проведением тестов,
сочетающая в себе собственные аппаратные средства прибора (опорные элементы,
систему измерения) и соответствующие программные средства, необходимые для
проведения тестов и съема первичной информации;
- подсистема математической обработки
данных, использующая исключительно вычислительные ресурсы и память
персонального компьютера.
1.4
Структурная организация аппаратного обеспечения тензометрического треморографа
Системное проектирование программного обеспечения
предполагает разработку программного продукта от общего к частному, путем
разбиения системы на иерархические единицы. Результатом этого этапа
проектирования является модульная структура, причем разработанные модули должны
удовлетворять следующим требованиям:
модули должны быть автономными и иметь независимые
функциональные характеристики;
ссылка на модуль должна выполняться по идентификатору, в том
числе и из других модулей;
окончив работу модуль должен возвращать управление тому
модулю, которым был вызван;
модуль должен иметь один вход и выход;
модуль должен быть компактным, чтобы обеспечивалось его
легкое восприятие и анализ результатов.
Оптимальная форма представления результата по модульной схеме
представления структуры ПО - схемы прохождения данных.
На основании функциональной спецификации тензометрического
треморографа была разработана структурная организация прибора в виде
функциональной (структурной) схемы.
Представленная на рисунке 4 структурная схема содержит
условное изображение функциональных блоков (модулей) прибора и связей между
ними.
Структурная схема тензометрического треморографа
Тензометрический треморограф содержит опорный элемент 1,
измерительный элемент 2 в виде мостовой схемы из тензорезисторов, одна
диагональ которой подключена к источнику опорного напряжения, а вторая
диагональ - к входам дифференциального усилителя 3 для преобразования
механического усилия в электрический сигнал. Выход дифференциального усилителя
3 соединен с входом интегрирующего (сглаживающего) усилителя 4 и первым входом
дополнительного дифференциального усилителя 5, второй вход которого соединен с
выходом интегрирующего усилителя 4. Выходы всех усилителей подключены к
соответствующим входам аналого-цифрового преобразователя 6. Выход
многоканального аналого-цифрового преобразователя 6 соединен с вычислителем
(компьютером) 7. Преобразованные в вычислителе 7 сигналы поступают в качестве
сигналов обратной связи на соединенный с вычислителем 7 блок сенсорных
раздражителей (сигналов) 8, а результаты обработки выводятся на подключенные к вычислителю
7 средства отображения информации 9. Сигнал обратной связи 10 с выхода
интегрирующего усилителя 4 может непосредственно или через вычислитель
отображаться в блоке сенсорных раздражителей 8, воздействуя на человека 11,
который прикладывает усилие к опорному элементу 1. Количество измерительных
блоков 12 соответствует числу одновременно регистрируемых звеньев тела, а
выходы усилителей каждого из измерительных блоков подключены к соответствующим
входам аналого-цифрового преобразователя 6, в который также могут быть введены
по другим входам 13 аналого-цифрового преобразователя 6 физиологические
параметры (например, электрокардиограмма), характеризующие состояние
испытуемого или применяемые к нему воздействия.
Вычислительный блок выделен на рисунке пунктиром, так как
именно его функциям посвящена данная работа.
Работа тензометрического треморографа заключается в
следующем. Пациент 11 надавливает пальцами рук на опорные элементы 1
соответствующих измерительных блоков 12-а и 12-6 и по сигналам обратной связи
(например, визуальным) с блока сенсорных раздражителей 8 следит за уровнем
собственного усилия, удерживая в соответствии с инструкцией на определенном
уровне метки на экране, смещение которых пропорционально уровню произвольного
усилия, поступающего в блок сенсорных раздражителей 8 с выходов интегрирующих
усилителей 4 как сигналы обратной связи 10 отдельно для левой и правой руки.
Структурная схема программного обеспечения приведена на
рисунке 2.2.
Модуль "Биологическая обратная связь (БОС)"
отображает на экране компьютера динамику движений конечностей человека. Вход
модуля - результаты по произвольному изометрическому усилию подготовленные
информационные материалы, выход - визуальное отображение действий человека.
Модуль "Регистрация" обеспечивает фиксацию тензометрических
показателей анализируемого человека. Вход модуля - тензометрические
измерительные данные по произвольному и непроизвольному изометрическому
усилиям, выход - результаты измерений, приведенные к определенной
структурированной форме представления. Модуль "Результаты" реализует
отбор полученных результатов измерений непроизвольного усилия в соответствии с
установками измерений. Вход модуля - результаты регистрации непроизвольного
усилия, выход - структурированные результаты тензоизмерений.
Модуль "Статистическая обработка" реализует
статистическую обработку структурированных данных: спектр сигналов на различных
временных участках, области, определенные методом главных компонент, функции
корреляции и автокорреляции. Вход модуля - структурированные результаты
регистрации, выход - результаты статистической обработки.
Структурная схема программного обеспечения
1.5 Варианты
реализации
Были рассмотрены несколько вариантов реализации: MS Excel, язык программирования
С, MathLab.
Реализация в MS Excel
Достоинства
- легкое наглядное представление результата;
- относительно низкая стоимость.
Недостатки
- относительно низкая скорость вычисления;
- невозможность переноса программы
куда-либо.
Реализация на языке С
Достоинства:
- возможность относительно легкого переноса
программного кода для выполнения на контроллере без участия ПК;
- С - компилируемый язык, соответственно
быстрее исполняется программа;
- возможность написания программы без
использования больших библиотек;
- наличие бесплатного компилятора.
Недостатки
- сложность реализации программы и
наглядного представления результатов работы;
- скорость разработки.
Реализация в MathLab
Достоинства
- легкость разработки и представления
результата;
- возможность "выгрузки"
программного кода, ориентированного на конкретную систему.
Недостатки
- большая стоимость.
В предлагаемой ВКР выбран вариант реализации на языке C, как
язык, наиболее близкий к машинному уровню, обеспечивающий быстрое выполнение
программного кода и не требующий затрат на покупку компилятора.
2. Задачи
статистической обработки изометрических данных
В конечном итоге целью разработки анализатора является
автономный прибор на базе микроконтроллера, что налагает ограничения на
используемые алгоритмы и методы их реализации. Например, нецелесообразно
подключать библиотеку математических функций из-за ее объема, так как
использоваться будет только ее малая часть, которую можно реализовать
самостоятельно.
2.1 Алгоритм
преобразования Фурье
Для ускорения получения результата реализован алгоритм
быстрого преобразования Фурье (БПФ).
Были рассмотрены два варианта алгоритма [5]: рекурсивный и
итерационный. Опишем их на алгоритмическом языке и затем проанализируем.
Рекурсивный алгоритм:
_FFT (a):
n = a. length // nявляется степенью2n==1an= e2iπ/n
w =1 [0] = (a0, a2,
…, an-2) [1] = (a1, a3, …, an-1)
[0] = Recursive_FFT (a [0]) [1] = Recursive_FFT (a [1])k
= 0 to n/2 - 1k = y [0] k + wy [1] k
yk+ (n/2) = y [0] k
- wy [1] k
w = wwn
return y
Итерационный алгоритм:
Iterative_FFT (a)_Reverse_Copy (a,A)
n = a. length // n является степенью 2s = 1 to lg (n)= 2s
wm= e2iπ/n
for k = 0 to n - 1 by m= 1j = 0 to m/2 - 1= wA
[k+j+m/2]= A [k+j][k+j] = u+t[k+j+m/2] = u-t
w = wwm
Процедура Bit_Reverse_Copy размещает данные исходного массива a в массиве A в требуемом порядке, то
есть осуществляет обратную перестановку или реверс битов. Например, если массив
а= (0, 1, 2, 3, 4, 5, 6,7), то на выходе получится массив А = (0, 4, 2, 6, 1,
5, 3,7). Нумерация элементов меняется согласно правилу реверса битов.
Каждый алгоритм имеет вычислительную сложность θ (nlg (n)). При этом рекурсивный алгоритм разбивает исходный массив на
части и сохраняет их копии, а итерационный копирует массив только один раз в
требуемом порядке. Кроме того итерационный алгоритм требует вычисления lg (n) и округления его до
целого числа, что так же задерживает выполнение преобразования.
В результате для реализации был выбран рекурсивный алгоритм.
Так как алгоритм рассчитан на отделение вычислительного
модуля от ПК и библиотека математических вычислений не используется, возведение
экспоненты в комплексную степень реализованачерез вычисление синуса.
Для вычисления синуса разложим его в ряд Тейлора, для чего
реализуем так же функцию возведения числа в степень и вычисление факториала.
Для работы с комплексными числами понадобится класс,
реализующий их в алгебраической форме, а так же функции сложения, вычитания и
умножения комплексных чисел.
Так как алгоритм рассчитан на количество отсчетов равное
степени двойки, а на вход приходит кратное десяти, то перед началом
преобразования дополним [6] сигнал нулями.
Продемонстрируем работу алгоритма для вектора из 8 элементов.
Рекурсивный алгоритм БПФ для вектора из 8 элементов
2.2 Метод
корреляции
Корреляция [7] - величина, характеризующая взаимную
зависимость двух случайных величин.
Корреляция двух функций f(x) иh(x) определяется выражением
Где 
- комплексно сопряженная функция к f(x). Так как в нашем случае
обе входные функции действительны, то 
.
Автокорреляция - статическая взаимосвязь между двумя
случайными величинами из одного ряда, но взятых со сдвигом.
По графику автокорреляционной функции можно судить о
периодичности базовой функции, а следовательно, о ее частотных характеристиках.
2.3 Метод
главных компонент
Метод главных компонент (МГК) [8] трансформирует описание
стационарного стохастического процесса таким образом, что выходное пространство
редуцированного размера сохраняет наиболее важную информацию об исходном
процессе. Применяется для уменьшения размерности пространства и выделения
наиболее значимой информации.
Пусть 
- случайный вектор, а 
обозначает ожидаемое
(среднее) значение матрицы автокорреляции по всем векторам. При конечном
количестве векторов эту матрицу можно описать следующим образом\
где матрица данных 
образована
последовательностью векторов 
.
Обозначим 
собственные числа матрицы автокорреляции, 
- сопряженные с ними
собственные векторы. Собственные значения и векторы связаны зависимостью
Где 
Собственные числа симметричной неотрицательной
матрицы, какой является 
, являются рациональными и неотрицательными.
Упорядочим их в последовательности убывания и, в соответствии
с этой последовательностью упорядочим собственные вектора. Возьмем первые собственных значений. В
этом случае матрица преобразования будет иметь вид 
, при 
. Эта матрица определяет
линейное преобразование 
.
Вектор 
представляет собой вектор главных компонент,
имеющих наибольшее влияние на реконструкцию вектора данных 
С точки зрения машинной реализации метода самой сложной
задачей является нахождение собственных чисел и векторов матрицы. [9]Для
решения этой проблемы чаще всего применяетсяQR метод поиска собственных
векторов и собственных чисел.
Реализация
Продемонстрируем реализацию методов статистической обработки,
а так же функцию построения графиков для визуализации результата, чтения данных
из файла и главную функцию, через которую осуществляются вызовы всех методов.
Полный текст программы приведен в приложении А.
Преобразование Фурье
Функция вычисления быстрого преобразования Фурье принимает на
вход указатели на массив исходных данных, массив результата и размерность этих
массивов целым числом.
// БПФ
// Аргументы:
// mas - указатель на массив исходных данных
// n - размерность массива mas (явялется степенью 2)
// y - указатель на массив выходных данных разиерности n
void FFT (double * mas, int n, complex * y)
{pi = 3.1415926535;(n==1)
{= &complex (mas [0], 0);;
}wn (sinx (2*pi/n + pi/2), sinx (2*pi/n));w (1,
0);i = 0, j = n/2, k = 0;*a0 = (double *) malloc (sizeof (double) * (n/2)),
*a1 = (double *) malloc (sizeof (double) *
(n/2));(k < j-1)
{[k] = mas [i];[k] = mas [i+1];++;= i+2;
}*y0 = (complex *) malloc (sizeof (complex) *
(n/2));*y1 = (complex *) malloc (sizeof (complex) * (n/2));(&a0 [0], j,
&y0 [0]);(&a1 [0], j, &y1 [0]);(k = 0; k < j; k++)
{[k] = add (y0 [k], mul (w,y1 [k]));[k+j] = sub
(y0 [k], mul (w,y1 [k]));= mul (w, wn);
}
}
Структура complex имеетвидstruct complex
{(double real, double imag)
{= real;= imag;
}re;
doubleim;
};
Вычисление синуса реализовано через разложение в ряд Тейлора,
для чего так же реализованы функции вычисления степени и факториала числа, но в
виду их простоты приводить мы их не будем.
Корреляция
Метод вычисления корреляции принимает указатели на два
массива исходных данных, массив выходных данных и размерность этих массивов
целым числом.
// Корреляция
// Аргументы:
// vect_1 - указатель на первый массив исходных данных
// vect_2 - указатель на второй массив исходных данных
// corr - указатель на массив выходных данных
// n - размерность данных
// примечание: предполагается, что размерности vect_1 и
vect_2 // совпадают
void correlation (double *vect_1, double *vect_2,
double *corr, int n)
{i, j;(i = 0; i < n; i++)
for (j = 0; j<n - i - 1; j++) // т к в выборке всего n
// значений, а индекс i+j не должен становиться больше
// n, то второй цикл до n - i
corr [i] = corr [i] + vect_1 [j] *vect_2 [i+j]
/n;
}
Пример построения графика корреляционной функции
На рисунке 5.1 представлен пример построения графика
корреляционной функции двух случайных векторов по 512 отсчетов. График
нормирован.
Автокорреляция
Метод вычисления автокорреляции принимает в качестве
аргументов указатели на массив входных данных, массив результата и размерность
этих массивов целым числом.
// Автокорреляция
// Аргументы:
// vect - указатель на массив исходных данных
// n - размерность массива mas (явялется степенью 2)
// acorr - указатель на массив выходных данных
void autocorrelation (double * vect, int n,
double * acorr)
{i, j;(i = 0; i < n; i++)
for (j = 0; j<n - i - 1; j++) // т к в выборке всего
// n значений, а индекс i+j не должен становиться
// больше n, то второй цикл до n-i
acorr [i] = acorr [i] + vect [j] *vect [i+j] /n;
}
Метод главных компонент
Аргументами метода являются указатели на массив исходных
данных, массив выходных данных и размерности обоих массивов целыми числами. Для
вычисления матрицы автокорреляции вектор-строку исходных данных транспонируем,
делим на несколько векторов и записываем их последовательно в матрицу.
// Метод главных компонент
// Аргументы:
// vect - указатель на одномерный массив исходных данных
// y - указатель на одномерный массив результата
// n - размерность входных данных
// k - размерность результата
void PCA (double *vect, double *y, int n, int k)
{, j, m = n/k, h;
// так как vect - вектор-строка, а алгоритм расчитан на
// вектор=столбец, сформируем вектор х=transp (vect)
double **x = (double **) malloc (sizeof
(double*));[0] = (double *) malloc (sizeof (double) * n);(i = 0; i < n;
i++)[0] [i] = vect [i];**Xm = (double **) malloc (sizeof (double *) * n);**XT =
(double **) malloc (sizeof (double *) * m);(i = 0; i < n; i++)[i] = (double
*) malloc (sizeof (double) * m);(i = 0; i < m; i++)[i] = (double *) malloc
(sizeof (double) * n);**Rxx = (double **) malloc (sizeof (double *) * n);(i =
0; i < n; i++)[i] = (double *) malloc (sizeof (double) * n);
// сформируем матрицу Х исходных векторов
// размерность матрицы определим как 1000хn/1000
h = 0;(i = 0; i < m; i++)(j = 0; j < n;
j++)
{[j] [i] = x [0] [h];++;
}(&Xm [0], &XT [0], n, m);(&Xm [0],
&XT [0], n, m, &Rxx [0]);(i = 0; i < n; i++)(j = 0; j < n;
j++)[i] [j] = Rxx [i] [j] /m;**R = (double **) malloc (sizeof (double *) *
n);**Q = (double **) malloc (sizeof (double *) * n);(i = 0; i < n; i++)
{[i] = (double *) malloc (sizeof (double) *
n);[i] = (double *) malloc (sizeof (double) * n);
}(&Rxx [0], n, &R [0], &Q [0]);*eig =
(double *) malloc (sizeof (double) * n);(i = 0; i < n; i++)[i] = R [i]
[i];(&eig [0], &Q [0], n);**W = (double **) malloc (sizeof (double*) *
k);(i = 0; i < k; i++)[i] = (double *) malloc (sizeof (double) * n);(i = 0;
i < k; i++)(j = 0; j < n; j++)[i] [j] = Q [j] [i];(i = 0; i < k;
i++)(j = 0; j < n; j++)[i] = W [i] [j] *x [0] [j];(x);(Xm);(XT);
free (Rxx);(W);
}
Для реализации метода главных компонент были так же
реализованы умножение и транспонирование матриц, QR метод поиска собственных
чисел и векторов, для которого потребовалось вычисление квадратного корня и
сортировка собственных чисел и сопряженных с ними собственных векторов.
Рассмотрим реализацию только QR-алгоритма и алгоритма
сортировки, так как все остальные алгоритмы просты и не требуют детального
рассмотрения.
QR
- алгоритм
В качестве аргументов метод принимает указатель на матрицу,
для которой будут вестись вычисления, указатели на результирующие матрицы RиQ, а так же размерность
всех матриц. Так как алгоритм реализован для квадратных матрицnxn, то размерность матриц
можно передать одним целочисленным параметром.
// QR-алгоритм
// Аргументы:
// A - указатель на входную квадратную матрицу nxn
// n - размерность данных
// R - указатель на результирующую матрицу с СЧ матрицы А на
// главной диагонали разиерности nxn
// Q - указатель на результирующую матрицу СВ матрицы А
// разиерностиnxnQR (double ** A, int
n, double ** R, double ** Q)
{**E = (double **) malloc (sizeof (double *) *
n);**T = (double **) malloc (sizeof (double *) * n);i, j, k, g;(i = 0; i <
n; i++)
{[i] = (double *) malloc (sizeof (double) *
n);[i] = (double *) malloc (sizeof (double) * n);
}(i = 0; i < n; i++)
{(j = 0; j < n; j++)
{[i] [j] = A [i] [j];[i] [j] = 1;(i == j)[i] [j]
= 1;[i] [j] = 0;
}
}(k = 0; k < (n-1); k++)
{(i = k+1; i < n; i++)
{(j = 0; j < n; j++)
{(g = 0; g < n; g++)
{[j] [g] = E [j] [g];
}
}(R [i] [k]! = 0)
{[k] [k] = R [k] [k] /sqrtx (powx (R [k] [k],2)
+powx (A [i] [k],2));[i] [i] = T [k] [k];[k] [i] = R [i] [k] /sqrtx (powx (R
[k] [k],2) +powx (A [i] [k],2));[i] [k] = T [k] [i];
}(&T [0], &R [0], n, n, &R
[0]);(&T [0], &Q [0], n, n, &Q [0]);
}
}(E);(T);
}
Сортировка
Модуль сортировки осуществляет перестановку массива
собственных чисел матирцы в порядке убывания и, в соответствии с ним,
перестановку сопряженных с собственными числами собственных векторов матрицы.
То есть при перемене местами iиjсобственных чисел, соответственно меняются местами так же и iиjсобственные вектора.
Модуль принимает в качестве аргументов указатели на массив
собственных чисел, матрицу собственных векторов и размерность массива. Так как
вычисления собственных значений велись для квадратной матрицы nxn, то в общем случае у нее
nсобственных чисел и nn-размерных собственных
векторов.
// Сортировка СЧ и СВ
// Аргументы:
// eig - указатель на массив СЧ
// Q - указатель на массив СВ
// n - размерностьмассивовsort (double *eig, double
** Q, int n)
{= false;(i = 0; i < n-1; i++)(eig [i] <
eig [i+1])
{= true;= eig [i];[i] = eig [i+1];[i+1] = a;(j =
0; j < n-1; j++)
{= Q [i] [j];[i] [j] = Q [i+1] [j];[i+1] [j] = a;
}
}
}
}
2.4 Построение
графика функции
Так как стандартная графическая библиотека языка C не может
быть применена при переносе реализации с ПК, этот метод в любом случае будет
модифицирован, поэтому воспользуемся более удобной свободной графической
библиотекой FLTK [10] языка C++.
В этом случае модуль построения графика функции выглядит
следующим образом.
// Построение графиков
// Аргументы:
// vect - указатель на массив значений
// n - количество значений
#include<FL/FL. H>
#include "Simple_window. h"
#include "Graph. h"plot (double * vect,
int n)
{namespace Graph_lib;int xmax = 600; // размер окнаint ymax =
400;int x_orig = xmax/2; // точка (0,0) - центр окнаint y_orig = ymax/2;Point orig (x_orig,
y_orig);int r_min = - 1000; // диапазонint r_max = 1000;int n_noits = (r_max - r_min)
*n; // количество точек в
// диапазоне_scale = 30; // масштабирующий множительint y_scale =
30;p1 (100,100);_window win (p1, 600, 400, "Plot"); // Создание окнаint xlen =
xmax - 40; // создаем оси
constintylen = ymax - 40; // делаем их чуть меньше
// размераокна(Axis:: x, Point (20,
y_orig), xlen, xlen/x_scale, "x");
// созданиекоординатныхосей(Axis:: y, Point (x_orig,
ylen + 20), ylen, ylen/y_scale, "y");
// созданиекоординатныхосей. attach (x);
win. attach (y);max = vect [0];i;(i = 0; i <
n; i++)(vect [i] > max)= vect [i];* y1 = (int*) malloc (sizeof (int) * n);(i
= 0; i < n; i++)[i] = vect [i] *ymax/max;_polyline poly;t = r_min;dt =
(r_max - r_min) /n;(i = 0; i < n; i++)
{. add (Point (t, y1 [i]));= t+dt;
}. set_color (Color:: black);. attach (poly); // добавление poly в окно.
wait_for_button (); // вывод окна
return;
}
2.5 Чтение
данных из файла
Результаты каждого теста помещаются в отдельный файл, в
котором каждый новый отсчет пишется в новой строке.
Данный модуль производит чтение результатов проведенного
теста для дальнейшей обработки.
// Чтение из файла
// Аргументы:
// vect - указатель на вектор значений, считываемых из файла
#include <stdlib. h>readFile (double *
vect)
{i;("%s","filename:
");filename [200]; // ={"I: \\1. txt"};("%s",
&filename);* file;= fopen (filename, "r");
i = 0;(file! = NULL)
{
while (true)
{= (double *) realloc (vect, sizeof (double) *
(i+1));(vect! = NULL)
{[i] = fscanf (file,"%f",&vect
[i]); // читаемизфайла
// очередноечисло
}(vect [i] == EOF || vect == NULL)
{(file);i;
}++;
}i;
}
}
Главная функция
// Главная функция
int _tmain (int argc, _TCHAR* argv [])
{f0 = true, f1 = true;
while (f0)
{= true;
// ввод имени файла и данных из него
printf ("%s", "Please enter the
file name\n");* vect = (double *) malloc (sizeof (double));n = readFile
(vect);
// выбор метода обработкиnum;("%s",
"Please select the method\n");("%s", "1.
FFT\n");("%s", "2. PCA\n");("%s", "3.
Correlation\n");("%s", "4. Autocorrelation\n");("%s",
"5. To attach another file\n");("%s", "6.
Exit\n");(f1)
{("%s", "Number:
");("%d", &num);
if (num == 1)
{= n;(n& (n-1)) // если побитовое и дает 1, то n // не //
равно степени 2
{(i& (i-1))
{= (double *) realloc (vect, sizeof (double) *
(i+1));[i] = 0;
}
}* y = (complex*) malloc (sizeof (complex) *
i);(vect, i, y);* y1 = (double*) malloc (sizeof (double) * i);(int j = 0; j
< i; j++)[j] = sqrt (powx (y [j]. re,2) + powx (y [j]. im,2));(y1,
n);(y);(y1);
}if (num == 2)
{* rPCA = (double*) malloc (sizeof
(double));(vect, rPCA, n, 1000);(rPCA, 1000);(rPCA);
}if (num == 3)
{("%s","Please enter the name of
the second file\n");* vect_1 = (double *) malloc (sizeof
(double));(vect_1);* corr = (double *) malloc (sizeof (double) * n);(vect,
vect_1, corr, n);(corr, n);(vect_1);(corr);
}if (num == 4)
{* acorr = (double *) malloc (sizeof (double)
*n);(&vect [0], n, &acorr [0]);(acorr, n);
free (acorr);
}
else if (num == 5)= false;if (num == 6)
{= false;
f0 = false;
}
}
}0;
}
3.
Технико-экономическое обоснование
Расчет полных затрат при выполнении исследования, разработки
прибора или устройства, написания программы, начинается с составления
детализированного плана работ, которые необходимо выполнить на каждом этапе
проектирования.
Под проектированием будем понимать совокупность работ,
которые необходимо выполнить, чтобы решить поставленную в задачу.
Для расчета затрат на этапе проектирования необходимо
определить продолжительность каждой работы (начиная с составления технического
задания (ТЗ) и до оформления документации включительно).
Результаты сведем в таблицу.
Длительность этапа разработки
|
№
|
Наименование
работы
|
Длительность
работы, дни
|
|
|
|
|
|
|
1
|
Изучение
литературы
|
5
|
8
|
6,5
|
|
2
|
Разработка
алгоритма
|
5
|
7
|
6
|
|
3
|
Написание
программного кода
|
5
|
7
|
6
|
|
4
|
Поиск и
исправление ошибок
|
5
|
8
|
6,5
|
|
5
|
Составление
отчета
|
5
|
6
|
5,5
|
|
6
|
Сдача проекта
|
5
|
7
|
6
|
|
Итого
|
25
|
43
|
36,5
|
Для расчета заработной платы примем в качестве месячной
заработной платы студента заработную плату инженера - программиста, что в
среднем составляет 31500 руб/мес (информация взята с сайта hh.ru для человека
без опыта работы).
В этом случае дневная ставка при 21 рабочем дне в месяце
составит 
.
Так как для написания работы другие исполнители не
привлекались, то основная заработная плота составит
,
а дополнительная заработная плата составит
,
где 
- норматив дополнительной заработной платы (%).
При выполнении ВКР норматив дополнительной заработной платы принимается равным
14%.
Отчисления на социальные нужды составляют
,
где 
- норматив отчисленинй на социальные нужды (%),
составляющий в настоящее время 30%.
Оценка потребности в ресурсах устанавливается в натуральных и
стоимостных показателях.
Себестоимость по каждой из статей можно рассчитать
нормативным методом.
Сырье и материалы:
Стоимость сырья и материалов рассчитывается по формуле
,
Где 
- затраты на сырье и материалы (руб.);
- индекс вида сырья или материала;
- норма расхода 
-того материала на единицу продукции (ед.);
- цена приобретения единицы -го материала (руб. /ед.);
- норма транспортно-заготовительных расходов (%).
При выполнении расчетов в ВКР норму транспортно-заготовительных
расходов ( ) принимаем равной 10%.
Результаты сведем в таблицу.
Затраты на сырье и материалы
|
Изделие
|
Материал
|
Тип (профиль,
сорт, марка, размер)
|
Норма расхода
на изделие
|
Цена за
единицу, руб.
|
Сумма на
изделие, руб.
|
|
Бумага
|
Бумага
|
А4, пачка
|
100 листов из
500
|
150
|
33
|
|
Картридж
|
Пластик,
чернила
|
Samsung
MLT-D101S
|
100 страниц из
1500
|
2150
|
36
|
|
Итого
|
2300
|
69
|
Покупки крупных комплектующих и услуги сторонних организаций
не имели места. Для написания программы используются основные средства в виде
компьютера и принтера, поэтому необходимо учесть и включить в затраты
амортизационные отчисления по этим основным средствам. Амортизационные
отчисления по основному средству i за год определяются как:
,
где 
- амортизационные отчисления за год по - му
основному средству (руб.);
- первоначальная стоимость -го основного средства (руб.);
- годовая норма амортизации -го основного средства (%), равная 20%.
Для компьютера амортизационные отчисления составят 11 руб/день.
Для принтера - 3, 44 руб/день.
Для определения величины амортизационных отчислений по основным
средствам, используемым в процессе выполнения ВКР необходимо определить время,
в течение которого студент использует это основное средство
Компьютер использовался со 2 по 5 этапы работы (24 дня), принтер -
только в 5 (5,5 дней).
Тогда для компьютера общие амортизационные отчисления составят 264
руб, для принтера - 18, 92 руб.
На последнем этапе расчета затрат на выполнение ВКР необходимо
определить величину накладных расходов.
Для расчета совокупной величины затрат, связанных написанием
программы, все проведенные расчеты оформим в таблицу.
Смета затрат на ВКР
|
№ п/п
|
Наименование
статьи
|
Сумма, руб
|
|
1
|
Расходы на
оплату труда
|
62415
|
|
2
|
Отчисления на
социальные нужды
|
18724,5
|
|
3
|
Материалы
|
69
|
|
4
|
Амортизационные
отчисления
|
382,92
|
|
5
|
Накладные
расходы
|
24966
|
|
Итого затрат
|
106557,42
|
Заключение
Выпускная квалификационная работа содержит в себе описание,
программную реализацию и тестирование различных модулей статистической
обработки тензометрического треморографа, включая визуализацию данных.
Реализованы БПФ, метод корреляции и автокорреляции, метод
главных компонент на языке программирования С. Язык реализации позволит
продолжить дальнейшую разработку прибора, подразумевающую перенос основных
вычислений с ПК в аппаратную часть, без осложнений, связанных с переносом
реализации алгоритмов из одной среды в другую.
Метод визуализации данных реализован на языке C++, а не C,
что тем не менее не должно повлиять на дальнейшую разработку программного
обеспечения.
Список
использованных источников
1. Романов,
С.П. Метод тензометрии для количественной оценки тремора /С.П. Романов, А.Ф.
Якимовский, М.Г. Пчелин // Физиол. ж. им. И.М. Сеченова. 1996, 82 (2). - С
118-123.
2. Романов,
С.П. Характеристики изометрически регистрируемого тремора при поражении экстрапирамидной
системы / С.П. Романов, А.Ф. Якимовский, М.Г. Пчелин // Российский
физиологический журнал им. И.М. Сеченова. 1997.83 (3). - C 133-139.
. Буренева,
О.И. Анализатор актуализации состояния человека [Текст] / П.Н. Бондаренко, О.И.
Буренева, Н.М. Сафьянников // Сб. тр. регионального совещания "Актуальные
проблемы психологической безопасности", Санкт-Петербург июнь 2012 г. -
СПб: Изд-во "Своё Издательство", 2012. - С.101 - 104.
. Разин,
Н.Н. Треморографическая система для изменения формы и амплитуды колебаний
подвижных звеньев конечностей человека [Текст] / З.А. Алексанян, О.Б. Кострова,
Н.Н. Разин, С.П. Романов // Сб. научных тр. Пятого международного научного
конгресса "Нейробиотелеком-2012", Санкт-Петербург декабрь 2012 г. -
СПб: Политехника, 2012. - С. 197. - 200.
. Кормен,
Т. Алгоритмы: построение и анализ, 3-е изд. / Ч. Лейзерсон, Р. Ривест, К. Штайн
// Пер. с анг. - М.: ООО "И.Д. Виьямс" 2013. - 1328с.
. Солонина
А.И. Основы цифровой обработки сигналов: Курс лекций. Изд.2-е испр. и перпраб.
/ Д.А. Улахович, С.М. Арбузов, Е.Б. Соловьева // СПб.: БХВ-Петербург, 2005 -
768с.
. Гонсалес,
Р. Цифровая обработка изображений / Р. Вудс // Москва, Техносфера, 2006 - 1072
с.
. Осовский,
С. Нейронные сети для обработки информации // Пер. с польского И.Д. Рудинского.
- М.: Финансы и статистика, 2004 - 344 с.
. Бахвалов,
Н.С. Численные методы - 5-е изд. / Н.П. Жидков, Г.М. Кобельков // М.: БИНОМ.
Лаборатория знаний, 2007 - 636 с.
. Страуструп,
Б. Программирование: принципы и практика использования C++, испр. изд. // Пер.
с англ. - "И.Д. Вильямс", 2011 - 1248 с.
. Отчето
выполнении НИОКР по теме: "Исследование и разработка структурной
организации, алгоритмического обеспечения и измерительной части
тензометрического треморографа" // ООО "Тензотрем",
Санкт-Петурбург, 2012 - 121с.
Приложения
Приложение А
Полный текст программы реализации модулей статистической
обработки анализатора "Тензотрем"
// diplom. cpp: определяет точку входа для консольного
приложения.
#include "stdafx. h"
#include <stdio. h>
#include <malloc. h>
// Транспонирование матрицыtransponse (double ** x, double ** xt)
{(int i = 0; i<512; i++)(int j = 0; j<512;
j++)
xt [i] [j] = x [j] [i];
}
// структура complex для работы с комплексными числами
// re - действительная часть числа
// im - мнимая часть числа
typedef struct complex
{(double real, double imag)
{= real;= imag;
}re;
doubleim;
};
// Функции работы с комплексными числами
// умножениеmul (complex a, complex
b)
{return complex (a. re*b. re - a. im*b. im, a.
re*b. im + a. im*b. re); }
// сложениеadd (complex a, complex b)
{return complex (a. re + b. re, a. im + b. im); }
// вычитаниеsub (complex a, complex
b)
{return complex (a. re - b. re, a. im - b. im); }
// вычисление степени числаpowx (double
x, int k)
{res = 1;(int i = 0; i < k; i++)= res*x;res;
}
// вычисление факториала(intn)
{res = 1;(int i = 1; i < n; i++)= res*i;res;
}
// разложение синуса вряд тейлора(doublex)
{res = 0;(int i = 0; i < 100; i++)= powx (-1,
i-1) *powx (x, 2*i-1) /fact (2*i-1) + res;
returnres;
}
// БПФ
// Аргументы:
// mas - указатель на массив исходных данных
// n - размерность массива mas (является степенью 2)
// y - указатель на массив выходных данных размерности n
void FFT (double * mas, int n, complex * y)
{pi = 3.1415926535;(n==1)
{= &complex (mas [0], 0);;
}wn (sinx (2*pi/n + pi/2), sinx (2*pi/n));w (1,
0);i = 0, j = n/2, k = 0;*a0 = (double *) malloc (sizeof (double) * (n/2)),
*a1 = (double *) malloc (sizeof (double) *
(n/2));(k < j-1)
{[k] = mas [i];[k] = mas [i+1];++;= i+2;
}*y0 = (complex *) malloc (sizeof (complex) *
(n/2));*y1 = (complex *) malloc (sizeof (complex) * (n/2));(&a0 [0], j,
&y0 [0]);(&a1 [0], j, &y1 [0]);(k = 0; k < j; k++)
{[k] = add (y0 [k], mul (w,y1 [k]));[k+j] = sub
(y0 [k], mul (w,y1 [k]));= mul (w, wn);
}
}
// Корреляция
// Аргументы:
// vect_1 - указатель на первый массив исходных данных
// vect_2 - указатель на второй массив исходных данных
// corr - указатель на массив выходных данных
// n - размерность данных
// примечание: предполагается, что размерности vect_1 и
vect_2 // совпадают
void correlation (double *vect_1, double *vect_2,
double *corr, int n)
{i, j;(i = 0; i < n; i++)
for (j = 0; j<n - i - 1; j++) // т к в выборке всего n
// значений, а индекс i+j не должен становиться больше
// n, то второй цикл до n - i
corr [i] = corr [i] + vect_1 [j] *vect_2 [i+j]
/n;
}
// Автокорреляция
// Аргументы:
// vect - указатель на массив исходных данных
// n - размерность массива mas (явялется степенью 2)
// acorr - указатель на массив выходных данных
void autocorrelation (double * vect, int n,
double * acorr)
{i, j;(i = 0; i < n; i++)
for (j = 0; j<n - i - 1; j++) // т к в выборке всего
// n значений, а индекс i+j не должен становиться
// больше n, то второй цикл до n-i
acorr [i] = acorr [i] + vect [j] *vect [i+j] /n;
}
// умножение матриц размера nxm
// Аргументы:
// А - указатель на первый множитель
// В - указатель на второй множитель
// n, m - размерность матриц
// R - указатель на результат
void mulMatr (double ** A, double ** B, int n,
int m, double ** R)
{i, j, k;(i = 0; i < n; i++)(j = 0; j < m;
j++)(k = 0; k < m; k++)[i] [j] = R [i] [j] + A [i] [k] *B [k] [j];
}
// Вычисление корня числа
// Аргументы:
// x - число, квадратный корень которого надо вычислить
// Примечание: вычисление ведется методом Ньютона(doublex)
{i;res = 1;(i = 0; i < 100; i++)= (res +
x/res) /2;res;
}
// QR-алгоритм
// Аргументы:
// A - указатель на входную квадратную матрицу nxn
// n - размерность данных
// R - указатель на результирующую матрицу с СЧ матрицы А на
// главной диагонали размерности nxn
// Q - указатель на результирующую матрицу СВ матрицы А
// размерности nxnQR (double ** A, int
n, double ** R, double ** Q)
{**E = (double **) malloc (sizeof (double *) *
n);**T = (double **) malloc (sizeof (double *) * n);i, j, k, g;(i = 0; i <
n; i++)
{[i] = (double *) malloc (sizeof (double) *
n);[i] = (double *) malloc (sizeof (double) * n);
}(i = 0; i < n; i++)
{(j = 0; j < n; j++)
{[i] [j] = A [i] [j];[i] [j] = 1;(i == j)[i] [j]
= 1;[i] [j] = 0;
}
}(k = 0; k < (n-1); k++)
{(i = k+1; i < n; i++)
{(j = 0; j < n; j++)
{(g = 0; g < n; g++)
{[j] [g] = E [j] [g];
}
}(R [i] [k]! = 0)
{[k] [k] = R [k] [k] /sqrtx (powx (R [k] [k],2)
+powx (A [i] [k],2));[i] [i] = T [k] [k];[k] [i] = R [i] [k] /sqrtx (powx (R
[k] [k],2) +powx (A [i] [k],2));[i] [k] = T [k] [i];
}(&T [0], &R [0], n, n, &R
[0]);(&T [0], &Q [0], n, n, &Q [0]);
}
}(E);(T);
}
// Трансспонирование матрицы размера nxm
// Аргументы:
// x - указатель на исходную матрицу
// xt - указатель на результирующую матрицу
// n,m - размерность матриц
void transp (double **x, double **xt, int n, int
m)
{i, j;(i = 0; i < n; i++)(j = 0; j < m;
j++)[i] [j] = x [j] [i];
}
// Сортировка СЧ и СВ
// Аргументы:
// eig - указатель на массив СЧ
// Q - указатель на массив СВ
// n - размерность массивовsort (double *eig, double
** Q, int n)
{i, j;fl = true;a;(fl)
{= false;(i = 0; i < n-1; i++)(eig [i] <
eig [i+1])
{= true;= eig [i];[i] = eig [i+1];[i+1] = a;(j =
0; j < n-1; j++)
{= Q [i] [j];[i] [j] = Q [i+1] [j];[i+1] [j] = a;
}
}
}
// Метод главных компонент
// Аргументы:
// vect - указатель на одномерный массив исходных данных
// y - указатель на одномерный массив результата
// n - размерность входных данных
// k - размерность результата
void PCA (double *vect, double *y, int n, int k)
{, j, m = n/k, h;
// так как vect - вектор-строка, а алгоритм расчитан на
// вектор=столбец, сформируем вектор х=transp (vect)
double **x = (double **) malloc (sizeof
(double*));[0] = (double *) malloc (sizeof (double) * n);(i = 0; i < n;
i++)[0] [i] = vect [i];**Xm = (double **) malloc (sizeof (double *) * n);**XT =
(double **) malloc (sizeof (double *) * m);(i = 0; i < n; i++)[i] = (double
*) malloc (sizeof (double) * m);(i = 0; i < m; i++)[i] = (double *) malloc
(sizeof (double) * n);**Rxx = (double **) malloc (sizeof (double *) * n);(i =
0; i < n; i++)[i] = (double *) malloc (sizeof (double) * n);
// сформируем матрицу Х исходных векторов
// размерность матрицы определим как 1000хn/1000
h = 0;(i = 0; i < m; i++)(j = 0; j < n;
j++)
{[j] [i] = x [0] [h];++;
}(&Xm [0], &XT [0], n, m);(&Xm [0],
&XT [0], n, m, &Rxx [0]);(i = 0; i < n; i++)(j = 0; j < n;
j++)[i] [j] = Rxx [i] [j] /m;**R = (double **) malloc (sizeof (double *) *
n);**Q = (double **) malloc (sizeof (double *) * n);(i = 0; i < n; i++)
{[i] = (double *) malloc (sizeof (double) *
n);[i] = (double *) malloc (sizeof (double) * n);
}(&Rxx [0], n, &R [0], &Q [0]);*eig =
(double *) malloc (sizeof (double) * n);(i = 0; i < n; i++)[i] = R [i]
[i];(&eig [0], &Q [0], n);**W = (double **) malloc (sizeof (double*) *
k);(i = 0; i < k; i++)[i] = (double *) malloc (sizeof (double) * n);(i = 0;
i < k; i++)(j = 0; j < n; j++)[i] [j] = Q [j] [i];(i = 0; i < k;
i++)(j = 0; j < n; j++)[i] = W [i] [j] *x [0] [j];(x);(Xm);(XT);
free (Rxx);(W);
}
// Построение графиков
// Аргументы:
// vect - указатель на массив значений
// n - количество значений
#include<FL/FL. H>
#include "Simple_window. h"
#include "Graph. h"plot (double * vect,
int n)
{namespace Graph_lib;int xmax = 600; // размер окнаint ymax =
400;int x_orig = xmax/2; // точка (0,0) - центр окнаint y_orig = ymax/2;Point orig (x_orig,
y_orig);int r_min = - 1000; // диапазонint r_max = 1000;int n_noits = (r_max - r_min)
*n; // количество точек в
// диапазоне_scale = 30; // масштабирующий множительint y_scale =
30;p1 (100,100);_window win (p1, 600, 400, "Plot"); // Создание окнаint xlen =
xmax - 40; // создаем оси
constintylen = ymax - 40; // делаем их чуть меньше
// размера окна(Axis:: x, Point (20,
y_orig), xlen, xlen/x_scale, "x");
// создание координатных осей(Axis:: y, Point
(x_orig, ylen + 20), ylen, ylen/y_scale, "y");
// создание координатных осей. attach (x);.
attach (y);max = vect [0];i;(i = 0; i < n; i++)(vect [i] > max)= vect
[i];* y1 = (int*) malloc (sizeof (int) * n);(i = 0; i < n; i++)[i] = vect
[i] *ymax/max;_polyline poly;t = r_min;dt = (r_max - r_min) /n;(i = 0; i <
n; i++)
{. add (Point (t, y1 [i]));= t+dt;
}. set_color (Color:: black);. attach (poly); // добавление poly в окно.
wait_for_button (); // вывод окна
return;
}
// Чтение из файла
// Аргументы:
// vect - указатель на вектор значений, считываемых из файла
#include <stdlib. h>readFile (double *
vect)
{i;("%s","filename:
");filename [200]; // ={"I: \\1. txt"};("%s",
&filename);* file;= fopen (filename, "r");= 0;(file! = NULL)
{(true)
{= (double *) realloc (vect, sizeof (double) *
(i+1));(vect! = NULL)
{[i] = fscanf (file,"%f",&vect
[i]); // читаем
из файла
// очередное число
}(vect [i] == EOF || vect == NULL)
{(file);i;
}++;
}i;
}
}
// Главная функция_tmain (int argc, _TCHAR* argv [])
{f0 = true, f1 = true;
while (f0)
{= true;
// ввод имени файла и данных из него
printf ("%s", "Please enter the
file name\n");* vect = (double *) malloc (sizeof (double));n = readFile
(vect);
// выбор метода обработкиnum;("%s",
"Please select the method\n");("%s", "1.
FFT\n");("%s", "2. PCA\n");("%s", "3.
Correlation\n");("%s", "4.
Autocorrelation\n");("%s", "5. To attach another
file\n");("%s", "6. Exit\n");(f1)
{("%s", "Number:
");("%d", &num);
if (num == 1)
{= n;(n& (n-1)) // если побитовое и дает 1, то n // не //
равно степени 2
{(i& (i-1))
{= (double *) realloc (vect, sizeof (double) *
(i+1));[i] = 0;
}
}* y = (complex*) malloc (sizeof (complex) *
i);(vect, i, y);* y1 = (double*) malloc (sizeof (double) * i);(int j = 0; j
< i; j++)[j] = sqrt (powx (y [j]. re,2) + powx (y [j]. im,2));(y1,
n);(y);(y1);
}if (num == 2)
{* rPCA = (double*) malloc (sizeof
(double));(vect, rPCA, n, 1000);(rPCA, 1000);(rPCA);
}if (num == 3)
{("%s","Please enter the name of
the second file\n");* vect_1 = (double *) malloc (sizeof
(double));(vect_1);* corr = (double *) malloc (sizeof (double) * n);(vect,
vect_1, corr, n);(corr, n);(vect_1);(corr);
}if (num == 4)
{* acorr = (double *) malloc (sizeof (double)
*n);(&vect [0], n, &acorr [0]);(acorr, n);(acorr);
}if (num == 5)= false;if (num == 6)
{= false;= false;
}
}
}0;
}