CC BY
9ISSN G868-5886
НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2GG3, том ІЗ, № 2, c. 57-63
ОРИГИНАЛЬНЫЕ СТАТЬИ
УДК 543:[535+532]+ 621.391.1
©А. А. Евстрапов, А. Л. Буляница, Г. Е. Рудницкая, Б. Г. Беленький,
А. О. Петряков, В. Е. Курочкин
ОСОБЕННОСТИ ПРИМЕНЕНИЯ АЛГОРИТМОВ ЦИФРОВОЙ ФИЛЬТРАЦИИ ЭЛЕКТРОФОРЕГРАММ ПРИ АНАЛИЗЕ ВЕЩЕСТВ НА МИКРОЧИПЕ
В данной работе анализируются особенности применения алгоритмов цифровой фильтрации информативных сигналов, получаемых при анализе веществ с помощью микрофлюидной аналитической системы (МФАС), прототип которой создан в Институте аналитического приборостроения РАН при участии сотрудников Физико-технического института им. А.Ф. Иоффе РАН.
При этом сопутствующие вопросы обработки аналитической информации, связанные с механизмом формирования информативного аналитического сигнала, а также способами и средствами его регистрации остаются за рамками данной работы.
ВВЕДЕНИЕ
В Институте аналитического приборостроения РАН создан и апробирован прототип микрофлюидной аналитической системы (МФАС) для изучения и анализа биологических веществ. В прототипе реализована идеология "lab-on-a-chip" (или Micro Total Analytical System), в соответствии с которой все аналитические действия с микрообъемами пробы происходят на миниатюрном устройстве — микрофлюидном чипе, разделение пробы осуществляется электрофоретическим методом, а регистрация результатов производится с помощью детектора лазер-индуцированной флуоресценции. Составляющий основу МФАС стеклянный микрочип имеет рабочую длину канала порядка 35 мм и трапецеидальное сечение с размерами 30 и 50 мкм при глубине 10 мкм.
Особенностями измерения информативных сигналов являются: во-первых, съем данных
с низкой тактовой частотой, что требует оценивания параметров узких пиков; во-вторых, наличие высокочастотного шума достаточно большой интенсивности.
Применением цифровых фильтров (ЦФ) нижних частот удается выполнить две задачи: существенно подавить высокочастотный шум и расширить пик. Очевидно, что шум влияет на оценивание положения вершины пика и его амплитуды. Кроме того, необходимость оценивания площади, ограниченной пиком, характеризующей количество соответствующего анализируемого компонента, часто требует выделения достаточно широкого пика. Последнее крайне затруднительно при низкой частоте дискретизации. Использование фильтра нижних частот позволяет достичь подавления
высокочастотного шума с одновременной деформацией пика: смещение положения максимума и уменьшение его амплитуды, увеличение площади, ограниченной пиком, расширение границ пика и т.п.
Кроме того, необходимым этапом обработки информативного сигнала является компенсация нулевой (базовой) линии. В большинстве случаев базовая линия может быть аппроксимирована линейным трендом первого или второго (реже нулевого) порядка. При этом может быть осуществлена соответствующая компенсация тренда по Кендаллу—Стьюарту [1]. Анализ применения указанного алгоритма остается за рамками данной работы. Однако следует отметить, что дискретная передаточная характеристика алгоритма компенсации тренда соответствует передаточной характеристике ЦФ верхних частот. То есть комбинированное применение алгоритмов фильтрации высокочастотного шума и подавления базовой линии по существу требует использования полосового ЦФ.
В данной работе ограничимся анализом приме -нения ЦФ нижних частот второго порядка, использование которых не вносит существенных временных задержек и которые легко встраиваются в программу съема данных.
ВЫБОР ЦИФРОВОГО ФИЛЬТРА
За основу взят дискретный цифровой тангенс-ный фильтр Баттерворта второго порядка [2], передаточная характеристика которого имеет вид
H (z) =
kn (1 + z-1)2
і -1 -2
1 -ax z + a z
(1)
Все коэффициенты (1) однозначно определяются относительной частотой среза фильтра (П), т. е. частотой среза (/0), отнесенной к частоте дискретизации (/*): П = /0 / / * . Для сравнения приведем передаточную характеристику синусного фильтра Баттерворта нижних частот второго порядка, имеющую более простую форму, нежели (1),
Н (2) =
к0
1 -а,02 1 +а22 2
Здесь коэффициенты
*2 _ _ 2 а / =&о
2т
(2)
12=1
Интеграл (2) вычисляется на основе теории вычетов в точках, лежащих внутри единичного круга. В нашей ситуации вычеты должны вычисляться в точках г = 0 и г = ^12, где последние — корни знаменателя дроби (1) и лежат внутри еди-
Рис. 1. Амплитудно-частотная характеристика ЦФ нижних частот (Баттерворта, второго порядка)
ничного круга комплексной плоскости ^ Расчеты показывают, что в явном виде величина
кп (3 + а1 -а2)
к = а//а0
выражается как
к=
\<1 4* I 1_Л< 2 ^
также однозначно определяются частотой среза фильтра. Однако синусный фильтр не обеспечивает хорошего подавления шума в том случае, если частота среза не очень близка к нулю — уже при П = 0.10 и выше (см. рис. 1). Также заметим, что оба фильтра — тангенсный и синусный — представляют тип цифровых фильтров с бесконечной импульсной характеристикой (БИХ-фильтры).
Приняв в качестве модели высокочастотного шума модель белого шума [3] с дискретным изображением X(г) = Сто / г , можно оценить величину дисперсии шума на выходе фильтра () формулой
1 | Н (г)Н (г ч) —.
2(1 -а2)
Эта зависимость к от П, построенная по 10 равноотстоящим точкам в диапазоне от 0.01 до 0.40, весьма близка к линейной (выборочный коэффициент корреляции г = 0.9996).
МОДЕЛИ ПИКА. ВЛИЯНИЕ ФИЛЬТРАЦИИ НА ФОРМУ ПИКА
Рассмотрим две идеализированные модели пика: островершинную (линейную) и гладко вершинную (параболическую). Соответственно их уравнения будут (3а), (3б):
модель 1
х[к ] = <
модель 2
х[к ] =
к / т для к = 0,1,..., т,
2 - к / т для к = т + 1, т + 2,...,2т, 0 для к і [0,2т];
(3б)
1 - (1 - к / т)2 для к = 0,1,...,2т,
0 для к і [0,2т].
(3 а)
(3б)
Здесь к — номер такта измерения; параметр т можно считать полушириной пика; высота пика равна единице.
В зависимости от ширины пика (2т) и частоты среза фильтра (П) происходит искажение формы пика, выразившееся в трех эффектах: а) смещение максимума, б) уменьшение амплитуды пика, в) увеличение площади, ограниченной пиком. Все эти эффекты в основном ослабляются по мере увеличения частоты среза фильтра.
Степень уменьшения амплитуды иллюстрируется данными табл. 1 а и 1 б.
Табл. 1а. Зависимость амплитуды от параметров т и П при островершинных пиках
т Исходный П=0.04 П=0.06 П=0.09 П=0.12
6 1 0.722 0.740 0.813 0.877
8 1 0.727 0.794 0.868 0.913
10 1 0.758 0.839 0.899 0.931
15 1 0.838 0.901 0.933 0.954
20 1 0.885 0.927 0.950 0.965
25 1 0.911 0.941 0.960 0.972
Табл. 1б. Зависимость амплитуды от параметров т и П при параболических пиках
т Исходный П=0.04 П=0.06 П=0.09 П=0.12
6 1 1.019 0.933 0.955 0.987
8 1 0.941 0.948 0.986 1.000
10 1 0.935 0.971 0.997 1.000
15 1 0.972 0.998 1.000 1.000
20 1 0.994 1.001 0.999 1.000
25 1 1.000 1.000 1.000 1.000
Табл. 2б. Зависимость площади, ограниченной пиком, от параметров т и П при параболических пиках
т Исходный П=0.04 П=0.06 П=0.09 П=0.12
6 7.94 13.59 10.61 9.22 8.70
8 10.63 14.95 12.67 11.60 11.21
10 13.30 16.81 14.96 14.10 13.77
15 19.98 22.36 21.10 20.52 20.30
20 26.65 28.45 27.50 27.06 26.89
25 33.32 34.77 34.00 33.65 33.52
Табл. 2а. Зависимость площади, ограниченной пиком, от параметров т и П при островершинных пиках
т Исходный П=0.04 П=0.06 П=0.09 П=0.12
6 6 9.05 7.46 6.70 6.42
8 8 10.31 9.09 8.52 8.31
10 10 11.85 10.87 10.42 10.25
15 15 16.23 15.58 15.28 15.17
20 20 20.92 20.44 20.21 20.12
25 25 25.74 25.35 25.17 25.10
Динамика изменения площади под пиком (для обеих моделей) в зависимости от частоты среза фильтра и ширины пика иллюстрируется данными табл. 2а, 2б.
Наиболее сложной является оценка возможного смещения максимума пика. Очевидно, что это смещение должно быть связано с уровнем шума, в том числе и фильтрованного шума. Стационарное (методическое) смещение максимума может быть учтено и заложено в алгоритм оценивания времени удержания соответствующей анализируемой компоненты. Однако это время не должно зависеть от случайных параметров (от параметров шума и ширины пика, которые априорно неизвестны). Заведомо очевидно, что методическая задержка оценивания положения максимума составляет по крайней мере два шага (в соответствии с порядком цифрового фильтра).
ФУНКЦИОНАЛ КАЧЕСТВА. ПОСТАНОВКА ЗАДАЧИ ОПТИМИЗАЦИИ ВЫБОРА ЧАСТОТЫ СРЕЗА ФИЛЬТРА НИЖНИХ ЧАСТОТ
В качестве исходных предпосылок примем: априорно неизвестными предполагаются параметры пика (в частности, его относительная ширина), а также его модель (остро- или гладковершинная). Искажения сигнала, которые могут явиться следствием фильтрации и должны быть включены в функционал качества, были приведены выше: уменьшение амплитуды максимума пика, неучтенное смещение временного положения максимума (и его возможная вариация) и неучтенное увеличение площади, ограниченной пиком. Кроме того, в функционал качества следует включить величину к (см. ранее), которая, безусловно, влияет как на оценку величины площади, так и на величину вариации положения максимума, т. к. характеризует некомпенсированный (нефильтрованный) шум, искажающий информативный сигнал.
Таким образом, функционал качества должен иметь вид
Ф = ^1к + ^2 А + w3^ + w4V, (4)
где wi — весовые коэффициенты; А, ^ — величины, характеризующие соответственно искажение амплитуды и неучтенное искажение площади, и V — величина, интегрально характеризующая смещение положения максимума и вариации этого смещения.
Оценка величин А, ^ и V может быть проведена на основе модельных примеров при отсутствии входного шума и, как следствие, фильтрованного шума. Целью оптимизации будет Ф ^ шт по параметру П.
Зависимость положения максимума фильтрованного сигнала от ширины пика и от модели
Табл. 3а. Зависимость временного положения максимума пика от параметров т и П при островершинных пиках
т П=0.04 П=0.06 П=0.09 П=0.12
6 10 9 9 8
8 13 12 11 10
10 15 14 13 12
15 21 19 18 17
20 26 24 23 22
25 31 29 28 27
Выражение (5) расшифровывается следующим образом:
а) измерение на шаге i с учетом аддитивной случайной помехи может принимать любое значение t в пределах yi ± г ;
б) измерения, выполненные на других шагах (к = 1,..., п; к Ф i), должны быть меньше V,
в) события, описанные в пункте б, независимы и должны выполняться одновременно.
Вероятности соответствующих событий рассчитываются, основываясь на гипотезе равномерного закона распределения помехи:
Табл. 3б. Зависимость временного положения максимума пика от параметров т и П при параболических пиках
т П=0.04 П=0.06 П=0.09 П=0.12
6 9 9 9 8
8 12 12 11 10
10 15 14 13 12
15 21 19 17 17
20 26 24 22 22
25 31 29 27 27
г Ч7} .
Оценим вероятность того, что у{, i = 1,2,...,п будет максимумом. этого события есть
1 У1+г
измерение
Вероятность
Р(£< 2) =
1 для г > г,
(г + г) /(2г) для - г < г < г, 0 для г < -г.
представлена в таблицах 3а, 3б. Выявляется общая тенденция: при достаточно большой частоте среза (П > 0.06) смещение максимума пика в точности соответствует только временной задержке фильтра (два такта), т. к. используется фильтр второго порядка. По мере уменьшения частоты среза величина задержки увеличивается (до 6 тактов). Времена задержек для обеих моделей практически одинаковы. Т. е. этот параметр определяется исключительно частотой среза фильтра. Оценить возможные вариации максимума (и его амплитуды) можно, только задавшись величиной дисперсии входного высокочастотного шума. По существу требуется построение статистической модели фильтрованного шума.
Предлагается следующая статистическая модель: шум является независимой (для измерения в различные моменты времени) равномерно распределенной случайной величиной £ с размахом
Р( уг = шах) = — | П Р(^к < t - Ук Ж (5)
1,— 1
Очевидно, что различная дисперсия входного шума приведет (при фиксированной частоте среза фильтра) к различной дисперсии фильтрованного шума, различному размаху помехи £ и ряду распределений. Полагаем исходный шум имеющим дисперсию 0.03 и 0.10. При фильтрации с П = 0.04 дисперсии фильтрованного шума 7 будут соответственно 0.009 и 0.03; при П = 0.09 — 0.013 и 0.044.
На примере пика с полушириной т = 6 для обеих моделей приведем расчет вероятностей того, что измерение, произведенное на шаге i, будет максимумом. Результаты представлены в табл. 4а, 4б. Обозначения МО и Д даны математическому ожиданию и дисперсии временного положения максимума пика соответственно.
Как следует из данных таблиц:
а) уровень входного шума практически не влияет на оценку положения максимума, при этом возможные вариации этой оценки (дисперсия) существенно возрастают с ростом входного шума;
б) фактор смещения максимума существенно уменьшается с ростом параметра П.
Очевидно, что аналогичным образом могут быть табулированы модельные пики любой заданной ширины и при любом уровне входного шума. Заметим, что смещение максимума пика при достаточно больших П практически (с точностью до одного дополнительного такта) совпадает с задержкой, задаваемой фильтром Баттерворта второго порядка (т.е. двумя тактами). В то же время при малых П вариации параметров А, 5 и V для моделей 1 и 2 различны: вариации амплитуды и временное смещение максимума пика меньше у параболической модели, но возможные вариации
ОСОБЕННОСТИ ПРИМЕНЕНИЯ АЛГОРИТМОВ.. Табл. 4а. Расчет вероятностей нахождения максимума для острого пика с т = 6
; П = 0.04 i П = 0.09
Уi Р(у = тах) У; Р(У; = тах)
а2г = 0.009 о} = 0.030 о) = 0.013 о2г = 0.044
7 0.6138 0 0 7 0.7244 0 0.0206
8 0.6684 0 0.0197 8 0.8077 0.3831 0.4273
9 0.7065 0.0773 0.2222 9 0.8134 0.6169 0.4931
10 0.7222 0.7075 0.4285 10 0.7452 0 0.0590
11 0.7125 0.2152 0.2913 11 0.6215 0 0
12 0.6764 0 0.0382 12 — — —
М О 10.14 10.11 М О 8.62 8.59
Д 0.2735 0.7441 Д 0.2350 0.4011
Табл. 4б. Расчет вероятностей нахождения максимума для параболического пика с т = 6
; П = 0.04 ; П = 0.09
У; Р(У; = тах) У; Р(У = тах)
о} = 0.009 о} = 0.030 о2 = 0.013 о} = 0.044
7 0.9703 0 0.0225 7 0.8943 0 0.0615
8 1.0027 0.0526 0.1900 8 0.9472 0.3476 0.3726
9 1.0186 0.5847 0.3825 9 0.9548 0.6524 0.4531
10 1.0147 0.3626 0.3278 10 0.9088 0 0.1128
11 0.9866 0 0.0761 11 0.8052 0 0
М О 9.31 9.24 М О 8.65 8.62
Д 0.3278 0.8607 Д 0.2268 0.5849
оценки максимума меньше для островершинных пиков (модель 1). Особенно этот эффект выражен для пиков с большой шириной т > 15.
Основные выводы, следующие из анализа зависимостей, представленных в табл. 1а-4б, в значительной степени повторяют сделанные ранее:
— уровень шума не является фактором, смещающим максимум фильтрованного пика;
— по мере роста частоты среза (П) смещение максимума уменьшается и при П = 0.09 с хорошей точностью совпадает с методическим смещением, даваемым фильтром второго порядка (2 такта);
— эффект предыдущего пункта при больших П мало зависит от ширины пика т;
— на величину смещения мало влияет модель пика (остро- или гладковершинная), за исключением случаев очень узких пиков (менее 8); в этом случае при П = 0.04 смещение МО максимума при гладковершинной модели существенно меньше;
— при малых частотах среза фильтра само смещение достаточно мало связано (незначитель-
но возрастает) с шириной пика (при островершинной модели — от 4.14 при т = 6 до 5.77 при т = = 25); коэффициент корреляции между смещением максимума и шириной пика т (т = 6, 8, 10, 15, 20 и 25) равен 0.79, что свидетельствует о высокой корреляционной связи, но характер зависимости существенно нелинейный;
— по мере роста ширины пика т и увеличения амплитуды входного нефильтрованного шума растет дисперсия оценки положения максимума (для обеих моделей пика);
— как правило, для параболической (гладко-вершинной) модели пика дисперсия временного положения максимума выше.
МОДИФИКАЦИЯ ФУНКЦИОНАЛА КАЧЕСТВА. РЕШЕНИЕ ЗАДАЧИ ОПТИМИЗАЦИИ
Как следует из данных таблиц 1а-4б, на оценивание параметров А, Б и V параметр к существенно не влияет. Его величина косвенно влияет на пара-
метр V и непосредственно (но менее существенно в количественном плане) на оценки амплитуды пика и площади, ограниченной пиком. Последнее влияние, очевидно, связано с искажением каждого измерения вследствие аддитивной помехи. Благодаря центрированности помехи этот эффект, особенно при оценивании площади под широким пиком, будет минимальным.
Однако, несмотря на отсутствие значимого влияния величины к на параметры выявленных пиков, эта величина существенным образом влияет на обнаружение, выявление и локализацию пиков, поскольку обнаружение пика представляет собой выявление разладки в последовательности измерений (в форме скачкообразного изменения порядка или параметра положения тренда). Принятие решения может базироваться на критериях типа Фишера, чувствительных к масштабу случайной помехи.
При измерениях полагаем априорно неизвестными:
а) амплитуду входного шума (в диапазоне 0.03-
0.10);
б) относительную ширину (т) и характер (модель) пика.
Таким образом, в функционал качества должны быть внесены усредненные как по ширине, так и по шуму характеристики.
Пример поиска оптимальной частоты среза фильтра
Выбираем функционал качества в форме, аналогичной (4). Важнейшей проблемой является подбор весовых функций wi или так называемое "обучение" функционала. Это обучение выполняется на основании экспертных оценок, носящих субъективный характер. Как правило, при этом в качестве оптимума находится не граничная, а какая-либо промежуточная точка. Исключение может составить случай доминирования какого-либо одного слагаемого в выражении (4). Например, при необходимости уменьшения дисперсии фильтрованного шума без оценивания остальных характеристик фильтрованного информативного сигнала оптимум находится на минимально реализуемой частоте среза фильтра.
На рис. 2, а, б представлены фрагменты информативного сигнала, фильтрованные цифровым фильтром (1) с различными частотами среза.
Анализ кривых на рис. 2, а, б позволяет сделать следующие основные выводы:
— достаточно адекватной представляется модель пика малой ширины т;
— форма первого и, возможно, второго пика более соответствует островершинной модели; для третьего и четвертого пиков более адекватна гладковершинная модель; таким образом, выбор модели пика, скорее всего, должен быть не детермини-
стическим, а стохастическим (вероятностным);
— в обоих случаях качество фильтрации шума на нейтральной линии (отсчеты 400-550), по-видимому, достаточное, чтобы исключить ложное выявление пиков;
— не имея каких-либо дополнительных составляющих целевой функции (4), помимо названных выше, представляется, что результат фильтрации при П = 0.09 (рис. 2, б) должен быть ближе к оптимуму (т. е. величина функционала качества должна быть меньше, чем для случая рис. 2, а);
— последнее утверждение базируется на экспертной оценке авторов статьи и в силу этого субъективно. Однако, при согласии с ним, можно наложить определенные ограничения на весовые функции, входящие в функционал, что позволит частично "обучить" его.
Повторение подобного анализа для модельных пиков при большем числе частот среза позволит предложить достаточно обоснованные оценки весовых функций wi.
а
б
Рис. 2. Фильтрация фрагмента информативного сигнала.
а — ФНЧ с П = 0.04; б — ФНЧ с П = 0.09
ЗАКЛЮЧЕНИЕ
Данная работа не претендует на окончательное решение проблемы по созданию универсального оптимального алгоритма обработки информативных сигналов МФАС и, возможно, информативных сигналов приборов капиллярного электрофореза, высокоэффективной жидкостной хроматографии и аналогичных. В работе описаны теоретические основы для выбора оптимального алгоритма цифровой фильтрации, для чего предложен критерий оптимизации и многокомпонентная целевая функция. Использование разработанных алгоритмов в макете МФАС позволило существенно улучшить соотношение сигнал/шум и выявить неявные информационные пики. Однако дополнение алгоритмов обработки еще и процедурой компенсации базовой линии может существенно расширить возможности по выбору оптимального (по заданному критерию) алгоритма компенсации как высокочастотного шума, так и низкочастотного сигнала базовой линии.
Кроме того, на основе параметров прибора можно априорно более точно установить: а) вероятностно-статистическую модель пика (по принципу: с вероятностью Р — островершинная, 1-Р — гладковершинная); б) диапазон входного шума (исходя из параметров измерительных устройств, в частности разрядности АЦП и т.д.); в) априорно сузить диапазон возможных значений относительной ширины пика (т), прежде всего исходя из частоты дискретизации измерений и номенклатуры анализируемых веществ (суммарные электрофоретические и электроосмотические подвижности и т. д.).
Адаптация критерия оптимальной фильтрации информативных сигналов МФАС к более узким классам анализируемых веществ или/и конкрет-
ным приборам также является предметом дальнейших исследований авторского коллектива.
Благодарности
Авторы благодарят сотрудников Центра "Физика наногетероструктур" Физико-технического института им. А.Ф. Иоффе РАН П.С. Копьева, В.Л. Суханова, О.Ф. Позднякова и В.В. Филимонова за неоценимую помощь в работе, выразившуюся в создании стеклянных планарных микрочипов.
Работа выполнялась при поддержке МНТП "Вакцины нового поколения и медицинские диагностические системы будущего" в рамках НИР № ВНП-02-06 "Новые принципы детекции и разработка на их основе приборов для автоматизации лабораторно-диагностических методов исследования".
СПИСОК ЛИТЕРАТУРЫ
1. Кендалл М., Стьюарт А. Многомерный статистический анализ и временные ряды. М: Наука, 1976. С.505-509.
2. Верешкин А.Е., Катковник В.Я. Линейные цифровые фильтры и их реализация. М: Сов. Радио. 1973. 152 с.
3. Первозванский А.А. Курс теории автоматического управления: Учеб. пособ. М: Наука, 1986. 616 с.
Институт аналитического приборостроения РАН, Санкт-Петербург
Материал поступил в редакцию 11.03.2003.
CHARACTERISTIC FEATURES OF DIGITAL SIGNAL FILTERING ALGORITHMS AS APPLIED TO ELECTROPHORESIS ON A MICROCHIP
A. A. Evstrapov, A. L. Bulyanitsa, G. E. Rudnitskaya,
B. G. Belenkii, A. O. Petryakov, V. E. Kurochkin
Institute for Analytical Instrumentation RAS, Saint-Petersburg
The paper deals with characteristic features of digital signal filtering algorithms applied to substance analysis in a microfluidic analytical system (MFAS) developed at the Institute for Analytical Instrumentation RAS in association with Ioffe Physico-Technical Institute RAS. The accompanying problems of analytical data processing related to analytical signal forming mechanisms as well as tools and methods of signal detection are not consi dered.
