Smoothing of chromatographic signals by their approximation in the basis of Chebyshev-Hermite functions

Full Text

Введение

Одной из задач обработки хроматографических сигналов является сглаживание или подавление помехи. На практике сглаживание чаще всего реализуют различными цифровыми фильтрами. Проблемой цифровой фильтрации является тот факт, что частотный спектр пиков хроматограммы, как правило, частично совпадает со спектром помехи и вместе с подавлением помехи происходит ослабление полезной составляющей, приходится либо снижать сглаживающие свойства фильтра, либо пытаться восстановить потерянную при сглаживании информацию [1]. Фильтр Савицкого – Голея решает проблему ослабления полезной составляющей, но осуществляет сглаживание менее эффективно, причем степень сглаживания зависит от подхода к расчету параметров фильтра [2]. Фильтр Калмана имеет лучшие сглаживающие свойства, чем фильтр Савицкого – Голея, но подразумевает существенную зависимость от формы шума, выходной сигнал фильтра может иметь фазовые и прочие нелинейные искажения [3, 4]. Сплайн-фильтры эффективно подавляют помеху с учетом гладкости полезного сигнала и при сравнении превосходят фильтр Калмана, но имеют значительную зависимость погрешности от метода расчета коэффициентов сплайна [4–6]. Фильтр Винера также эффективно сглаживает пики хроматограммы, но вызывает осцилляции, обусловленные краевыми эффектами; трудностью являются одновременно хорошее сглаживание и низкие осцилляции [3]. Таким образом, актуальным является вопрос разработки нового математического метода сглаживания хроматографических сигналов.

В настоящей работе для создания такого метода предлагается воспользоваться аппроксимацией сигналов в базисе функций Чебышева – Эрмита (за рубежом просто Эрмита). При этом из коэффициентов аппроксимации в данном базисе возможно восстановить по другим, предварительно рассчитанным базисам не только сам сигнал [7], но и его производную различных порядков [7, 8], вейвлет-преобразование [9, 10] и т. п. Функции Чебышева – Эрмита находят широкое применение в науке и технике. Разложение по данным функциям используется для аппроксимации двухмерных [11] и одномерных [12] сигналов и аналитически заданных функций, обработки файлов мультимедиа форматов, в частности для выделения границ зашумленного изображения [13], биометрической идентификации по изображению радужки глаза человека [14], фильтрации и поиска изображения по базе данных, индексации говорящего по записи голоса [15], извлечения полезной медицинской информации из изображений компьютерной томографии [16] и распознавания речи [17]. Помимо этого, данные функции используются для построения банка фильтров [18], диагностики состояния оборудования [19] и здоровья человека [20], обнаружения полезных составляющих в сигнале [21], а также для решения дифференциальных уравнений [22, 23].

Аппроксимация сигналов в базисе функций Чебышева – Эрмита

Функции Чебышева – Эрмита определяются следующим образом:

${\phi }_n\left(x\right)=\frac{1}{{\alpha }_n}\cdot e^{-\frac{x^2}{2}}\cdot H_n\left(x\right)$, (1)

где ${\text{α}}_n=\sqrt{2^{n}n!\sqrt{\text{π}}}$ – нормирующая константа; $H_n\left(x\right)$ – полином Эрмита $n$-го порядка;

$H_n\left(x\right)={\left(-1\right)}^n{e}^{x^2}\frac{d^n{e}^{-x^2}}{d{x}^n}.$ (2)

Набор функций ${\text{φ}}_n\left(x\right)$ (1) задает базис для кодирования и декодирования сигналов. Представление сигнала $f\left(x\right)$ с помощью $N+1$ функций разложения (алгоритм кодирования) заключается в нахождении коэффициентов $c_n$. В силу локализации функций Чебышева – Эрмита нахождение данных коэффициентов можно свести к формуле:

$c_n=\underset{-{\tau }_n}{\overset{{\tau }_n}{\int }}f\left(x\right)\cdot {\text{φ}}_n\left(x\right)dx, n=\mathrm{0,1,}\dots ,N$, (3)

где $\left[-{\tau }_n,{\tau }_n\right]$ – отрезок, на котором локализована функция ${\text{φ}}_n\left(x\right)$.

Алгоритм декодирования с помощью $N+1$ коэффициентов, в результате которого восстанавливается аппроксимированный сигнал, выражается формулой:

$\widehat{f}\left(x\right)=\sum _{n=0}^N{c}_n\cdot {\phi }_n\left(x\right).$ (4)

Описанный алгоритм кодирования-декодирования обладает сглаживающим свойством [18] и позволяет получить не только сглаженный сигнал, но также и его сглаженную производную k порядка при наличии соответствующего базиса. Формируя различные базисы перехода от коэффициентов разложения (4) к различным преобразованиям, можно получить универсальный аппарат для быстрого получения также коэффициентов непрерывного вейвлет-преобразования, устранения аппаратной функции прибора и т. д. Для восстановления необходимого преобразования сигнала требуется заменить в (4) ${\phi }_n\left(x\right)$ на соответствующее преобразование данной функции и вычислить полученную формулу относительно рассчитанных для исходного сигнала коэффициентов разложения (4) [7].

Одной из отличительных черт функций Чебышева – Эрмита является их локализация на временной оси. Данное обстоятельство может создать трудности при кодировании сигналов большой длительности – даже при простой их конфигурации (6 и менее точек перегиба) для кодирования в рассматриваемом базисе потребуются функции высоких порядков, вычисление которых затруднительно; трудность вычисления только нарастает с ростом порядка функции.

Очевидным решением является масштабирование функций Чебышева – Эрмита на некоторый коэффициент $\gamma$. Помимо этого наилучших результатов кодирования можно достичь, если положение максимума в сигнале будет совпадать с началом координат базисных функций, что обусловлено наибольшим перекрытием областей локализации исследуемого сигнала и базисных функций. Следовательно, требуется ввести и коэффициент сдвига для данных функций, в настоящей работе обозначенный как x₀. Тем самым формируется новый базис функций ${\phi }_n\left(\frac{x-x_0}{\gamma }\right)$. Для восстановления нормировки базиса был рассчитан коэффициент $\sqrt{1/\gamma }$, на который умножается и результат кодирования, и результат декодирования. Тем самым можно вынести множитель $1/\gamma$ в выражение для вычисления коэффициентов кодирования, которое при этом следует записать:

$c_n=\frac{1}{\gamma }\underset{-{\tau }_n}{\overset{{\tau }_n}{\int }}f\left(x\right){\phi }_n\left(\frac{x-x_0}{\gamma }\right)dx$. (5)

С учетом описанных масштабных коэффициентов был разработан алгоритм аппроксимации сигналов в базисе функций Чебышева – Эрмита (рис. 1). Помимо вычисления коэффициентов масштаба и сдвига в ходе данного алгоритма происходит деление сигнала на фрагменты, что снижает погрешность аппроксимации [24].

 

Рис. 1. Алгоритм аппроксимации

 

Для решения различных задач первичной обработки в хроматографии используются сглаженные производные хроматографических сигналов. Как было сказано выше, при аппроксимации сигналов в базисе функций Чебышева – Эрмита можно аппроксимировать также преобразования этого сигнала; для получения производной требуется вычислить производную функций Чебышева – Эрмита, получив выражение

$\frac{d^k{\phi }_n\left(\frac{x-x_0}{\gamma }\right)}{d{x}^k}=\frac{1}{{\gamma }^k}\sum _{m=0}^{\min \left(k,n\right)}\left[C_k^m{\left(-1\right)}^{k-m}\sqrt{\frac{2^{m}n!}{\left(n-m\right)!}}H{e}_{k-m}\left(\frac{x-x_0}{\gamma }\right){\phi }_{n-m}\left(\frac{x-x_0}{\gamma }\right)\right],$ (6)

которое задает базис производной k порядка n-й функции Чебышева – Эрмита с учетом масштаба $\gamma$ и сдвига $x_0$ и может быть использовано для восстановления k-й сглаженной производной исходного сигнала по рассчитанным заранее коэффициентам разложения в базисе Чебышева – Эрмита, если вынести за знак суммы нормировочный коэффициент:

$\frac{d^k\widehat{f}\left(x\right)}{d{x}^k}=\frac{1}{{\gamma }^k}\sum _{n=0}^N{c}_n\cdot \frac{d^k{\phi }_n\left(\frac{x-x_0}{\gamma }\right)}{d{\left(\frac{x-x_0}{\gamma }\right)}^k},$ (7)

где $\frac{d^k\widehat{f}\left(x\right)}{d{x}^k}$ – оценка k-й сглаженной производной сигнала $f\left(x\right); c_n$ – коэффициенты кодирования (5) сигнала $f\left(x\right)$ в базисе функций Чебышева – Эрмита.

Подавление случайной помехи

Аппроксимация сигналов в базисе функций Чебышева – Эрмита обладает сглаживающими свойствами, что позволяет устранять случайную помеху. На рис. 2 представлен пример сглаживания реального хроматографического сигнала c предварительным его делением на фрагменты по алгоритму рис. 1. Начальный и конечный фрагменты классифицируются как шум, поэтому на рис. 2 данные области обнулены.

 

Рис. 2. Пример сглаживания хроматографического сигнала

 

Для оценки эффективности устранения помехи относительно описанных в первой главе методов было произведено сравнение результата декодирования сглаженного сигнала по коэффициентам базиса функций Чебышева – Эрмита со сглаженным сигналом, полученным при использовании наиболее распространенного на практике метода сглаживания – фильтра скользящего среднего:

$L_k\left(N\right)=\frac{1}{N}\sum _{k=-\mathrm{0,5}\left(N-1\right)}^{\mathrm{0,5}\left(N-1\right)}\varsigma \left(i,k\right)S_{i-k},$ (8)

где

$\varsigma \left(i,k\right)=\left\{\begin{array}{l}\mathrm{1,}\text{ если }\left(i-k>0\right)\vee \left(i-k\le M\right),\\ \text{иначе }\mathrm{0,}\end{array}\right.$ (9)

N – порядок фильтра (только нечетное число); $\Delta x$ – интервал дискретизации сигнала.

При исследованиях был задан сигнал S, полученный дискретизацией пика, заданного выражением

$f\left(t\right)=A{e}^{-\frac{{\left(t-\mu \right)}^2}{2{\sigma }^2}},$ (10)

где A – амплитуда пика; µ – время удерживания пика; σ – среднеквадратичная ширина пика [25].

На полученный сигнал накладывалась сгенерированная псевдослучайная помеха с нормальным законом распределения и различным СКО помехи σ. При разных σ достигалось разное отношение сигнал/шум (ОСШ):

$\text{ОСШ}=20\lg \left(\frac{S_{\text{СКЗ}}}{n_{\text{СКЗ}}}\right),$ (11)

где $S_{\text{СКЗ}}$ – среднеквадратическое значение сигнала; $n_{\text{СКЗ}}$ – среднеквадратическое значение помехи.

Оценка погрешности сглаженного сигнала производилась по выражению

$\delta =\max \left(\frac{{\widehat{S}}_i-S_i}{\max \left(S_i\right)}\right)\cdot 100\%,$ (12)

где $\widehat{S}$ – сглаженный сигнал; $S$ – исходный сигнал до наложения помехи; ${\widehat{S}}_i-S_i$ – разность i-х значений сглаженного и исходного сигналов.

На рис. 2 данное сравнение представлено в виде графиков, где а – сглаживание с применением цифрового фильтра N порядка, б – путем восстановления сигнала по коэффициентам кодирования N порядка (максимальный порядок функции) в базисе функций Чебышева – Эрмита.

По рис. 3, б можно видеть, что при больших значениях ОСШ цифровой фильтр малого порядка обеспечивает низкую погрешность $\left(\delta <1\%\right)$, но со снижением ОСШ повышение порядка не дает снижения погрешности, поскольку, как было сказано в первой главе, вместе с помехой происходит ослабление полезной составляющей. Для сглаживания в базисе функций Чебышева – Эрмита видна иная картина – по достижении минимально необходимого порядка (в данном случае – $N=8$) обеспечивается $\delta <1\%$ и со снижением ОСШ погрешность растет значительно медленнее, чем при использовании цифрового фильтра.

 

Рис. 3. Сравнение сглаживания с цифровым фильтром

 

Поскольку помимо самих сигналов сглаживание требуется (и в большей степени) при вычислении их производных, для того же S была оценена погрешность вычисления производных первого и второго порядка относительно численно определенной производной [26] по сглаженному цифровым фильтром сигналу. Для оценки погрешности использовалось выражение (12) с заменой сигнала на его производную (как в случае исходного, так и в случае сглаженного сигнала). Результаты представлены на рис. 4, где а и в – соответственно первая и вторая производные, полученные численно по сглаженному цифровому сигналу, б и г – восстановленные по базису (7) сглаженные производные.

По данному рисунку можно видеть, что применение цифрового фильтра не позволяет обеспечить погрешность $\delta <50\%$ для обеих производных даже при больших значениях ОСШ, а с его уменьшением ситуация меняется только в сторону ухудшения. Увеличение порядка фильтра не дает существенных результатов, позволяя вернуться лишь к начальному значению погрешности $\delta \approx 50\%$. Применение сглаживания в базисе функций Чебышева – Эрмита дает иные результаты: при достижении минимально необходимого порядка (в данном случае – $N=8$) обеспечивается $\delta <1\%$ для обеих производных, но при снижении ОСШ погрешность возрастает до $\delta \approx 10\%$.

 

Рис. 4. Сравнение погрешности восстановления производной

 

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

$\delta =\frac{\left|ex\left[{\widehat{S}}^{\left(n\right)}\right]-ex\left[S^{\left(n\right)}\right]\right|}{ex\left[S^{\left(n\right)}\right]}\cdot 100\%,$

где $ex\left[\cdot \right]$ – функция для вычисления положения экстремума; ${\widehat{S}}^{\left(n\right)}$ – восстановленная по функциям Чебышева – Эрмита производная n порядка; $S^{\left(n\right)}$ – полученная при цифровой фильтрации производная n порядка.

В качестве экстремумов рассматриваются наиболее применяемые на практике точки: положение прохода через ноль для первой производной и положение минимума второй производной [1]. Результаты представлены на рис. 5.

 

Рис. 5. Сравнение погрешности восстановления экстремумов производных

 

Можно видеть, что применение цифровой фильтрации позволяет сохранить положение экстремумов при достижении определенного N (при котором согласно рис. 5 возрастает погрешность аппроксимации производной) с $\delta <5\%$ для первой производной и с $\delta <10\%$ для второй. Восстановление сглаженной производной в базисе функций Чебышева – Эрмита при достижении N=8 обеспечивает для большинства значений ОСШ близкую к нулю погрешность, которая на малых его значениях растет до $\delta <3\%$ для первой производной и $\delta <5\%$ для второй. Таким образом, подавление помехи в используемом базисе обеспечивает снижение погрешности относительно традиционного метода.

Погрешность восстановления одиночных пиков

Для оценки погрешности представленного метода были аппроксимированы пики реальных хроматограмм. Хроматограммы были извлечены из программы Хроматэк Аналитик 3.1. Все они были получены на хроматографе Кристалл-5000 с пламенно-ионизационным детектором.

Из полученных хроматограмм были извлечены 9 наиболее характерных пиков, конфигурация которых наиболее часто встречается в реальных хроматограммах (рис. 6). Многие пики на практике близки по форме к гауссиане, но с небольшой разницей в крутизне фронтов; при выборе пиков для обработки были собраны наиболее часто встречаемые асимметричные конфигурации, которые, как правило, являются одиночными. Для оценки зависимости погрешности отличия пика от симметричной гауссианы (как наиболее простой для обработки модели) был использован метод моментов, на основе которого вычислялись коэффициенты асимметрии ${\gamma }_1$ и эксцесса ${\gamma }_2$ [27]:

$m_0=\Delta x\sum _{i=1}^M{S}_i, m_1=\frac{\Delta x^2}{m_0}\sum _{i=1}^M{S}_i\left(i-1\right), m_n=\frac{\Delta x}{m_0}\sum _{i=1}^M{S}_i{\left[\left(i-1\right)\Delta x-m_1\right]}^n, n>\mathrm{1,}$

${\gamma }_1=\frac{m_3}{m_2^{3/2}}, {\gamma }_2=\frac{m_4}{m_2^2}-3.$

 

Рис. 6. Исследуемые реальные пики

 

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

Результаты приведены в таблице. Закон распределения помехи в пустых областях исходных сигналов совпадает с законом распределения остатка. Математическое ожидание равно нулю во всех случаях, поэтому в таблице не приводится. Почти для всех рассматриваемых сигналов помеха описывается нормальным законом, для других – законом Коши или логистическим законом, которые отличны от нормального плавностью фронтов, но приближенно являются сопоставимыми с нормальным. Гистограммы распределения помехи для некоторых сигналов приведены на рис. 7.

 

Разделение симметричных пиков

Пик	γ1	γ2	СКО шума сигнала	СКО остатка	Закон распределения шума сигнала
№ 1	0,10	0,42	1,72∙10-2	1,47∙10-2	Нормальный
№ 2	–0,69	–3,59	1,54∙10-2	1,76∙10-2	Нормальный
№ 3	0,79	1,02	1,50∙10-2	1,59∙10-2	Нормальный
№ 4	0,96	0,41	1,97∙10-2	1,57∙10-2	Нормальный
№ 5	0,97	0,78	1,71∙10-2	1,44∙10-2	Нормальный
№ 6	1,34	3,83	0,91∙10-1	1,14∙10-1	Нормальный
№ 7	1,53	3,38	1,39∙10-2	1,79∙10-2	Нормальный
№ 8	1,86	3,85	1,55∙10-2	3,18	Коши
№ 9	2,37	8,78	4,70∙10-2	6,35∙10-1	Логистический

 

Рис. 7. Гистограммы распределения помехи: а, б – для пика № 8, распределение Коши; в, г – для остатка хроматограммы, нормальное распределение; д, е – для пика № 9, логистическое распределение

 

Выводы

Анализируя результаты, приведенные в таблице, можно видеть, что СКО шума совпадает порядком с СКО остатка во всех случаях, кроме пиков 8 и 9. Поскольку согласно рис. 6 помеха имеет существенно меньшую амплитуду по сравнению с полезным сигналом, можно положить, что совпадение порядков СКО шума и остатка говорит о низкой погрешности сглаживания. Следовательно, конфигурации пиков, подобных номерам 8 и 9, не могут быть аппроксимированы с достаточно низкой погрешностью, что является ограничением представленного метода. Соответственно предельным значением асимметрии для аппроксимированных одиночных пиков можно считать наибольшую асимметрию в таблице, не считая параметров пиков 9 и 10, ${\gamma }_{1\max }=\mathrm{1,53}$, а влияния γ₂ на качество аппроксимации не выявлено. С учетом сравнения качества сглаживания сигнала при использовании предложенного аппроксимационного метода и цифрового фильтра можно заключить, что аппроксимационный подход позволяет совместить высокую степень сглаживания и низкую погрешность.

About the authors

Andrey V. Bochkarev

Samara State Technical University

Email: bochkarevpubliciit@gmail.com
ORCID iD: 0000-0003-1615-5659

Postgraduate Student

Russian Federation, 244, Molodogvardeyskaya st., Samara, 443100

Raukhat T. Sayfullin

Samara State Technical University

Author for correspondence.
Email: bochkarevpubliciit@gmail.com
SPIN-code: 4594-4229

Professor

Russian Federation, 244, Molodogvardeyskaya st., Samara, 443100

References

Supplementary files