Адаптивная фильтрация сетевой помехи электрокардиосигнала

using Random

"""
    generate_ecg(; fs=1000.0, duration=180.0, hr=75.0, amplitude=1.0,
                   variable_hr=true, hr_sigma=3.0, phase0=0.0, seed=1234)

Формирует модельный электрокардиосигнал как сумму несимметричных
гауссовых компонент P, Q, R1, R2, S, ST и T.

Параметры:
- `fs` – частота дискретизации, Гц;
- `duration` – длительность сигнала, с;
- `hr` – средняя частота сердечных сокращений, уд/мин;
- `amplitude` – коэффициент масштаба амплитуды;
- `variable_hr` – включение вариабельности ЧСС;
- `hr_sigma` – стандартное отклонение ЧСС, уд/мин;
- `phase0` – начальное положение внутри кардиоцикла, с;
- `seed` – начальное значение генератора случайных чисел.

Возвращает:
- `t` – временная ось, с;
- `ecg` – модельный ЭКС, мВ.
"""
function generate_ecg(; fs=1000.0, duration=180.0, hr=75.0, amplitude=1.0,
    variable_hr=true, hr_sigma=3.0, phase0=0.0, seed=1234)

    fs > 0 || error("Частота дискретизации должна быть больше нуля.")
    duration > 0 || error("Длительность сигнала должна быть больше нуля.")
    40 <= hr <= 180 || error("ЧСС должна находиться в диапазоне 40–180 уд/мин.")
    hr_sigma >= 0 || error("Стандартное отклонение ЧСС не может быть отрицательным.")

    N = Int(round(duration * fs))
    Ts = 1 / fs
    t = collect(0:N-1) .* Ts

    # Параметры компонент P, Q, R1, R2, S, ST и T
    A = [0.11, -0.10, 0.95, 0.02, -0.18, 0.00, 0.20]
    mu = [0.08, 0.14, 0.16, 0.17, 0.19, 0.22, 0.27]
    sigma_l = [0.020, 0.008, 0.006, 0.005, 0.008, 0.020, 0.025]
    sigma_r = [0.025, 0.008, 0.008, 0.006, 0.010, 0.020, 0.035]
    smooth = 0.001

    ecg = zeros(Float64, N)
    hr_current = Float64(hr)
    phase = clamp(Float64(phase0), 0.0, 60 / hr_current)
    rng = MersenneTwister(seed)

    for n in 1:N
        for k in eachindex(A)
            d = phase - mu[k]
            transition = 1 / (1 + exp(-clamp(d / smooth, -60.0, 60.0)))
            sigma = sigma_l[k] + (sigma_r[k] - sigma_l[k]) * transition
            ecg[n] += A[k] * exp(-d^2 / (2 * sigma^2))
        end

        ecg[n] *= amplitude

        rr = 60 / hr_current
        phase_next = phase + Ts

        if phase_next >= rr
            phase = 0.0
            variable_hr && (hr_current = clamp(hr + hr_sigma * randn(rng), 40.0, 180.0))
        else
            phase = phase_next
        end
    end

    return t, ecg
end