Биорадар: Часть 1. Анализ сигнала радарного и ЭКГ сигналов

using MAT
using PyCall
skimage = pyimport("skimage");
function load_data(filepath::String; offset::Int=0)
    data = matread(filepath)
    
    radar_i = vec(data["radar_i"])[offset+1:end]
    radar_q = vec(data["radar_q"])[offset+1:end]
    fs_radar = data["fs_radar"]
    ecg_data1 = data["tfm_ecg1"]
    ecg_data2 = data["tfm_ecg2"]
    fs_ecg = data["fs_ecg"]
    
    return radar_i, radar_q, fs_radar, ecg_data1, ecg_data2, fs_ecg
end

function EllipseReconstruction(i, q)
    model = skimage.measure.EllipseModel();
    stackData = hcat(i, q);
    model.estimate(stackData)
    xc, yc, a, b, theta = model.params;
    IQ_meas = Complex.(i, q);
    I = real.(IQ_meas) .- xc;
    Q = imag.(IQ_meas) .- yc;
    phi  = -theta;

    Arot = [[ cos(phi),  sin(phi)],
    [-sin(phi),  cos(phi)]]

    IQ1 = hcat(I, Q) * Arot;
    IQ1 = mapreduce(permutedims, vcat, IQ1)

    I2 = IQ1[:, 1] ./ a;
    Q2 = IQ1[:, 2] ./ b;
    circle = Complex.(I2, Q2);
    I_comp = real.(circle);
    Q_comp = imag.(circle);


    sigma = angle.(circle);
    lambda_ = 12.5e-3;
    range_ = unwrap(sigma) ./ (4 * pi) * lambda_;
    return I_comp, Q_comp, sigma, range_
end
function Filter(data, lowcut, highcut, fs, order)
    responsetype = Bandpass(lowcut, highcut)
    designmethod = Butterworth(order)
    filter = digitalfilter(responsetype, designmethod; fs=fs)
    return filtfilt(filter, data)
end

function getVitalSignals(radar_dist, fs_radar, order)
    radar_respiration = Filter(radar_dist, 0.05, 0.5, fs_radar, order)
    radar_pulse = Filter(radar_dist, 1.0, 6.0, fs_radar, order)
    radar_heartsound = Filter(radar_dist, 10.0, 80.0, fs_radar, order)

    radar_heartsound_denoise = skimage.restoration.denoise_wavelet(
        radar_heartsound,
        wavelet="db10",
        mode="soft",
        wavelet_levels=5,
        method="VisuShrink",
        rescale_sigma=true
    )

    return radar_respiration, radar_pulse, radar_heartsound, radar_heartsound_denoise
end


function plot_ecg_with_peaks(ecg_cleaned, waves, rpeaks; range_start=90000, range_end=100000)
    # 1. Извлекаем подмножества
    r_subset = filter(x -> range_start ≤ x ≤ range_end, rpeaks["ECG_R_Peaks"])
    r_local = r_subset .- (range_start - 1)

    p_subset = filter(x -> range_start ≤ x ≤ range_end, waves["ECG_P_Peaks"])
    p_local = p_subset .- (range_start - 1)

    q_subset = filter(x -> range_start ≤ x ≤ range_end, waves["ECG_Q_Peaks"])
    q_local = q_subset .- (range_start - 1)

    s_subset = filter(x -> range_start ≤ x ≤ range_end, waves["ECG_S_Peaks"])
    s_local = s_subset .- (range_start - 1)

    t_subset = filter(x -> range_start ≤ x ≤ range_end, waves["ECG_T_Peaks"])
    t_local = t_subset .- (range_start - 1)

    ecg_part = ecg_cleaned[range_start:range_end]

    # 2. Рисуем график
    p = plot(
        layout = (1, 2),
        size   = (900, 400),
        theme  = :default,
        framestyle = :box,
        grid  = true
    )

    # Левый: ECG + пики
    plot!(p, ecg_part;
        subplot = 1,
        color   = :black,
        lw      = 1.5,
        label   = "ECG",
        xlabel  = "Отсчёт",
        ylabel  = "Амплитуда",
        title   = "Фрагмент с пиками"
    )
    scatter!(p, r_local, ecg_cleaned[r_subset]; subplot = 1, marker = :circle,   mc = :red,    ms = 5, label = "R-пики")
    scatter!(p, p_local, ecg_cleaned[p_subset]; subplot = 1, marker = :pentagon, mc = :green,  ms = 6, label = "P-пики")
    scatter!(p, q_local, ecg_cleaned[q_subset]; subplot = 1, marker = :diamond,  mc = :blue,   ms = 5, label = "Q-пики")
    scatter!(p, s_local, ecg_cleaned[s_subset]; subplot = 1, marker = :utriangle,mc = :purple, ms = 5, label = "S-пики")
    scatter!(p, t_local, ecg_cleaned[t_subset]; subplot = 1, marker = :dtriangle,mc = :orange, ms = 5, label = "T-пики")

    # Правый: просто ECG
    plot!(p, ecg_part;
        subplot = 2,
        color   = :black,
        lw      = 1.2,
        label   = "",
        xlabel  = "Отсчёт",
        title   = "Чистый участок ЭКГ"
    )

    return p
end