Адаптивная фильтрация сетевой помехи электрокардиосигнала
Адаптивная фильтрация сетевой помехи электрокардиосигнала
В демонстрационном примере рассматривается подавление гармонической сетевой помехи частотой 50 Гц в модельном электрокардиосигнале. Адаптивный алгоритм формирует текущую оценку периодической помехи и вычитает её из входного сигнала. Отдельно анализируются временная реализация ЭКС и изменение его амплитудного спектра после обработки.
Настройка отображения
Функция engee.clear() выполняет очистку рабочего пространства:
engee.clear()
Для визуализации результатов доступны два режима отрисовки графиков:
- для статической быстрой отрисовки используется
gr(); - для интерактивной визуализации с возможностью масштабирования и просмотра значений используется
plotlyjs().
Выберите один из режимов, закомментировав второй. В данном примере используется интерактивный режим.
#gr()
plotlyjs()
Установка необходимых библиотек
Для формирования вариабельности частоты сердечных сокращений используется библиотека Random. Функция 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
Формирование модельного электрокардиосигнала
Для исследования используется модельный ЭКС с управляемой частотой сердечных сокращений. Форма кардиоцикла задаётся суммой несимметричных гауссовых компонент P, Q, R1, R2, S, ST и T. Подход является развитием ранее представленного примера сообщества Engee «Моделирование электрокардиосигнала».
В режиме вариабельной ЧСС в начале очередного кардиоцикла задаётся новое значение частоты сердечных сокращений, после чего рассчитывается соответствующий RR-интервал.
Подключение генератора ЭКС
Функция формирования модельного электрокардиосигнала вынесена в отдельный файл ecg_generator.jl. Файл должен находиться в одной директории с интерактивным скриптом.
include(joinpath(@__DIR__, "ecg_generator.jl"))
Параметры сигнала
fs = 1000.0
duration = 180.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
)
println("Частота дискретизации: ", fs, " Гц")
println("Длительность записи: ", duration, " с")
println("Средняя заданная ЧСС: ", hr, " уд/мин")
Отобразим первые 10 с сигнала. На временной реализации видна последовательность кардиоциклов с небольшой вариабельностью RR-интервалов.
idx = t .<= 10.0
plot(t[idx], ecg[idx],
xlabel="Время, с",
ylabel="Амплитуда, мВ",
label="ЭКС",
title="Модельный электрокардиосигнал",
legend_position=:outertopright)
Амплитудный спектр ЭКС
Для оценки частотного состава выполняется быстрое преобразование Фурье и формируется односторонний амплитудный спектр. Расчёт выполняется по всей сформированной записи.
Функция расчёта амплитудного спектра
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_ecg, A_ecg = spect_amp(ecg, fs)
idx_f = f_ecg .<= 100.0
plot(f_ecg[idx_f], A_ecg[idx_f],
xlabel="Частота, Гц",
ylabel="Амплитуда, мВ",
label="ЭКС",
title="Амплитудный спектр модельного ЭКС",
legend_position=:outertopright)
Моделирование сетевой помехи
Сетевую наводку представим гармоническим сигналом
где – амплитуда помехи, – частота сетевой составляющей, – начальная фаза. В примере используется частота 50 Гц.
f0 = 50.0
A_net = 0.50
phi_net = 0.35
net = A_net .* sin.(2π .* f0 .* t .+ phi_net)
ecg_noisy = ecg .+ net
idx = t .<= 10.0
p = plot(t[idx], ecg[idx],
xlabel="Время, с",
ylabel="Амплитуда, мВ",
label="Исходный ЭКС",
title="Моделирование сетевой помехи 50 Гц",
legend_position=:outertopright)
plot!(p, t[idx], ecg_noisy[idx], label="ЭКС с сетевой помехой")
Адаптивное подавление сетевой помехи
Пусть – входной ЭКС, содержащий полезную составляющую и сетевую помеху, а – текущая оценка помехи. Для запуска рекуррентного алгоритма первые две оценки принимаются равными входным отсчётам:
Для гармонической составляющей с номинальной частотой следующий отсчёт можно определить по двум предыдущим значениям. Предварительная оценка помехи рассчитывается как
Здесь – частота дискретизации.
После прогнозирования определяется невязка между соседними отсчётами остаточного сигнала:
По знаку невязки прогнозируемое значение корректируется на фиксированный инкремент :
После уточнения оценки сетевой составляющей выходной сигнал определяется выражением
Таким образом, на каждом шаге фильтр прогнозирует очередной отсчёт периодической помехи, уточняет его по невязке и вычитает полученную оценку из входного ЭКС.
Инкремент определяет скорость адаптации. При увеличении фильтр быстрее подстраивается к текущей помехе, однако возрастает остаточная ошибка. Одновременно изменяется диапазон отклонения фактической частоты помехи от номинальной частоты , в котором сохраняется эффективная работа алгоритма:
где – амплитуда сетевой помехи.
Величина характеризует допустимое отклонение частоты помехи относительно . При увеличении этот диапазон расширяется. При увеличении амплитуды помехи без изменения , наоборот, диапазон уменьшается. Например, при и Гц получаем Гц.
Функция адаптивной фильтрации
Реализуем адаптивный алгоритм подавления сетевой помехи в виде отдельной функции.
function adaptive_filter(x, fs; f0=50.0, d=0.002)
N = length(x)
e = zeros(Float64, N)
y = zeros(Float64, N)
e[1] = x[1]
e[2] = x[2]
y[1] = x[1] - e[1]
y[2] = x[2] - e[2]
c = cos(2π * f0 / fs)
for n in 2:N-1
e_pred = 2 * c * e[n] - e[n-1]
mismatch = (x[n+1] - e_pred) - (x[n] - e[n])
if mismatch > 0
e[n+1] = e_pred + d
elseif mismatch < 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, e
end
Фильтрация сетевой помехи
Для выбранной амплитуды помехи используется инкремент мВ. Такое значение позволяет на временной реализации наблюдать начальную настройку оценки помехи, после чего фильтр переходит в установившийся режим.
d = 0.002
ecg_filt, net_est = adaptive_filter(ecg_noisy, fs; f0=f0, d=d)
idx = t .<= 10.0
p = plot(t[idx], ecg[idx],
xlabel="Время, с",
ylabel="Амплитуда, мВ",
label="Исходный ЭКС",
title="Результат адаптивной фильтрации",
legend_position=:outertopright)
plot!(p, t[idx], ecg_filt[idx], label="После фильтрации")
В начале записи оценка формируется из первых двух отсчётов входного сигнала и ещё не соответствует сетевой составляющей. Поэтому на выходе наблюдается переходный процесс. По мере последовательной коррекции на величину оценка приближается к помехе 50 Гц, после чего остаточные высокочастотные колебания уменьшаются.
Спектр до и после фильтрации
Для оценки эффективности подавления сетевой помехи сравним амплитудные спектры ЭКС до и после адаптивной фильтрации. Спектральный анализ выполняется по всей длительности сформированной записи.
f_noisy, A_noisy = spect_amp(ecg_noisy, fs)
f_filt, A_filt = spect_amp(ecg_filt, fs)
idx_f = f_noisy .<= 100.0
p = plot(f_noisy[idx_f], A_noisy[idx_f],
xlabel="Частота, Гц",
ylabel="Амплитуда, мВ",
label="До фильтрации",
title="Амплитудный спектр до и после фильтрации",
legend_position=:outertopright)
plot!(p, f_filt[idx_f], A_filt[idx_f], label="После фильтрации")
После адаптивной обработки амплитуда узкополосной составляющей 50 Гц уменьшается. При этом спектральные компоненты модельного ЭКС вне области сетевой помехи сохраняются.
Использованные источники
- Фролов А. В., Дик С. К., Лихачевский Д. В., Давыдов М. В., Чечик Н. М. Цифровая обработка биомедицинских сигналов и изображений. Минск: БГУИР, 2016.
Заключение
В примере рассмотрено адаптивное подавление гармонической сетевой помехи частотой 50 Гц в модельном электрокардиосигнале. Алгоритм прогнозирует следующий отсчёт периодической помехи по двум предыдущим оценкам и уточняет его на фиксированную величину в соответствии со знаком невязки.
В начале обработки наблюдается переходный процесс, связанный с настройкой оценки помехи по входному сигналу. После установления оценки сетевой составляющей высокочастотные колебания на временной реализации уменьшаются, что также подтверждается снижением амплитуды компоненты 50 Гц в спектре обработанного ЭКС.