Коррекция фазовых искажений

Данная статья была опубликована в журнале Вестник МГТУ. Сер. Приборостроение. 1993. № 4. с.40-53
В статье рассматриваются редко учитываемые фазовые искажения, вносимые в сигнал аналоговыми цепями предварительных усилителей. На примере двух конкретных приборов: реографа и фотоплетизмографа, проанализированы источники возникновения таких искажений, показаны возможные последствия фазовых искажений для результатов цифровой обработки сигналов, предложены пути их устранения с помощью преобразования Фурье.

Рекомендуем также обратить внимание на статьи:
Непрерывная амплитудно-фазо-частотная цифровая коррекция аппаратных искажений реографического сигнала в режиме ON-LINE
Средства и методы измерения параметров кровообращения. Реокардиомониторные системы PDF (~670 Kb)
Breath waves filtration in bioimpedance cardiography.

УДК 534.6
Беляев К.Р., Морозов А.А.

Коррекция фазовых искажений и обработка биомедицинских сигналов.

Изложены проблемы графических методов обработки биомедицинских диагностических сигналов, связанные с фазо-частотными характеристиками первичных измерительных преобразователей. Показаны практические пути реализации процедур фазовой коррекции сигналов, их области применения, приведен обзор вычислительных методов для решения задач восстановления сигналов. Показаны результаты моделирования реальных медицинских диагностических преобразователей. Работа может быть полезна специалистам в области создания автоматизированных систем обработки медико-биологической информации.
Создание неинвазивных диагностических медицинских приборов и систем является перспективным и активно развивающимся направлением. Для получения стабильных, воспроизводимых результатов измерения физических параметров биообъекта необходимо стандартизировать условия измерения, и устранить систематические ошибки, связанные с амплитудно-фазочастотными характеристиками (АФЧХ) измерительных преобразователей биомедицинских датчиков. В рассчитываемых параметрах состояния пациента эти искажения могут давать, как показали исследования, погрешность до 50% (например, при контурном анализе импедансной реовазограммы). С амплитудно-фазочастотные искажения можно устранять апостериорной обработкой полученного сигнала. Такая обработка предполагает большой объем вычислений, поэтому ее использование возможно в рамках диагностического комплекса, включающего ЭВМ.

Низкочастотная диагностическая аппаратура имеет в своем составе пассивные или активные цепи фильтрации сигналов, параметры которых выбираются исходя из противоречивых требований обеспечения помехоустойчивости, необходимой избирательности и минимально возможных искажений формы сигналов. При этом в лучшем случае нормируются амплитудно-частотные характеристики (АЧХ) фильтрующих звеньев, а фазовые характеристики в результате модификации аппаратуры претерпевают существенные изменения. В данной работе сделана попытка обощения и анализа существующих методов устранения фазочастотных искажений и предложены практические способы решения задач восстановления сигналов.

Благодаря успехам, достигнутым в области радиоэлектроники, вычислительной техники существенно изменилась и продолжает модифицироваться аппаратура для клинических исследований и регистрации, а также методы анализа сигналов, при этом значительно расширились границы применения используемых методик. В частности, в последнее время из-за распространения ВИЧ-инфекции, значительно возрос интерес к неинвазивным методам исследования гемодинамики. Эти методы основаны, главным образом, на регистрации пульсовых колебаний давления в сосудах и кровенаполнения тканей и находят широкое применение при изучении внутрисердечной гемодинамики, сократительной функции миокарда, параметров артериального давления, состояния сосудистой системы и регионарного кровообращения. Уже давно стали традиционными такие методы клинического исследования гемодинамики как сфигмография, флебография, плетизмография, осциллография; быстро развиваются сравнительно новые методы - импедансная реоплетизмография, фотоэлектроплетизмография, тахоосциллография, ангиотензография, исследование производных кривых пульса и др. [1]. Исследование кривых пульса позволяет получить представление о ряде параметров кровообращения. Его целью в зависимости от задач является изучение:

  • сердечной деятельности и центральной гемодинамики;
  • функционального состояния сосудистой системы, в том числе параметров артериального и венозного давления, упругого состояния (тонуса) артерий, соотношений между сердечным выбросом и периферическим сопротивлением;
  • регионарного кровообращения, в том числе проходимости периферических сосудов, величин перфузионного давления, объемной и линейной скорости кровотока.
  • Для осуществления этих целей предложено множество методик, среди которых необходимо в первую очередь различать методики, основанные на регистрации перемещений стенки сосуда (артерии или вены) и изменений объема тканей.

    Для оценки параметров состояния пациента по пульсовым кривым необходимо быть уверенным, что полученные результаты достоверны, сопоставимы и могут быть воспроизведены многократно. Таким образом естественно возникает требование стандартизации условий получения данных для анализа.

    Следует отметить существование ряда однотипных медицинских приборов, например, реографов, предназначенных для регистрации изменения электрического сопротивления, которое связывают с пульсовым кровенаполнением, в которых, в целях фильтрации сигнала и повышения помехоустойчивости применяются различные технические решения. Так как нормированию подвергают лишь АЧХ устройства, а фазовые характеристики не оговариваются, то полученные в идентичных условиях сигналы при использовании реографов с различными фазочастотными характеристиками (ФЧХ) будут различаться, что повлечет изменение вычисляемых показателей, например при проведении контурного анализа (таблица). Алгоритмы расчета контурных показателей взяты из работы [2].

    Изменение реографических гемодинамических показателей после проведения фазовой коррекции

    Показатель до коррекции после отличие(%)
    Амплитуда (Ом) 0,022 0,023 5
    Время изгнания (с) 0,18 0,18 0
    Дикротический индекс (%) 15,8 24,5 55
    Диастолический индекс (%) 28,7 35,6 24
    Временной показатель сосудистого тонуса 0,88 0,75 15
    Амплитудный показатель сосудистого тонуса 0,41 0,39 5
    Максимальная скорость быстрого наполнения (Ом/с) 0,13 0,15 15
    Средняя скорость медленного наполнения (Ом/с) 0,17 0,18 6
    Показатель эластичности (%) 33,8 34,8 3
    Ударный объем кровотока (по Кубичеку, мл) 2,53 2,44 4
    Ударный объем кровотока (по Найбору, мл) 1,79 1,63 9
    Минутный объем кровотока (по Кубичеку, мл/мин) 189,7 185,5 2
    Минутный объем кровотока (по Найбору, мл/мин) 134,1 123,7 8

    Примечание: даже довольно хорошая ФЧХ реографа все же искажает некоторые индексы более чем на 50%

    При решении проблемы фазочастотных искажений существуют 2 пути:

    1. зафиксировать какую-либо схему реографа, что фактически останавливает их техническое совершенствование;
    2. использовать алгоритм приведения выходной кривой (или рассчитываемых параметров) к какой-либо стандартной ФЧХ прибора, при этом естественным является использование в качестве стандартной ФЧХ тождественно равной нулю для всех частот (т.е. прибор не искажает фазу сигнала). Поскольку в оптике хорошо разработаны алгоритмы математического восстановления искаженных изображений при известных характеристиках передающих систем, то при использовании данного подхода мы можем опираться на уже готовый математический аппарат. Целью данной работы является попытка обобщить результаты применения эффективных методов восстановления сигналов, регистрируемых при помощи различной аппаратуры, к неискаженному виду, при этом акцент будет делаться в первую очередь на компенсацию именно фазовых искажений, хотя все рассматриваемые методы справедливы и для компенсации амплитудных.
    Для решения задач такого рода предлагается использовать методы спектрального анализа. Исходными данными при этом являются получаемые на выходе прибора сигналы, АЧХ и ФЧХ прибора и допущения о линейности преобразования сигнала прибором и известности полосы частот, в которой лежат информативные гармоники спектра сигнала. При этом обрабатываемый сигнал подвергается дискретизации с достаточным темпом (для реографических сигналов, например, 100 Гц) и квантованию с разрядностью, обеспечивающей требуемую точность представления сигнала (в медицинских задачах обычно 12 двоичных разрядов). Для фазовой коррекции сигнала необходимо определение входного, т.е. неискаженного, сигнала, а применение дискретного преобразования Фурье позволяет совместно с основной решать задачу дополнительной цифровой фильтрации диагностических сигналов.

    Сравнительно низкочастотные диагностические сигналы, тем не менее, приводят к необходимости обрабатывать большие объемы данных, что заставляет принимать во внимание время выполнения программ и применять эффективные алгоритмы дискретного преобразования Фурье.

    Алгоритмы коррекции фазовых искажений

    Существуют два принципиально различающихся способа устранения искажений, вносимых реальным регистрирующим прибором в наблюдаемый сигнал. Первый состоит в том, что при помощи искусных технических решений совершенствуют конструкцию реального прибора, однако по мере улучшения характеристик неизбежно возрастают сложность и стоимость аппаратуры. К тому же иногда характеристики сигнала (например, близкая по частоте помеха) не позволяют улучшать передающие свойства прибора.

    Другим способом устранения искажений является апостериорная обработка сигнала с целью уменьшения влияния реальных характеристик прибора (сведение их к "идеальному" виду). Вопрос о редукции к идеальному прибору, не искажающему полезный входной сигнал был поставлен еще Рэлеем в 1871 г. в связи с некоторыми задачами спектроскопии и состоял в следующем: можно ли после опыта так обработать функцию, наблюдаемую на выходе некоторого реального прибора, чтобы полностью восстановить сигнал, поступивший на его вход. В такой постановке задача относится к классу некорректно поставленных, и в общем случае не может быть решена, однако на практике при соблюдении некоторых требований (ограниченность спектра рассматриваемого сигнала по частоте, отсутствие в рассматриваемом диапазоне частот нулей в передаточной функции прибора и т.д.) она может быть решена с достаточной точностью и иметь практическое значение.

    Различные методы подобного восстановления сигнала хорошо развиты в области обработки изображений (см., например, работу [3]), однако в других областях в настоящее время не получили широкого распространения из-за необходимости цифровой обработки сигнала и большого объема требуемых вычислений (в частности, нам неизвестны медицинские приборы, реализующие такую коррекцию). Тем не менее существующий уровень компьютерной техники позволяет реализовать рассматриваемый алгоритм за приемлемое время (так обработка 1024 отсчетов на IBM PC/AT 80386 без сопроцессора занимает около 5 секунд и может быть ускорена за счет ряда специальных приемов). К тому же, как показывает анализ литературы, большинство исследователей уделяют внимание коррекции только амплитудных, но не фазовых искажений сигналов. В дальнейшем будет показана важность коррекции именно фазовых характеристик сигнала.

    Как правило, реальные приборы с достаточной степенью точности являются линейными, и сигнал f(t) на их выходе описывается интегральным уравнением:
     

    (1)

    где h(t) - импульсная переходная характеристика прибора, s(t) - сигнал, поступающий на вход прибора. Это уравнение является уравнением свертки, и может быть записано как
     

    (2)

    где F(w), H(w), S(w) - фурье-образы функций f(t), h(t) и s(t), т.е. соответственно спектр отклика, передаточная функция системы и спектр входного сигнала. Из уравнения (2) следует, что по крайней мере формально, можно записать
     

    (3)

    откуда обратным преобразованием Фурье находим:
     

    (4)

    Формула (4) справедлива, если f(t)  L2(-,), h(t)  L1(-,), F(w)/H(w)  L2(-,), и сохраняет силу, если H(w) обращается в нуль на множестве нулевой меры.

    На практике применение формулы (4) сопряжено с рядом сложностей:

    1. В действительности H(w) отлична от нуля лишь в конечном интервале частот. Радиотехнические приборы разрабатываются таким образом, что H(w)  0 вне частотного диапазона сигнала, поэтому естественной становится коррекция в диапазоне частот, ограниченном частотными характеристиками сигнала и зануление всех частот, где спектральная плотность мощности шума превышает сигнал. Если нуль частотной характеристики попадает в полосу частот сигнала, то некоторые способы разрешения данной проблемы приведены, например, в работе [2].

    2. Присутствие в сигнале аддитивной помехи. Выходной сигнал прибора f(t) известен всегда приближенно, и содержит в себе случайную помеху v(t):

    f(t) = fт(t) + v(t),

    тогда согласно выражению (3) спектр восстановленного сигнала:

    S(w) = Sт(w) + V(w)/H(w),

    где V(w) - спектр помехи. Спектральная плотность случайной ошибки восстановления (дисперсия решения задачи) равна [3]:
     

    (5)

    где Rv(w)= |V(w)|2 - спектральная плотность мощности помехи. Реальная помеха содержит компоненту белого шума, и ее спектральная мощность стремится к конечному пределу при . Поэтому интеграл (5) может не сходиться либо быть недопустимо большим. В этом случае можно рекомендовать ограничение частотного диапазона коррекции с подавлением ненужных частот спектра.

    3. Погрешности, связанные с квантованием уровня сигнала, округлением при машинном представлении действительных чисел, и т.д. Оценка ошибок восстановления, обусловленных данными погрешностями рассматривается в литературе (например, [3]) и должна в каждом случае оцениваться отдельно. В наших задачах сколько-нибудь заметного пpоявления этих факторов замечено не было.

    4.Граничные эффекты. Наиболее неприятным фактом при коррекции искажений сигналов являются краевые эффекты. На практике сигнал известен не на всей временной оси, а в диапазоне (-Т, Т) (рис. 1,а). Поэтому для формулы фильтрации в пространственной области (1) возникает вопрос о доопределении сигнала за границами рассматриваемого интервала. В простейшем случае считают, что за границами (-Т, Т) сигнал равен своему математическому ожиданию (рис.1,б). Также используется различная экстраполяция (напр. считают, что сигнал вне (-Т, Т) равен своему значению на границах интервала).


    Рис 1. Варианты доопределения сигналов за границами выборки

    При цифровой фильтрации (особенно при использовании аппарата преобразований Фурье), часто исходный сигнал считают периодически продолженным с периодом, равным длине исходной последовательности. Так как начало и окончание записи сигнала расположены в произвольные (независимые от сигнала) моменты времени, то длина исследуемого интервала не равна целому числу периодов сигнала, что приводит к несовпадению амплитудного значения (и значения производных амплитуды) сигнала в начальной и конечной точках и, как следствие, к появлению в спектре высокочастотных гармоник, а после фильтрации - к возникновению осцилляций на краях последовательности (эффект Гиббса). Некоторые исследователи рекомендуют в этом случае использование линейной аппроксимации между начальной и конечной точками (рис. 1,в, длина дополнения должна быть не менее удвоенной длительности импульсной характеристики фильтра) или четного продолжения исходной последовательности (рис. 1,г). Проведенные численные эксперименты показали, что хотя эти способы уменьшают высокочастотные выбросы спектра, однако сильно искажают его низкочастотную область. Более удачные результаты при использовании дискретного преобразования Фурье (ДПФ) получаются при использовании так называемых "весовых окон". Обычно они применяются при оценке спектра по малой выборке: предварительно сигнал умножают на быстро спадающую функцию, например, на кривую Гаусса:
     

    (6)

    где 2T - длительность выборки, | t | < T, а - параметр (наилучшие результаты получаются при 2.5 < а < 3), затем проводят Фурье-преобразование для оценки спектра. Другие весовые окна рассматриваются, например, в работе [4]. Их применимость в нашем случае определяют два фактора: скорость спада боковых лепестков и малые вычислительные погрешностями. Окно Ханна  дает лучшие результаты в смысле отсутствия боковых лепестков, но на границах его значения близки к нулю, поэтому после фильтрации часть информации на концах интервала будет потеряна.

    Исходный сигнал можно рассматривать как бесконечный, умноженный на прямоугольное весовое окно, что в частотной области эквивалентно свертке с Фурье-образом этого окна, функцией sinc(x) = sin(x)/x. В результате истинный спектр как бы "размывается" по частоте (боковые лепестки имеют уровень -13.3 дБ), что приводит

    • к неограниченности полученного спектра (даже если в качестве исходного сигнала взят sin(t));
    • к тому, что гармоника вносит ненулевой вклад во весь спектр сигнала (т.е. их невозможно разделить при фильтрации).
    Предварительное умножение на функцию Гаусса приводит к тому, что свертка производится с ее Фурье-образом, а не с sinc(x), таким образом, за счет некоторого уширения самой спектральной полосы, сильно подавляются ее боковые лепестки (до уровня -45..-55 дБ в зависимости от параметра а), и спектр более точно аппроксимирует спектр реального бесконечного сигнала. На рис.2,а - представлен вычисленный спектр синусоидального сигнала с использованием прямоугольного окна (сплошная линия) и окна Гаусса (теоретически, спектр должен был бы иметь вид дельта-функции). После фильтрации, разумеется, надо вернуть сигнал к первоначальному виду, т.е. разделить на функцию Гаусса. Здесь может возникнуть единственная тонкость, связанная с округлением чисел при машинной обработке.






    Рис. 2 Спектр мощности выборки синусоидального сигнала (сплошная линия) и спектр того же сигнала с использованием окна Гаусса (пунктирная линия) (а);
    видно расширение спектральной полосы и подавление боковых лепестков; исходная выборка (сплошная линия), она же после фильтрации (пунктирная), она же после фильтрации с применением окна Гаусса (штрих-пунктирная) (б)

    Использование прямоугольного окна при фильтрации вносит заметные искажения на протяжении примерно трех периодов сигнала с каждой стороны, в то время как окно Гаусса позволяет уменьшить влияние краев выборки до 0.5 периода На рис.2,б представлен синусоидальный сигнал (сплошная линия), он же после проведения узкополосной фильтрации с применением окна Гаусса (штрих-пунктирная линия) и без применения окна Гаусса (штриховая линия). Вторым моментом, связанным с граничными эффектами, является то, что импульсная передаточная характеристика корректирующего фильтра должна быть конечной во времени, и ее длина должна быть много меньше . На практике это означает невозможность реализации прямоугольной АЧХ фильтра. Поскольку именно такую характеристику считают "идеальной", можно рекомендовать следующий алгоритм построения фильтра. В частотной области определяются амплитудно-фазочастотные характеристики требуемого фильтра, затем при помощи обратного фурье-преобразования вычисляется его импульсная характеристика, которая обрезается до нужной длины (чем короче - тем меньше влияние границ выборки, но и тем ниже порядок получающегося фильтра), и, для повышения гладкости АФЧХ, умножается на функцию Гаусса (6). После этого, вычисляя прямое фурье-преобразование, получим сглаженные АЧХ и ФЧХ нашего фильтра. В работе [5] приведен интересный способ фильтрации, когда помимо восстановления сигнала производится усиление его гармоник пропорционально соотношению сигнал/шум для каждой частоты:

    где H(w) и H(w)* - передаточная функция прибора, и ее комплексно сопряженная величина, Rs(w) и Rv(w) - спектральные плотности мощности входного сигнала и шума соответственно.

    Этот способ предполагает знание спектральных характеристик шума, и, вообще говоря, искажает спектр входного сигнала.

    Вычислительные алгоритмы. Поскольку при обработке сигналов приходится искать компромисс между точностью и скоростью работы математических алгоритмов, то и мы в своем изложении будем уделять основное внимание именно этим двум факторам.

    Коррекция исходного сигнала есть не что иное как фильтрация с передаточной функцией 1/H(w). Синтез подобного фильтра аналоговыми методами в общем случае невозможен, то единственным общим способом его реализации является апостериорная цифровая обработка.

    Цифровой сигнал представляет собой значения реального аналогового сигнала, взятые через равные промежутки времени. Необходимую частоту дискретизации (время между отсчетами) определяют из теоремы Котельникова (она должна быть больше 2Fмакс, где Fмакс - максимальная присутствующая в спектре сигнала частота). На практике используют частоту, превышающую частоту Котельникова (Найквиста) в 2-5 раз. Фильтрация в цифровой области имеет свои особенности, преимущества и недостатки, и сходна с аналоговой фильтрацией. Простейшим и часто используемым методом цифровой фильтрации является прямое вычисление свертки входного сигнала с импульсной передаточной характеристикой (ИПХ) фильтра по формуле (для нерекурсивного фильтра):

    которая представляет собой дискретный аналог формулы (1), здесь Nипх - количество отсчетов импульсной передаточной характеристики. Для вычислений по этой формуле требуется произвести Nc Nипх умножений и столько же сложений (Nc - число отсчетов сигнала), что при длинной ИПХ занимает больше времени, чем реализация такого же фильтра в частотной области.

    Согласно формуле (2) свертка сигналов во временной области аналогична перемножению их спектров. Аналогом преобразования Фурье для дискретных сигналов является дискретное преобразование Фурье (ДПФ):
     

    (7)

    Его главное отличие от непрерывного преобразования Фурье - цикличность, т.е. и сигнал, и его спектр являются периодическими с периодом N. Наиболее употребительные свойства ДПФ:

      1. A*N-S  -- ДПФ комплексно-сопряженного сигнала а;
      2. AS = A*N-S  -- ДПФ действительного сигнала;
      3. AS = AN-S  -- ДПФ четно-симметричного сигнала;
      4. AS = -AN-S  -- ДПФ нечетно-симметричного сигала;
    При непосредственных вычислениях по формуле (7) требуется произвести N2 умножений и сложений комплексных чисел, поэтому на практике пользуются быстрым преобразованием Фурье (БПФ), в котором за счет некоторого дополнительного расхода машинной памяти и использования рекурсии число операций уменьшено до N logN, что при больших выборках дает существенный выигрыш во времени. Существуют различные алгоритмы БПФ, рассматриваемые в специальной литературе (например, [3]), суть которых в том, что исходная последовательность делится на несколько меньших частей, для каждой из которых вычисляется ДПФ (время вычисления которого меньше, т.к. N12+N22 < (N1+N2)2), а затем по найденным спектрам вычисляется спектр исходной последовательности.

    Кроме алгоритма БПФ, двукратного выигрыша по времени выполнения можно достичь, применяя совмещенные алгоритмы ДПФ. Для последовательностей вещественных чисел коэффициенты Фурье с номерами k и N-k являются комплексно-сопряженными числами (свойство 2). Поэтому можно либо преобразовывать две последовательности одновременно, либо разбить исходную последовательность на две, одновременно преобразовать их, а потом пересчитать результаты на всю последовательность в целом. Используя первый способ, из последовательностей Ak и Bk создают последовательность Ck, такую, что: Ck = Ak + i Bk и вычисляют ее ДПФ (Fr). ДПФ исходных последовательностей находят по соотношениям:

    учитывая, что AN-r = Ar, BN-r = Br.

    При втором способе последовательность A длины 2N разбивается на две подпоследовательности: четных элементов (Аr_even) и нечетных (Аr_odd). Далее вычисления ведутся как и в первом способе, а затем по формуле

    находят ДПФ исходной последовательности. Если в исходной последовательности или в ее ДПФ заведомо много нулей (например, некоторые гармоники в процессе фильтрации должны быть обнулены), то дополнительной экономии времени можно достичь просто их не вычисляя.

    Поскольку обычно реальные сигналы вводятся в ЭВМ через АЦП, т.е. являются целыми числами, то естественно для них вести обработку в целочисленной области (которая в 3-10 раз по скорости превышает обработку действительных чисел). Однако в этом случае встает вопрос о точности проводимых преобразований, что ограничивает применимость этих методов для коррекции сигналов. Здесь не приведены численные эксперименты, а показанные ниже методы указаны для полноты обзора, поскольку в литературе встречаются редко.

    Наиболее просто использование целых чисел реализуется в методе прямой свертки. Найденные коэффициенты просто домножаются на некоторое число (чтобы получить необходимую точность), а в процессе свертки полученное отфильтрованное значение делится на то же число.

    В работах [3, 6] приводится метод квантованного БПФ. Суть его в том, что значения синуса и косинуса квантуются, т.е. принимают дискретные значения, например -1, 0, 1 при трех уровнях и -1, -1/2, 0, 1/2, 1 при пяти уровнях (вообще говоря, это разложение сигнала по другим базисным функциям). В этом случае умножения либо вообще отсутствуют, либо заменяются сдвигами. Для фильтрации изображений этот метод дает вполне приемлемые результаты.

    Существует еще ряд базисов (функции Уолша, Хаара, и т.д.), разложение по которым привлекает внимание тем, что операции можно вести с целыми числами, однако мы не исследовали их применимость для цифровой коррекции сигналов.

    Обобщение. После проведенных исследований мы остановились на следующей схеме вычислений.

    1. Предварительно измеряем или вычисляем передаточную функцию прибора H(w). Вычисляем передаточная функция восстанавливающего фильтра 1/H(w) (если предполагается прямая свертка во временной области, то затем вычисляем импульсную характеристику (обратное ДПФ от передаточной характеристики), с учетом приведенных выше замечаний.
    2. Исходный сигнал умножаем на функцию окна.
    3. Проводим фильтрацию-восстановление (во временной или в частотной области).
    4. Отфильтрованный сигнал делим на функцию окна.
    Вопрос о том в какой области проводить фильтрацию-восстановление - во временной или в частотной, решается исходя из требуемой точности и скорости обработки. В частотной области можно достичь наивысшей точности при восстановлении и наиболее крутых фронтов фильтров при приемлемом времени вычислений (как уже упоминалось - примерно 5 с на IBM PC/AT для двух последовательностей в 1024 отсчета при одновременной обработке), однако неудобство состоит в необходимости сначала ввести всю последовательность, а затем ее обработать, к тому же при использовании БПФ ограничения накладываются на число точек в последовательности и появляются довольно жесткие требования к объему машинной памяти. Во временной области три последних недостатка отсутствуют, но сильно возрастает время вычислений, которое можно сократить за счет выбора короткой ИПХ (ориентировочно меньше 0.1N, где N - длина исходной последовательности) и проведения операций с целыми числами, однако при этом сильно проигрываем в точности.

    Для более тонкого понимания эффектов, возникающих при цифровой фильтрации-восстановлении сигналов, можно рекомендовать изучение работ [7, 8].

    Пример практического использования алгоритмов. На рис.3 представлен фрагмент схемы серийно выпускаемого реографа 2РГ. Амплитудно-модулированный полезный сигнал поступает через трансформатор TV1 на вход двухполупериодного выпрямителя выполненного на операционном усилителе DA1. Сглаженное напряжение усиливается инструментальным дифференциальным усилителем на прецизионных операционных уилителях DA2 и DA3 с цепями балансировки и регулятором уровня выходного сгнала. Каскад на второй половине сдвоенного операционного усилителя DA3 выполняет функции фильтра нижних частот с частотой среза 30 Гц (считается что информативный спектр реографических сигналов локализован в полосе частот от 0,5 Гц до 30 Гц). Выходной сигнал этого каскада представляет собой реографический сигнал, который далее поступает на скорректированный дифференциатор (поз. DA4) и пропускается через дополнительный активный фильтр нижних частот второго порядка (поз. DA5). В результате формируется сигнал "Диф. рео", по которому при анализе вычисляют такие показатели, как время изгнания, более точно определяется ударный объем крови.
     








    Рис.3 Фрагмент принципиальной электрической схемы серийно выпускаемого реографа 2РГ.

    В области низких частот частотная характеристика реографа определяется наличием единственного пассивного RC-фильтра высоких частот установленного после детектора сигнала, котороый определяет нижнюю границу полосы пропускания реографа 0,5 Гц. Именно от параметров фильтров высоких частот в основном зависит характер поведения фазовой характеристики в полосе пропускания и связанные с ней искажения формы сигнала. Отметим, что для неискаженной передачи сигналов должно выполняться условие линейности фазовой характеристики цепей прохождения сигнала или постоянства группового времени задержки при равенстве фазы k на нулевой частоте (k - целое число). Для фильтров нижних частот эти условия достаточно легко выполнить с приемлемой точностью применяя, например, активные фильтры Бесселя. Для фильтров высоких частот обычно стараются найти компромисс между величиной фазовых искажений сигнала и необходимой степенью подавления низкочастотных гармоник в составе сигнала, сдвигая частоту среза фильтра высоких частот в низкочастотную область. При этом учитывают еще один критерий качества фильтра - время установления. Для реографических сигналов существенное влияние оказывают артефакты дыхания с типичной частотой 0,3 Гц для которых желательно обеспечить достаточное ослабление, в этом случае сдвиг частоты среза фильтра высоких частот в низкочастотную область применять нельзя.
     




    Рис.4 Схема замещения цепей преобразования сигнала реографа 2РГ



    Рис.5 Результаты моделирования цепей преобразования сигнала реографа 2РГ

    На рис.4 представлена схема замещения цепей преобразования сигнала реографа 2РГ (канал "Рео"). При анализе характеристик применяли систему Micro Cap III. Рис.5 содержит результаты моделирования цепей преобразования, которые полностью подтверждены непосредственными измерениями. По приведенным графикам можно сделать следующие выводы:

    1. Полоса пропускания реографа 2РГ лежит в пределах от 0,17 Гц до 32 Гц по уровню -3 дБ (ширина полосы 7,5 октав).
    2. Частота среза удалена от нижней граничной частоты информативного спектра сигнала на 2 октавы.
    3. В основной полосе частот информативного спектра от 1 Гц до 10 Гц групповое время задержки не является постоянной величиной, что определяет наличие фазовых искажений сигнала.
    4. Разность группового времени задержки на границах полосы пропускания составляет 488 мс, при 527 мс на низкочастотной и 45 мс на высокочастотной границах полосы пропускания.
    5. Разность сдвига фаз на границах полосы пропускния составляет 292 град.
    6. Амплитудно - частотная характеристика в полосе пропускания постоянна в пределах 0..-3 дБ.
    Обобщая эти выводы можно рекомендовать использование метода фазовой коррекции сигналов, принимая амплитудно - частотную характеристику постоянной для информативной области спектра реографических сигналов.

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

    Принимая во внимание частотный состав реографического сигнала (типичный вид которого показан на рис.6) и полосу пропускания реографа , определяется Fд - частота дискретизации, которая должна быть больше, по крайней мере, 60 Гц. Для упрощения аппаратно-программных решений, процедуры восстановления сигнала и уменьшения ошибок дискретизации, частота дискретизции выбрана равной 100 Гц.
     








    Рис.6 Типичный реографический сигнал (голень)

    Длительность сигнала Тс определяется из условий близости текущего спектра к идеальному спектру бесконечно длинной выборки и ограничена способностью алгоритмов обработки преобразовывать большие объемы данных и временем выполнения программы. Из условия близости спектра, ограниченной во времени выборки сигнала, к идеальному спектру необходимо Тс выбрать такой, чтобы в ней помещалось 3-4 периода дыхательных волн, присутствующих в составе реографического сигнала и представляющих собой наиболее низкочастотные гармоники спектра сигнала.

    Из физиологических диапазонов частоты дыхания следует что длительность сигнала Тс должна быть не менее 9..12 секунд, при Fд=100 Гц получается порядка 1000 отсчетов сигнала. Такие объемы данных, подлежащих преобразованию, требуют применения эффективных методов фурье-анализа, например быстрого преобразования фурье с прореживанием по времени. Такие алгоритмы накладывают специфические ограничения на количество отсчетов сигнала N, которое должно выбираться как N=2k , где k=1,2,3,... - целое число. Отсюда видно, что выбрав Тс=10 сек, получаем k=10 и N=1024 отсчета для Fд=100 Гц, что удовлетворяет перечисленным выше условиям. В выборке такой длительности в среднем находятся десять кардиоциклов и три периода дыхания.
     







    Рис.7. Фрагмент принципиальной электрической схемы фотоплетизмографа

    На рис.7 представлен фрагмент принципиальной электрической схемы фотоплетизмографа, представляющий цепи усиления и фильтрации сигнала. Для получения сигналов пульсового кровенаполнения используется диодная оптопара с открытым оптическим каналом АОД111А в 14-выводном металло-керамическом корпусе с планарным расположением выводов, состоящая из излучающего диода GaAlAs, работающего в инфракрасном диапазоне, изготовленного по эпитаксиальной технологии, и двух кремниевых эпитаксиальных p-i-n фотоприемников. Для регистрации сигнала задействован один фотоприемник, причем проведенный анализ способов включения фотоприемника указывает на целесообразность использования фотодиодного режима (т.е. при обратном смещении фотодиода ). Выбор именно этой серийно выпускаемой оптопары обусловлен относительной простотой изготовления датчиков на ее основе, технологичностью и приемлимыми массогабаритными характеристиками конструкции в целом и удобством использования при исследованиях периферической гемодинамики.

    ИК-фотодиод используется в фоторезисторном включении, составляя одно плечо измерительного моста. Операционный усилитель DA1 является преобразователем тока в напряжение. Второй каскад выполнен на микросхеме DA3 и реализует высокочастотную фильтрацию сигнала с частотой среза 0.6 Гц и крутизной склона АЧХ 12 дБ/окт. Коэффициент передачи сигнала в полосе пропускания близок к единице. Второй каскад (DA4) используется в активном фильтре нижних частот с частотой среза 12 Гц при единичной передаче в полосе пропускания. С целью подавления электромагнитной наводки от промышленной сети используется двойной Т-образный мост позволяющий получить дополнительное подавление на частоте 50 Гц не менее 30 дБ. На выходе устройства установлен линейный усилитель на мс. DA5 с коэффициентом усиления 40 дБ. При нормальных условиях амплитуда выходных сигналов лежит в пределах -5В..+5В.

    Преимуществами данного включения являются высокая линейность преобразования, возможность получения достаточно больших по амплитуде полезных сигналов, низкий суммарный шум (тракта приемник - усилитель). Основным недостатком данного режима является сравнительно большие величины сопротивления входных цепей на которых развивается значительная ЭДС электромагнитной помехи.

    Так, при длине соединительного кабеля один метр уровень сетевых наводок относительно уровня полезного сигнала составлял 15..20 дБ , а в таких условиях необходимо применение высококачественных помехоподавляющих фильтров сложной структуры. Типичный спектр фотоплетизмографического сигнала находится в пределах 0,7 - 15 Гц, поэтому синтезируемые структуры фильтров должны обеспечивать для удовлетворительного подавления электромагнитной наводки промышленной частоты затухание 40-50 дБ к частоте 50 Гц.

    Длительные наблюдения за видом полезного сигнала пульсового кровенаполнения, результаты спектрального анализа и многочисленные вариации параметров фильтрующих звеньев позволяют сделать вывод что информативные составляющие сигнала находятся в полосе частот примерно от 0,7 Гц (основная частота при ЧСС=40 уд/мин) до 8..12 Гц; высшие гармоники определяются свойствами сосудистой системы и придают сигналу характерную для реоплетизмографических данных форму, т.е. фотоплетизмографический сигнал качественно подобен сигналу реоплетизмографическому.

    Предыдущие замечания позволяют сформировать основные требования к выбору типа фильтров для подавления помех, представленных главным образом электромагнитными наводками на высокоомные входные цепи и артефактами движения и вибрации. Помехи от вибрации, лежащие в диапазоне частот полезного сигнала, неотделимы принципиально от полезного сигнала и лишь при наличии некоторой дополнительной информации о помехе, на уровне применения специальных алгоритмов фильтрации можно ослабить ее влияние на результат измерения. По схемотехнической реализации фильтров можно рекомендовать конфигурацию ИНУН или фильтры Саллена-Кея, а также индивидуально настраиваемые фильтры в виде двойного Т-образного моста. Очевидной становится противоречивость требовния эффективного подавления помех и проблемы фазовых искажений формы сигналов, что достаточно хорошо иллюстрирует практическую ценность процедур фазовой коррекции сигналов.

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

    • широким диапазоном активной области, что позволяет уменьшать величину напряжения питания;
    • достаточно высоким выходным сопротивлением, что позволяет получить от фотоприемника почти максимально возможный сигнал;
    • малые шумы в сравнении с аналогичным сопротивлением большой величины;
    • управляемость, что дает возможность организовать достаточно эффективные обратные следящие, частотно-зависимые связи;
    • простота реализации и малые размеры.
    С помощью следящей связи можно скомпенсировать нежелательный квазипостоянный фототок от засветки который изменяет свою величину при изменении условий измерения в несколько раз и инерционную температурную составляющую фототока а одновременно с этим добиться большого усиления полезного сигнала. Такая связь делает поведение всего устройства адаптивным к условиям измерения. При использовании указанных методов уже на выходе первого усилителя (повторителя) удается получить сигналы пульса амплитудой в размахе до 1В и более. Следует отметить еще одно положительное следствие описанных принципов - возможность применения однополярного источника питания; а в случае применения низковольтных ОУ и КМОП-цифровых схем работоспособность измерителя сохраняется при Uпит от 3 до 18В.
     





    Рис.8 Схема замещения цепей преобразования сигнала фотоплетизмографа
    Рис.9 Результаты моделирования цепей преобразования сигнала фотоплетизмографа

    На рис.8 представлена схема замещения усилительного канала фотоплетизмографа для цепей усиления и фильтрации сигнала. При анализе характеристик применялась система Micro Cap III. Рис.9 содержит результаты моделирования канала усиления. По приведенным графикам можно сделать следующие выводы:

    1. Полоса пропускания фотоплетизмографа лежит в пределах от 0,17 Гц до 32 Гц по уровню -3 дБ (ширина полосы 7,5 октав).
    2. Частота среза удалена от нижней граничной частоты информативного спектра сигнала менее чем на октаву.
    3. В основной полосе частот информативного спектра - от 1 Гц до 10 Гц групповое время задержки не является постоянной величиной, что определяет наличие фазовых искажений сигнала.
    4. Разность группового времени задержки на границах полосы пропускания составляет 538 мс, при 574 мс на низкочастотной и 36 мс на высокочастотной границах полосы пропускания.
    5. Разность сдвига фаз на границах полосы пропускния составляет 234 град.
    6. Амплитудно - частотная характеристика в полосе пропускания постоянна в пределах 0..-3 дБ.
    Обобщая эти выводы так же можно рекомендовать использование процедуры фазовой коррекции сигналов, принимая амплитудно - частотную характеристику постоянной для информативной области спектра фотоплетизмографических сигналов.

    Что касается вопроса повышения помехоустойчивости, то здесь можно использовать совокупность следующих методов :

    • тщательное экранирование входных цепей и выносного датчика;
    • уменьшение длины соединительных проводов;
    • питание излучающего диода от генератора стабильного тока;
    • надежная развязка цепей приемника от помех по питанию;
    • расположение предварительного усилителя в едином корпусе с датчиком пульса.
    Последний способ представляется наиболее эффективным, т.к. предварительный буферный усилитель можно расположить (в виде заказной гибридной МС) в одном корпусе с оптопарой. Кроме того, использование оптопары ИК-диапазона существенно уменьшает влияние посторонних паразитных засветок. Еще одним эффективным способом борьбы с наводками является высокочастотная модуляция излучения с последующим детектированием полезного сигнала, но применение этого метода существенно усложняет схему.

    Описанные методы и результаты частично используются в действующем образце компьютерной системы индивидуальной электромагнитной терапии и активной диагностики и в системе для неинвазивной дифференциальной диагностики кровообращения конечности, разработанных на кафедре "Биомедицинские технические системы и устройства" МГТУ им.Н.Э.Баумана


    Литература.

    1. Палеев Н.Р., Каевицер И.М. Атлас гемодинамических исследований в клинике внутренних болезней (бескровные методы).- М.: Медицина, 1975 г.- 240 стр.
    2. Гуревич М.И. Импедансная реоплетизмогрфия.- Киев: Наукова думка, 1982 г. - 186 с.
    3. Ярославский Л.П. Цифровая обработка сигналов в оптике и голографии.- М.: Радио и связь, 1987 г.- 296 с.
    4. Марпл С.Л. Цифровой спектральный анализ и его приложения.- М.: Мир, 1990 г. - 584 с.
    5. Василенко Г.И. Теория восстановления сигналов. - М.: Советское радио, 1979 г. - 272 с.
    6. Ярославский Л.П. Введение в цифровую обработку изображений.- М.: Советское радио, 1979 г.- 312 с.
    7. Хемминг Р.В. Цифровые фильтры. - М.: Недра, 1987 г.- 224с.
    8. Жуков А.И. Метод Фурье в вычислительной математике.- М.Ж Наука, 1992.- 176 с.
    Comments