Влияние функции невязки на адаптивную фильтрацию ЭКГ
Влияние функции невязки на адаптивную фильтрацию ЭКГ
Функция невязки определяет, как адаптивный фильтр корректирует оценку сетевой помехи. В примере рассматриваются невязка по текущему остатку и невязка по разности соседних остатков. На модельном электрокардиосигнале (ЭКС) показано, как эти способы влияют на восстановление формы сигнала и чувствительность фильтра к смещению изолинии.
Настройка отображения
Функция engee.clear() выполняет очистку рабочего пространства:
engee.clear()
Для визуализации результатов доступны два режима отрисовки графиков:
- для статической быстрой отрисовки используется
gr(). - для интерактивной визуализации с возможностью масштабирования и просмотра значений используется
plotlyjs().
Выберите один из режимов, закомментировав второй. В данном примере используется интерактивный режим.
#gr()
plotlyjs()
Установка необходимых библиотек
Для формирования вариабельности частоты сердечных сокращений используется библиотека Random. Перед подключением проверяется её наличие. если библиотека отсутствует, выполняется Pkg.add(). Функция fft из библиотеки EngeeDSP применяется для расчёта амплитудного спектра ЭКС.
let
installed_packages = collect(x.name for (_, x) in Pkg.dependencies() if x.is_direct_dep)
list_packages = ["Random"]
for pack in list_packages
pack in installed_packages || Pkg.add(pack)
end
end
using Random
import EngeeDSP.Functions: fft
Формирование модельного электрокардиосигнала
Используется исходный файл ecg_generator.jl, расположенный рядом со скриптом. Функция формирует ЭКС как сумму несимметричных гауссовых компонент кардиоцикла. Частота сердечных сокращений немного изменяется от цикла к циклу. Фиксированное значение seed обеспечивает воспроизводимость.
Сигнал формируется один раз на всём интервале 60 с. Параметры view_start и view_duration определяют только отображаемый участок и не меняют длительность расчёта.
include(joinpath(@__DIR__, "ecg_generator.jl"))
fs = 1000.0 # Частота дискретизации, Гц
duration = 60.0 # Длительность моделирования, с
hr = 75.0 # Средняя частота сердечных сокращений, уд/мин
hr_sigma = 3.0 # Стандартное отклонение ЧСС, уд/мин
phase0 = 0.16 # Начальное положение внутри кардиоцикла, с
t, ecg = generate_ecg(
fs=fs, # Частота дискретизации, Гц
duration=duration, # Длительность записи, с
hr=hr, # Средняя ЧСС, уд/мин
amplitude=1.0, # Коэффициент масштаба амплитуды
variable_hr=true, # Включение вариабельности ЧСС
hr_sigma=hr_sigma, # Стандартное отклонение ЧСС, уд/мин
phase0=phase0, # Начальное положение внутри кардиоцикла, с
seed=1234 # Зерно генератора для воспроизводимости
)
view_start = 0.0 # Начало отображаемого участка, с
view_duration = 5.0 # Длительность отображаемого участка, с
idx = (t .>= view_start) .& (t .< view_start + view_duration)
plot(t[idx], ecg[idx],
xlabel="Время, с",
ylabel="Амплитуда, мВ",
label="ЭКС",
title="Модельный электрокардиосигнал",
ylims=(-0.5, 1.25),
legend_position=:outertopright
)
На осциллограмме видны повторяющиеся кардиоциклы. Запись начинается вблизи R-зубца, поэтому первые отсчёты полезного сигнала заметно влияют на начальный переходный процесс фильтров.
Моделирование сетевой помехи
К полезному ЭКС добавим гармоническую помеху:
Здесь мВ и Гц. Полученный массив входного сигнала будет общим для обоих фильтров.
f0 = 50.0 # Номинальная частота сетевой помехи, Гц
A_net = 0.5 # Амплитуда сетевой помехи, мВ
d = 0.002 # Шаг адаптивной коррекции, мВ
net = A_net .* sin.(2π .* f0 .* t)
ecg_noisy = ecg .+ net
p = plot(t[idx], ecg[idx],
xlabel="Время, с",
ylabel="Амплитуда, мВ",
label="Исходный ЭКС",
title="ЭКС с сетевой помехой",
ylims=(-1.0, 1.5),
legend_position=:outertopright
)
plot!(p, t[idx], ecg_noisy[idx], label="ЭКС с помехой")
Сетевые колебания накладываются на зубцы ЭКС и особенно хорошо заметны между кардиоциклами. Адаптивный фильтр должен оценить эту гармонику и вычесть её из входа.
Принцип адаптивной фильтрации
Гармонический сигнал можно прогнозировать по двум предыдущим отсчётам. Поэтому предварительная оценка очередного значения помехи рассчитывается как
После прогноза определяется невязка. Сравним два варианта:
- Текущий остаток: .
- Разность остатков: .
В первом случае адаптация реагирует на абсолютное значение остатка, включая полезный ЭКС. Во втором постоянная составляющая остатка сокращается, а медленные изменения ослабляются. Однако быстрые изменения QRS-комплекса влияют на оба варианта. Кроме того, постоянное смещение меняет начальные условия, поэтому разностный вариант не является полностью независимым от него.
В обоих случаях прогноз корректируется по знаку невязки:
Шаг мВ определяет скорость подстройки. Большой шаг ускоряет коррекцию, но усиливает колебания оценки и реакцию на полезный сигнал. Начальные оценки задаются как в исходном примере: , . Поэтому первые два отсчёта выхода равны нулю.
Рекуррентный прогноз и разностная невязка соответствуют формулам (4.3) и (4.4) пособия Фролова и соавторов, с. 49. Коррекция оценки и вычитание помехи описаны там формулами (4.5) и (4.6). В коде сохранён один шаг коррекции на отсчёт, как в исходном демопримере Engee.
Фильтр с невязкой по текущему остатку
Функция adaptive_filter_direct() использует текущий остаток:
Она возвращает очищенный ЭКС, оценку помехи и невязку. Для первых двух отсчётов невязка не определена и обозначается NaN.
# f0 - частота настройки, Гц; d - шаг коррекции, мВ.
function adaptive_filter_direct(x, fs; f0=50.0, d=0.002)
N = length(x)
e = zeros(Float64, N)
y = zeros(Float64, N)
mismatch = fill(NaN, N)
for k in 1:min(2, N)
e[k] = x[k]
y[k] = x[k] - e[k]
end
c = cos(2π * f0 / fs)
for n in 2:N-1
e_pred = 2 * c * e[n] - e[n-1]
mismatch[n+1] = x[n+1] - e_pred
if mismatch[n+1] > 0
e[n+1] = e_pred + d
elseif mismatch[n+1] < 0
e[n+1] = e_pred - d
else
e[n+1] = e_pred
end
y[n+1] = x[n+1] - e[n+1]
end
return (y=y, e=e, mismatch=mismatch)
end
Фильтр с невязкой по разности остатков
Функция adaptive_filter_difference() сравнивает остатки соседних отсчётов:
Прогноз, шаг коррекции, начальные условия и возвращаемые величины совпадают с первым вариантом.
# f0 - частота настройки, Гц; d - шаг коррекции, мВ.
function adaptive_filter_difference(x, fs; f0=50.0, d=0.002)
N = length(x)
e = zeros(Float64, N)
y = zeros(Float64, N)
mismatch = fill(NaN, N)
for k in 1:min(2, N)
e[k] = x[k]
y[k] = x[k] - e[k]
end
c = cos(2π * f0 / fs)
for n in 2:N-1
e_pred = 2 * c * e[n] - e[n-1]
mismatch[n+1] = (x[n+1] - e_pred) - (x[n] - e[n])
if mismatch[n+1] > 0
e[n+1] = e_pred + d
elseif mismatch[n+1] < 0
e[n+1] = e_pred - d
else
e[n+1] = e_pred
end
y[n+1] = x[n+1] - e[n+1]
end
return (y=y, e=e, mismatch=mismatch)
end
Результат фильтрации
Применим оба варианта к одному входному массиву. Частота настройки, шаг и начальные условия совпадают. На графике показаны первые 5 с, включая начальный переходный процесс.
y_direct, e_direct, mismatch_direct = adaptive_filter_direct(ecg_noisy, fs; f0=f0, d=d)
y_difference, e_difference, mismatch_difference = adaptive_filter_difference(ecg_noisy, fs; f0=f0, d=d)
p = plot(t[idx], ecg[idx],
xlabel="Время, с",
ylabel="Амплитуда, мВ",
label="Исходный ЭКС",
title="Сравнение восстановленных ЭКС",
legend_position=:outertopright
)
plot!(p, t[idx], y_direct[idx], label="Текущий остаток")
plot!(p, t[idx], y_difference[idx], label="Разность остатков")
В начале записи оценка помехи содержит полезную составляющую первых отсчётов. По мере адаптации эта ошибка уменьшается. Для просмотра другого участка измените view_start и повторно выполните ячейки, начиная с формирования ЭКС. длительность моделирования останется равной 60 с.
p = plot(t[idx], net[idx],
xlabel="Время, с",
ylabel="Амплитуда, мВ",
label="Истинная помеха",
title="Оценка сетевой помехи",
legend_position=:outertopright
)
plot!(p, t[idx], e_direct[idx], label="Текущий остаток")
plot!(p, t[idx], e_difference[idx], label="Разность остатков")
Сопоставление с известной помехой показывает начальный переходный процесс и отклонения оценок. Для детального просмотра отдельных периодов можно увеличить нужный участок средствами интерактивного графика.
Сравнение функций невязки
Построим временные зависимости невязок двух фильтров. Они позволяют увидеть, как текущий остаток и разность остатков реагируют на участки ЭКС с медленными и быстрыми изменениями. Это сравнение поясняет работу коррекции, а точность восстановления оценивается отдельно по ошибке ЭКС.
p = plot(t[idx], mismatch_direct[idx],
xlabel="Время, с",
ylabel="Невязка, мВ",
label="Текущий остаток",
title="Два способа формирования невязки",
legend_position=:outertopright
)
plot!(p, t[idx], mismatch_difference[idx], label="Разность остатков")
Разностная невязка имеет меньший вклад медленных компонент. Сравнивать её амплитуду с прямой невязкой как показатель качества нельзя: эти величины определены по-разному. Качество оценивается по ошибке восстановления ЭКС.
error_direct = y_direct .- ecg
error_difference = y_difference .- ecg
p = plot(t[idx], error_direct[idx],
xlabel="Время, с",
ylabel="Ошибка, мВ",
label="Текущий остаток",
title="Ошибка восстановления ЭКС",
legend_position=:outertopright
)
plot!(p, t[idx], error_difference[idx], label="Разность остатков")
# Максимум по модулю ошибки после начального переходного процесса.
idx_after_transient = t .>= 2.0 # Интервал оценки ошибки, 2-60 с
println("Без смещения изолинии. Максимальная абсолютная ошибка на интервале 2-60 с, мВ")
println("Текущий остаток: ", round(maximum(abs.(error_direct[idx_after_transient])), digits=5))
println("Разность остатков: ", round(maximum(abs.(error_difference[idx_after_transient])), digits=5))
Без смещения изолинии оба фильтра подавляют сетевые колебания. Максимальные ошибки восстановления оцениваются на интервале 2-60 с, без начального участка переходного процесса. Всплески ошибки около QRS-комплекса отражают реакцию адаптации на быстрые изменения ЭКС.
Амплитудные спектры
Используется функция spect_amp() из исходного примера. Рассчитывается односторонний амплитудный спектр с нормировкой на число отсчётов. Спектры строятся по всем 60 с. ограничение отображения до 5 с относится только к временным графикам.
function spect_amp(x, fs)
N = length(x)
X = fft(x)
K = div(N, 2) + 1
A = abs.(X[1:K]) ./ N
iseven(N) ? (A[2:end-1] .*= 2) : (A[2:end] .*= 2)
f = collect(0:K-1) .* fs ./ N
return f, A
end
f, A_clean = spect_amp(ecg, fs)
_, A_noisy = spect_amp(ecg_noisy, fs)
_, A_direct = spect_amp(y_direct, fs)
_, A_difference = spect_amp(y_difference, fs)
idx_f = f .<= 100.0 # Верхняя граница отображаемого спектра, Гц
p = plot(f[idx_f], A_clean[idx_f],
xlabel="Частота, Гц",
ylabel="Амплитуда, мВ",
label="Исходный ЭКС",
title="Спектры до и после фильтрации",
legend_position=:outertopright
)
plot!(p, f[idx_f], A_noisy[idx_f], label="ЭКС с помехой")
plot!(p, f[idx_f], A_direct[idx_f], label="Текущий остаток")
plot!(p, f[idx_f], A_difference[idx_f], label="Разность остатков")
Уменьшение линии 50 Гц демонстрирует подавление наводки. Однако полезный ЭКС тоже содержит компоненты в этой области, поэтому спектральный график следует рассматривать вместе с ошибкой относительно известного чистого сигнала.
Фильтрация ЭКС со смещением изолинии
Добавим к тому же модельному ЭКС постоянное смещение:
Вход обоих фильтров равен . Параметр offset задаёт постоянное смещение в мВ. Сетевая помеха остаётся прежней. Фильтры должны подавить её и сохранить смещение полезного ЭКС.
offset = 0.3 # Постоянное смещение изолинии ЭКС, мВ
ecg_offset = ecg .+ offset
ecg_offset_noisy = ecg_offset .+ net
offset_direct = adaptive_filter_direct(ecg_offset_noisy, fs; f0=f0, d=d)
offset_difference = adaptive_filter_difference(ecg_offset_noisy, fs; f0=f0, d=d)
p = plot(t[idx], ecg_offset[idx],
xlabel="Время, с",
ylabel="Амплитуда, мВ",
label="Эталонный ЭКС",
title="Смещение изолинии +0,3 мВ",
legend_position=:outertopright
)
plot!(p, t[idx], offset_direct.y[idx], label="Текущий остаток")
plot!(p, t[idx], offset_difference.y[idx], label="Разность остатков")
При смещении ЭКС у фильтра с текущим остатком сохраняются заметные колебания. Разностный вариант лучше сохраняет форму сигнала относительно смещённого эталона. Это связано с тем, что постоянная составляющая сокращается при вычислении разности остатков.
Для проверки результата ниже показана ошибка относительно ecg_offset. Добавленное смещение входит в эталонный сигнал и не считается помехой.
error_offset_direct = offset_direct.y .- ecg_offset
error_offset_difference = offset_difference.y .- ecg_offset
p = plot(t[idx], error_offset_direct[idx],
xlabel="Время, с",
ylabel="Ошибка, мВ",
label="Текущий остаток",
title="Ошибка восстановления ЭКС со смещением",
legend_position=:outertopright
)
plot!(p, t[idx], error_offset_difference[idx], label="Разность остатков")
println("Со смещением изолинии +0,3 мВ. Максимальная абсолютная ошибка на интервале 2-60 с, мВ")
println("Текущий остаток: ", round(maximum(abs.(error_offset_direct[idx_after_transient])), digits=5))
println("Разность остатков: ", round(maximum(abs.(error_offset_difference[idx_after_transient])), digits=5))
При смещении изолинии у прямого варианта сохраняются выраженные колебания ошибки. Разностный вариант слабее реагирует на постоянную составляющую остатка. Максимальные ошибки приведены для одного и того же интервала 2-60 с.
Влияние выбора функции невязки
Графики показывают, что способ формирования невязки меняет характер остаточной ошибки. Текущий остаток учитывает абсолютное значение полезного сигнала и может направлять коррекцию даже при точном прогнозе помехи. Разность остатков сокращает постоянную составляющую, но быстрые изменения QRS-комплекса продолжают влиять на адаптацию.
При смещении изолинии различие между вариантами становится заметнее. Максимальная абсолютная ошибка на интервале 2-60 с характеризует наибольшее отклонение восстановленного ЭКС от эталона после исключения начального участка переходного процесса. Этот показатель не описывает среднюю ошибку и все особенности сохранения формы сигнала.
Полученные результаты показывают влияние двух критериев при заданных параметрах ЭКС и фильтров. Для выбора критерия в другой задаче необходимо учитывать форму полезного сигнала, характер помехи и требования к допустимым искажениям.
Использованные источники
- Фролов А. В., Дик С. К., Лихачевский Д. В., Давыдов М. В., Чечик Н. М. Цифровая обработка биомедицинских сигналов и изображений. Минск: БГУИР, 2016.
- Адаптивная фильтрация сетевой помехи электрокардиосигнала - сообщество Engee.