The effect of the residual function on adaptive ECG filtering
The discrepancy function determines how the adaptive filter adjusts the network interference estimate. The example considers the discrepancy between the current remainder and the discrepancy between the difference of neighboring residues. The model electrocardiosignal (ECS) shows how these methods affect the restoration of the waveform and the sensitivity of the filter to the displacement of the isoline.
Display Settings
Function engee.clear() cleans up the workspace:
engee.clear()
Two graph rendering modes are available to visualize the results.:
- for static fast rendering, it is used
gr(). - for interactive visualization with the ability to scale and view values, it is used
plotlyjs().
Select one of the modes by commenting out the second one. This example uses interactive mode.
# gr()
plotlyjs()
Installing the necessary libraries
A library is used to generate heart rate variability. Random. It is checked for its presence before connecting. if the library is missing, it is executed Pkg.add(). Function fft It is used from the EngeeDSP library to calculate the amplitude spectrum of the ECS.
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
Formation of a model electrocardiosignal
The source file is being used ecg_generator.jl, located next to the script. The function forms the ECS as the sum of the asymmetric Gaussian components of the cardiocycle. The heart rate varies slightly from cycle to cycle. Fixed value seed provides reproducibility.
The signal is generated once over the entire 60 s interval. Parameters view_start and view_duration only the displayed area is determined and the calculation duration is not changed.
include(joinpath(@__DIR__, "ecg_generator.jl"))
fs = 1000.0 # Sampling rate, Hz
duration = 60.0 # Duration of the simulation, s
hr = 75.0 # Average heart rate, beats/min
hr_sigma = 3.0 # Standard deviation of heart rate, beats/min
phase0 = 0.16 # The initial position inside the heart cycle, with
t, ecg = generate_ecg(
fs=fs, # Sampling rate, Hz
duration=duration, # Recording time, s
hr=hr, # Average heart rate, beats/min
amplitude=1.0, # Amplitude scale factor
variable_hr=true, # Enabling heart rate variability
hr_sigma=hr_sigma, # Standard deviation of heart rate, beats/min
phase0=phase0, # The initial position inside the heart cycle, with
seed=1234 # Generator grain for reproducibility
)
view_start = 0.0 # The beginning of the displayed section, from
view_duration = 5.0 # The duration of the displayed section, with
idx = (t .>= view_start) .& (t .< view_start + view_duration)
plot(t[idx], ecg[idx],
xlabel="Time, from",
ylabel="Amplitude, mV",
label="EX",
title="Model electrocardiosignal",
ylims=(-0.5, 1.25),
legend_position=:outertopright
)
The oscillogram shows repeating heart cycles. Recording begins near the R-wave, so the first samples of the useful signal noticeably affect the initial transient of the filters.
Simulation of network interference
To the useful EX, we add harmonic interference.:
Here mV and Hz. The resulting input array will be common to both filters.
f0 = 50.0 # Rated frequency of network interference, Hz
A_net = 0.5 # Amplitude of network interference, mV
d = 0.002 # Adaptive correction step, mV
net = A_net .* sin.(2π .* f0 .* t)
ecg_noisy = ecg .+ net
p = plot(t[idx], ecg[idx],
xlabel="Time, from",
ylabel="Amplitude, mV",
label="The original EX",
title="Problems with network interference",
ylims=(-1.0, 1.5),
legend_position=:outertopright
)
plot!(p, t[idx], ecg_noisy[idx], label="Problems with interference")
Network fluctuations are superimposed on the teeth of the ECS and are especially noticeable between the cardiocycles. The adaptive filter should evaluate this harmonic and subtract it from the input.
The principle of adaptive filtering
The harmonic signal can be predicted from the previous two counts. Therefore, a preliminary estimate of the next interference value is calculated as
After the forecast, the discrepancy is determined. Let's compare the two options:
- Current balance: .
- Balance difference: .
In the first case, the adaptation reacts to the absolute value of the remainder, including the useful EX. In the second case, the constant component of the remainder is reduced, and slow changes are weakened. However, rapid changes in the QRS complex affect both options. In addition, a constant offset changes the initial conditions, so the difference option is not completely independent of it.
In both cases, the forecast is adjusted based on the discrepancy sign.:
Step mV determines the speed of adjustment. A large step accelerates the correction, but increases fluctuations in the assessment and the reaction to a useful signal. The initial grades are set as in the original example: , . Therefore, the first two output counts are zero.
The recurrent forecast and the difference discrepancy correspond to the formulas (4.3) and (4.4) manuals by Frolov et al., p. 49. Estimation correction and interference subtraction are described there by formulas (4.5) and (4.6). The code retains one correction step per count, as in the original demo example Engee.
A filter with a discrepancy based on the current balance
Function adaptive_filter_direct() uses the current balance:
It returns the cleared EX, interference score, and discrepancy. For the first two samples, the discrepancy is not defined and is denoted by NaN.
# f0 is the tuning frequency, Hz; d is the correction step, mV.
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
A filter with a residual difference
Function adaptive_filter_difference() compares the remains of neighboring samples:
The forecast, correction step, initial conditions, and return values are the same as the first option.
# f0 is the tuning frequency, Hz; d is the correction step, mV.
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
Filtering result
Let's apply both options to the same input array. The tuning frequency, pitch, and initial conditions are the same. The graph shows the first 5 seconds, including the initial transition.
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="Time, from",
ylabel="Amplitude, mV",
label="The original EX",
title="Comparison of recovered EX",
legend_position=:outertopright
)
plot!(p, t[idx], y_direct[idx], label="Current balance")
plot!(p, t[idx], y_difference[idx], label="Balance difference")
At the beginning of the recording, the interference estimate contains a useful component of the first samples. As you adapt, this error decreases. To view another section, change view_start and re-execute the cells, starting from the formation of the EX. The simulation duration will remain 60 seconds.
p = plot(t[idx], net[idx],
xlabel="Time, from",
ylabel="Amplitude, mV",
label="The true hindrance",
title="Network interference assessment",
legend_position=:outertopright
)
plot!(p, t[idx], e_direct[idx], label="Current balance")
plot!(p, t[idx], e_difference[idx], label="Balance difference")
A comparison with the known interference shows the initial transition process and the deviations of the estimates. To view individual periods in detail, you can enlarge the desired area using the interactive graph.
Comparison of discrepancy functions
Let's construct the time dependencies of the residuals of the two filters. They allow you to see how the current balance and the difference in balances react to ECS sections with slow and fast changes. This comparison explains how the correction works, and the accuracy of the recovery is evaluated separately based on the error of the EX.
p = plot(t[idx], mismatch_direct[idx],
xlabel="Time, from",
ylabel="Discrepancy, mV",
label="Current balance",
title="Two ways to form a discrepancy",
legend_position=:outertopright
)
plot!(p, t[idx], mismatch_difference[idx], label="Balance difference")
The difference discrepancy has a smaller contribution of the slow components. Its amplitude cannot be compared with the direct discrepancy as an indicator of quality: these values are defined in different ways. The quality is assessed by the error of the EX recovery.
error_direct = y_direct .- ecg
error_difference = y_difference .- ecg
p = plot(t[idx], error_direct[idx],
xlabel="Time, from",
ylabel="Error, mV",
label="Current balance",
title="EX Recovery Error",
legend_position=:outertopright
)
plot!(p, t[idx], error_difference[idx], label="Balance difference")
# The maximum modulo error after the initial transient.
idx_after_transient = t .>= 2.0 # Error estimation interval, 2-60 seconds
println("Without displacement of the isoline. Maximum absolute error in the range of 2-60 s, mV")
println("Current balance: ", round(maximum(abs.(error_direct[idx_after_transient])), digits=5))
println("Balance difference: ", round(maximum(abs.(error_difference[idx_after_transient])), digits=5))
Without shifting the isoline, both filters suppress network fluctuations. The maximum recovery errors are estimated in the range of 2-60 seconds, without the initial phase of the transition process. The error spikes near the QRS complex reflect the adaptation response to rapid EC changes.
Amplitude spectra
The function is used spect_amp() from the original example. A one-sided amplitude spectrum is calculated, normalized by the number of samples. The spectra are plotted for all 60 s. The display limit to 5 s applies only to time graphs.
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 # The upper limit of the displayed spectrum, Hz
p = plot(f[idx_f], A_clean[idx_f],
xlabel="Frequency, Hz",
ylabel="Amplitude, mV",
label="The original EX",
title="Spectra before and after filtration",
legend_position=:outertopright
)
plot!(p, f[idx_f], A_noisy[idx_f], label="Problems with interference")
plot!(p, f[idx_f], A_direct[idx_f], label="Current balance")
plot!(p, f[idx_f], A_difference[idx_f], label="Balance difference")
Reducing the 50 Hz line demonstrates suppression of interference. However, the useful EX also contains components in this region, so the spectral graph should be considered together with the error relative to the known pure signal.
Filtration of the EX with an offset of the isoline
Let's add a constant offset to the same model value:
The input of both filters is equal to . Parameter offset sets a constant offset in mV. The network interference remains the same. Filters should suppress it and preserve the offset of the useful signal.
offset = 0.3 # Constant displacement of the isoline EX, mV
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="Time, from",
ylabel="Amplitude, mV",
label="The reference EX",
title="The displacement of the isoline is +0.3 mV",
legend_position=:outertopright
)
plot!(p, t[idx], offset_direct.y[idx], label="Current balance")
plot!(p, t[idx], offset_difference.y[idx], label="Balance difference")
When the EX is shifted, noticeable fluctuations remain in the filter with the current remainder. The difference version preserves the waveform better relative to the offset reference. This is because the constant component is reduced when calculating the remainder difference.
To check the result, the error is shown below regarding ecg_offset. The added offset is included in the reference signal and is not considered interference.
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="Time, from",
ylabel="Error, mV",
label="Current balance",
title="Error restoring an EX with an offset",
legend_position=:outertopright
)
plot!(p, t[idx], error_offset_difference[idx], label="Balance difference")
println("With an offset of the isoline +0.3 mV. Maximum absolute error in the range of 2-60 s, mV")
println("Current balance: ", round(maximum(abs.(error_offset_direct[idx_after_transient])), digits=5))
println("Balance difference: ", round(maximum(abs.(error_offset_difference[idx_after_transient])), digits=5))
When the isoline is shifted, the direct version retains pronounced error fluctuations. The difference variant reacts less strongly to the constant component of the remainder. The maximum errors are given for the same interval of 2-60 s.
Influence of the choice of the residual function
The graphs show that the way the discrepancy is formed changes the nature of the residual error. The current remainder takes into account the absolute value of the useful signal and can guide the correction even with an accurate prediction of interference. The difference in residues reduces the constant component, but rapid changes in the QRS complex continue to affect adaptation.
When the isoline is shifted, the difference between the variants becomes more noticeable. The maximum absolute error in the range of 2-60 s characterizes the largest deviation of the restored EC from the reference after the exclusion of the initial phase of the transient process. This indicator does not describe the average error and all the features of preserving the waveform.
The results obtained show the influence of two criteria with the specified parameters of the EX and filters. To select a criterion in another task, it is necessary to take into account the shape of the useful signal, the nature of the interference, and the requirements for acceptable distortion.
Used sources
- Frolov A.V., Dik S. K., Likhachevsky D. V., Davydov M. V., Chechik N. M. Digital processing of biomedical signals and images. Minsk: BGUIR, 2016.
- Adaptive filtering of network interference of an electrocardiosignal - community Engee.