data = load("30_4_2025_EMG+NOpiezo.mat")
emg = data.data;
emgBipol = zeros(size(emg,1), 8);
emgBipol(:,1) = emg(:,2) - emg(:,1);
emgBipol(:,2) = emg(:,4) - emg(:,3);
emgBipol(:,3) = emg(:,6) - emg(:,5);
emgBipol(:,4) = emg(:,8) - emg(:,7);
emgBipol(:,5) = emg(:,10) - emg(:,9);
emgBipol(:,6) = emg(:,13) - emg(:,12);
emgBipol(:,7) = emg(:,13) - emg(:,11);
emgBipol(:,8) = emg(:,15) - emg(:,14);

label = {'PP3', 'PP2', 'PP1', 'ZP3', 'ZP2', 'ZP1A', 'ZP1B', 'PL2'}


figure;
for i = 1:8
    subplot(8,1,i);
    plot(emgBipol(:,i));
    ylabel(label{i})
end
xlabel('Časové vzorky');
grid on;
figure;
plot(emgBipol(:,2))


%


 start_of_walk = 44537
end_of_walk = 127071

walk = emgBipol(start_of_walk:end_of_walk, :);

for i = 1:8
    figure;
    plot(walk(:,i));
    ylabel(label{i})
    xlabel('vzorky')
end
xlabel('Časové vzorky');
grid on;
walk = walk(15800:23300,:);
size(walk)
for i = 1:8
    figure;
    plot(walk(:,i));
    ylabel(label{i})
end

%


fs = 1000;               % vzorkovací frekvence
freq_band = [30 250];    % frekvenční rozsah pro wavelet filtraci
labels = {'P5', 'P3A', 'P3B', 'P1A', 'P1B', 'L3', 'P6', 'P4', 'P2'};  % názvy kanálů

[b, a] = butter(4, 3/(fs/2), 'low');  % dolní propust pro vyhlazení obálky

num_channels = size(walk, 2);
envelopes_all = cell(1, num_channels);  % výstupní obálky

for i = 1:num_channels
    sig = walk(:,i); 

    % Wavelet filtr
    fb = cwtfilterbank(SignalLength=length(sig), SamplingFrequency=fs, ...
        FrequencyLimits=freq_band);

    sig_cwt = wt(fb, sig);
    sig_filtered = icwt(sig_cwt, 'SignalMean', mean(sig));
    env = abs(hilbert(sig_filtered));
    env_smooth = filtfilt(b, a, env);
    envelopes_all{i} = env_smooth;


    figure;
    subplot(3,1,1); plot(sig); title(['EMG signál – ', label{i}]); grid on;
    subplot(3,1,2); plot(sig_filtered); title('Wavelet filtrovaný'); grid on;
    subplot(3,1,3); plot(env_smooth, 'r', 'LineWidth', 1.5);
    title('Obálka (wavelet + Hilbert)'); xlabel('Časové vzorky'); grid on;
end

%



fs = 1000;  
num_channels = length(envelopes_all);
period_lengths = zeros(1, num_channels);

for i = 1:num_channels
    sig = envelopes_all{i};
    sig = sig(:);
    [acf, lags] = xcorr(sig - mean(sig), 'coeff');
    center = ceil(length(acf)/2);
    lags_right = lags(center+1:end);
    acf_right = acf(center+1:end);

    min_shift = 900;
    search_range = find(lags_right > min_shift);
    [pks, locs] = findpeaks(acf_right(search_range));

    if isempty(pks)
        warning("Chybí vrchol v autokorelaci pro kanál %s", label{i});
        period_lengths(i) = NaN;
    else
        period = lags_right(search_range(locs(1)));
        fprintf("Kanál %s – perioda: %.3f s\n", label{i}, period / fs);

        period_lengths(i) = period;
    end

    figure;
    plot(lags_right, acf_right);
    hold on;
    xline(period, '--r', 'Perioda');
    title(['Autokorelační funkce – ' label{i}]);
    xlabel('Zpoždění (vzorky)');
    ylabel('Korelace');
    grid on;

end

valid_periods = period_lengths(~isnan(period_lengths) & period_lengths > 0);
mean_period = round(mean(valid_periods));
fprintf('Použitá průměrná perioda: %d vzorků (%.2f s)\n', mean_period, mean_period/fs);

mean_single_steps = cell(1, num_channels);

for i = 1:num_channels
    sig = envelopes_all{i};
    sig = sig(:);

    num_segments = floor(length(sig) / mean_period);
    sig_cut = sig(1 : num_segments * mean_period);
    segments = reshape(sig_cut, mean_period, num_segments)';

    mean_segment = mean(segments, 1);
    mean_single_steps{i} = mean_segment;

 
    time_axis = (0:length(mean_segment)-1) / fs;


    figure;
    hold on;
    h_segments = gobjects(size(segments, 1), 1);
    for j = 1:size(segments, 1)
        h_segments(j) = plot(time_axis, segments(j, :), 'Color', [0.6 0.6 1 0.4]);  % průhledná modrá
    end

    if all(isnan(mean_segment))
        warning("Průměrná obálka pro %s obsahuje jen NaN!", label{i});
        h_mean = NaN;
    elseif all(mean_segment == 0)
        warning("Průměrná obálka pro %s je nulová!", label{i});
        h_mean = NaN;
    else
        h_mean = plot(time_axis, mean_segment, 'Color', [1 0 0], 'LineWidth', 2);
        uistack(h_mean, 'top');
    end

    title(['Průměrný krok – ' label{i}]);
    xlabel('Čas [s]');
    ylabel('Amplituda obálky');

    if isgraphics(h_mean)
        legend([h_segments(1), h_mean], {'Jednotlivé kroky', 'Průměr'});
    else
        legend(h_segments(1), 'Jednotlivé kroky');
    end

    grid on;
end
figure;
numRows = ceil(sqrt(num_channels));
numCols = ceil(num_channels / numRows);

for i = 1:num_channels
    sig = envelopes_all{i};
    sig = sig(:);

    num_segments = floor(length(sig) / mean_period);
    sig_cut = sig(1 : num_segments * mean_period);
    segments = reshape(sig_cut, mean_period, num_segments)';

    mean_segment = mean(segments, 1);
    mean_single_steps{i} = mean_segment;

    time_axis = (0:length(mean_segment)-1) / fs;

    subplot(numRows, numCols, i);
    hold on;


    for j = 1:size(segments, 1)
        plot(time_axis, segments(j, :), 'Color', [0.6 0.6 1 0.3]);
    end

    if ~all(isnan(mean_segment)) && ~all(mean_segment == 0)
        plot(time_axis, mean_segment, 'Color', [1 0 0], 'LineWidth', 2);
    end

    title(label{i});
    xlabel('Čas [s]');
    ylabel('Amplituda');
    grid on;
end

sgtitle('Průměrné obálky a segmenty – všechny kanály');

figure;
hold on;
time_axis = (0:mean_period-1) / fs;
for i = 1:num_channels
    if isempty(mean_single_steps{i})
        continue;
    end
    plot(time_axis, mean_single_steps{i}, 'DisplayName', label{i});
end
title('Zarovnané průměrné obálky (společná perioda)');
xlabel('Čas [s]');
ylabel('Amplituda obálky');
legend('show');
grid on;

figure;
hold on;
time_axis = (0:mean_period-1) / fs;
for i = 1:num_channels
    if isempty(mean_single_steps{i})
        continue;
    end


    obalka = mean_single_steps{i};
    obalka_norm = obalka / max(obalka);

    plot(time_axis, obalka_norm, 'DisplayName', label{i});
end
title('Zarovnané průměrné obálky (normalizované, společná perioda)');
xlabel('Čas [s]');
ylabel('Normalizovaná amplituda obálky');
legend('show');
grid on;



%


disp(period_lengths);
disp(period_lengths / fs); 

%


upper_threshold_factor = 0.55;
lower_threshold_factor = 0.8;

for i = 1:numel(mean_single_steps)
    if isempty(mean_single_steps{i})
        continue;
    end

    sig = mean_single_steps{i};
    sig = sig(:);
    time_axis = (0:length(sig)-1) / fs;

    threshold_up = upper_threshold_factor * max(sig);
    threshold_down = lower_threshold_factor * threshold_up;

    y = zeros(1, length(sig));
    y_old = 0;
    for j = 2:length(sig)
        if y_old == 0
            y(j) = sig(j) > threshold_up;
        else
            y(j) = sig(j) > threshold_down;
        end
        y_old = y(j);
    end
    if y(end) == 1
        y(end+1) = 0;
    end
    active_segments = find(diff([0 y 0]));

    figure;
    hold on;
    plot(time_axis, sig, 'r', 'LineWidth', 2, 'DisplayName', 'Průměrná obálka');
    yline(threshold_up, '--k', 'DisplayName', 'Horní práh');
    yline(threshold_down, ':k', 'DisplayName', 'Dolní práh');


    for k = 1:2:length(active_segments)-1
        if active_segments(k+1) > length(time_axis)
            continue;
        end
        x_start = time_axis(active_segments(k));
        x_end = time_axis(active_segments(k+1));
        segment_indices = active_segments(k):active_segments(k+1);
        segment_signal = sig(segment_indices);
        peak_height = max(segment_signal);

        plot(time_axis(segment_indices(segment_signal == peak_height)), ...
             peak_height, 'ko', 'MarkerFaceColor', 'k');

        rectangle('Position', [x_start, 0, x_end - x_start, peak_height], ...
                  'EdgeColor', [0.2, 0.8, 0.2], 'LineWidth', 2, 'LineStyle', '--');
    end

    title(['Průměrná obálka – ' label{i}]);
    xlabel('Čas [s]');
    ylabel('Amplituda');
    legend('show', 'Location', 'northeastoutside');
    grid on;
end

idx_PP1  = find(strcmp(label, 'PP1'));
idx_PP2  = find(strcmp(label, 'PP2'));
idx_PP3  = find(strcmp(label, 'PP3'));
idx_ZP1A = find(strcmp(label, 'ZP1A'));
idx_ZP1B = find(strcmp(label, 'ZP1B'));
idx_ZP2  = find(strcmp(label, 'ZP2'));
idx_ZP3  = find(strcmp(label, 'ZP3'));

PP1  = mean_single_steps{idx_PP1};
PP2  = mean_single_steps{idx_PP2};
PP3  = mean_single_steps{idx_PP3};
ZP1A = mean_single_steps{idx_ZP1A};
ZP1B = mean_single_steps{idx_ZP1B};
ZP2  = mean_single_steps{idx_ZP2};
ZP3  = mean_single_steps{idx_ZP3};

period_samples = period_lengths(idx_PP1);
half_period_time = period_samples / 2 / fs;

time_axis = (0:length(PP1)-1) / fs;

L_PP1  = PP1;
L_PP2  = PP2;
L_PP3  = PP3;
L_ZP1A = ZP1A;
L_ZP1B = ZP1B;
L_ZP2  = ZP2;
L_ZP3  = ZP3;

figure;
hold on;

% Pravá noha
plot(time_axis, PP1,  'b', 'DisplayName', 'PP1');
plot(time_axis, PP2,  'g', 'DisplayName', 'PP2');
plot(time_axis, PP3,  'm', 'DisplayName', 'PP3');
plot(time_axis, ZP1A, 'c', 'DisplayName', 'ZP1A');
plot(time_axis, ZP1B, 'r', 'DisplayName', 'ZP1B');
plot(time_axis, ZP2,  'k', 'DisplayName', 'ZP2');
plot(time_axis, ZP3,  '-', 'Color', [1 0.5 0], 'DisplayName', 'ZP3');

% Levá noha – posuneme časovou osu
plot(time_axis + half_period_time, L_PP1,  '--b', 'DisplayName', 'PL1');
plot(time_axis + half_period_time, L_PP2,  '--g', 'DisplayName', 'PL2');
plot(time_axis + half_period_time, L_PP3,  '--m', 'DisplayName', 'PL3');
plot(time_axis + half_period_time, L_ZP1A, '--c', 'DisplayName', 'ZL1A');
plot(time_axis + half_period_time, L_ZP1B, '--r', 'DisplayName', 'ZL1B');
plot(time_axis + half_period_time, L_ZP2,  '--k', 'DisplayName', 'ZL2');
plot(time_axis + half_period_time, L_ZP3,  '--', 'Color', [1 0.5 0], 'DisplayName', 'ZL3');

% Levá noha – posun do mínusu
plot(time_axis - half_period_time, L_PP1,  '--b', 'HandleVisibility', 'off');
plot(time_axis - half_period_time, L_PP2,  '--g', 'HandleVisibility', 'off');
plot(time_axis - half_period_time, L_PP3,  '--m', 'HandleVisibility', 'off');
plot(time_axis - half_period_time, L_ZP1A, '--c', 'HandleVisibility', 'off');
plot(time_axis - half_period_time, L_ZP1B, '--r', 'HandleVisibility', 'off');
plot(time_axis - half_period_time, L_ZP2,  '--k', 'HandleVisibility', 'off');
plot(time_axis - half_period_time, L_ZP3,  '--', 'Color', [1 0.5 0], 'HandleVisibility', 'off');


title('Průměrné obálky – pravá a odhad levé nohy');
xlabel('Čas [s]');
ylabel('Amplituda obálky');
legend('show', 'Location', 'northeastoutside');
grid on;
xlim([0,1.2]);

%


upper_threshold_factor = 0.55;
lower_threshold_factor = 0.8;
fs = 1000;
time_ref = (0:length(PP1)-1) / fs;
half_period = time_ref(end) / 2;
t_limit = [0, 1.2];

shift_signal = @(sig, shift_time) deal(time_ref + shift_time, sig);

AVERAGE = {};
label = {};
time_axes = {};
style = {};

sig_list = {PP1, PP2, PP3, ZP1A, ZP1B, ZP2, ZP3};
names = {'PP1', 'PP2', 'PP3', 'ZP1A', 'ZP1B', 'ZP2', 'ZP3'};
for i = 1:length(sig_list)
    [t, s] = shift_signal(sig_list{i}, 0);
    AVERAGE{end+1} = s;
    time_axes{end+1} = t;
    label{end+1} = names{i};
    style{end+1} = '-';
end

odhad_list = {PP1, PP2, PP3, ZP1A, ZP1B, ZP2, ZP3};
names_odhad = {'L_PP1 (odhad)', 'L_PP2 (odhad)', 'L_PP3 (odhad)', 'L_ZP1A (odhad)', 'L_ZP1B (odhad)', 'L_ZP2 (odhad)', 'L_ZP3 (odhad)'};
for i = 1:length(odhad_list)
    [t, s] = shift_signal(odhad_list{i}, +half_period);
    AVERAGE{end+1} = s;
    time_axes{end+1} = t;
    label{end+1} = names_odhad{i};
    style{end+1} = '--';
end

names_zpet = {'L_PP1 (odhad zpětný)', 'L_PP2 (odhad zpětný)', 'L_PP3 (odhad zpětný)', 'L_ZP1A (odhad zpětný)', 'L_ZP1B (odhad zpětný)', 'L_ZP2 (odhad zpětný)', 'L_ZP3 (odhad zpětný)'};
for i = 1:length(odhad_list)
    [t, s] = shift_signal(odhad_list{i}, -half_period);
    AVERAGE{end+1} = s;
    time_axes{end+1} = t;
    label{end+1} = names_zpet{i};
    style{end+1} = '--';
end

figure;
hold on;
colors = lines(length(AVERAGE));
used_labels = {};

for i = 1:length(AVERAGE)
    t = time_axes{i};
    sig = AVERAGE{i};
    mask = (t >= t_limit(1)) & (t <= t_limit(2));
    t_cut = t(mask);
    s_cut = sig(mask);

    raw_label = label{i};

    clean_label = erase(raw_label, {' (odhad)', ' (odhad zpětný)'});

    clean_label = strrep(clean_label, 'L_PP', 'PL');
    clean_label = strrep(clean_label, 'L_ZP', 'ZL');
    
    clean_label = strtrim(clean_label); 

    if ~ismember(clean_label, used_labels)
        plot(t_cut, s_cut, style{i}, 'Color', colors(i,:), 'LineWidth', 1.2, 'DisplayName', clean_label);
        used_labels{end+1} = clean_label;
    else
        plot(t_cut, s_cut, style{i}, 'Color', colors(i,:), 'LineWidth', 1.2, 'HandleVisibility', 'off');
    end

    threshold_up = upper_threshold_factor * max(s_cut);
    threshold_down = lower_threshold_factor * threshold_up;
    y = zeros(1, length(s_cut));
    y_old = 0;
    for j = 2:length(s_cut)
        if y_old == 0
            y(j) = s_cut(j) > threshold_up;
        else
            y(j) = s_cut(j) > threshold_down;
        end
        y_old = y(j);
    end

    y_ext = [y, 0];
    t_cut_ext = [t_cut, t_cut(end) + 1/fs];
    seg_idx = find(diff([0 y_ext]) ~= 0);

    for k = 1:2:length(seg_idx)-1
        if seg_idx(k+1) > length(s_cut)
            continue;
        end
        x_start = t_cut_ext(seg_idx(k));
        x_end = t_cut_ext(seg_idx(k+1));

        segment = s_cut(seg_idx(k):seg_idx(k+1)-1);
        peak_val = max(segment);

        rectangle('Position', [x_start, 0, x_end - x_start, peak_val], ...
                  'EdgeColor', colors(i,:), 'LineStyle', '--', 'LineWidth', 1.5);
    end
end

title('Průměrné obálky všech elektrod (sjednocené odhady)');
xlabel('Čas [s]');
ylabel('Amplituda');
legend('show', 'Location', 'northeastoutside');
xlim(t_limit);
grid on;

title('Průměrné obálky všech elektrod s aktivními segmenty');
xlabel('Čas [s]');
ylabel('Amplituda');
legend('show', 'Location', 'northeastoutside');
xlim(t_limit);
grid on;

%


upper_threshold_factor = 0.55;
lower_threshold_factor = 0.8;
fs = 1000;
figure;

num_signals = length(AVERAGE);
active_times = cell(num_signals, 1);

clean_names = cellfun(@(lbl) erase(lbl, {' (odhad)', ' (odhad zpětný)'}), label, 'UniformOutput', false);
clean_names = strrep(clean_names, 'L_PP', 'PL');
clean_names = strrep(clean_names, 'L_ZP', 'ZL');
clean_names = strtrim(clean_names);
unique_names = unique(clean_names);

for i = 1:num_signals
    t = time_axes{i};
    sig = AVERAGE{i};

    mask = (t >= t_limit(1)) & (t <= t_limit(2));
    t_cut = t(mask);
    s_cut = sig(mask);

    % Detekce segmentů
    threshold_up = upper_threshold_factor * max(s_cut);
    threshold_down = lower_threshold_factor * threshold_up;

    y = zeros(1, length(s_cut));
    y_old = 0;
    for j = 2:length(s_cut)
        if y_old == 0
            y(j) = s_cut(j) > threshold_up;
        else
            y(j) = s_cut(j) > threshold_down;
        end
        y_old = y(j);
    end

    y_ext = [y, 0];
    t_cut_ext = [t_cut, t_cut(end) + 1/fs];
    seg_idx = find(diff([0 y_ext]) ~= 0);

    if mod(length(seg_idx), 2) == 0
        active_times{i} = [t_cut_ext(seg_idx(1:2:end))', t_cut_ext(seg_idx(2:2:end))'];
    else
        warning("Počet hranic segmentů není sudý pro signál %s", label{i});
    end
end

numPlots = length(unique_names);
numRows = ceil(sqrt(numPlots));
numCols = ceil(numPlots / numRows);

for u = 1:numPlots
    subplot(numRows, numCols, u);
    hold on;
    name = unique_names{u};
    title(name, 'Interpreter', 'none');
    xlabel('Čas [s]');
    ylabel('Amplituda');

    max_val = 0;

    for i = 1:num_signals
        if ~strcmp(clean_names{i}, name)
            continue;
        end

        t = time_axes{i};
        sig = AVERAGE{i};
        segs = active_times{i};

        % Styl čáry podle typu
        if contains(label{i}, '(odhad zpětný)')
            ls = '--'; 
        elseif contains(label{i}, '(odhad)')
            ls = '--';
        else
            ls = '-';
        end

 
        color = colors(i, :);

        % Výřez do limitu
        mask = (t >= t_limit(1)) & (t <= t_limit(2));
        t_cut = t(mask);
        s_cut = sig(mask);

        plot(t_cut, s_cut, ls, 'Color', color, 'LineWidth', 1.2);
        max_val = max(max_val, max(s_cut));

        % Výplň aktivních segmentů
        for k = 1:size(segs,1)
            x1 = segs(k,1);
            x2 = segs(k,2);
            rectangle('Position', [x1, 0, x2 - x1, max(s_cut)], ...
                      'EdgeColor', color, 'LineStyle', '--', 'LineWidth', 1.2);
        end
    end

    ylim([0, max_val*1.1]);
    xlim(t_limit);
    grid on;
end

sgtitle('Sjednocené signály s aktivními segmenty');

%


label_cleaned = label;
for i = 1:length(label)
    if contains(label{i}, '(odhad zpětný)')
        label_cleaned{i} = strrep(label{i}, ' (odhad zpětný)', ' (odhadovaný)');
    elseif contains(label{i}, '(odhad)')
        label_cleaned{i} = strrep(label{i}, ' (odhad)', ' (odhadovaný)');
    end
end

clean_names = cellfun(@(lbl) erase(lbl, {' (odhad)', ' (odhad zpětný)', ' (odhadovaný)'}), label, 'UniformOutput', false);
clean_names = strrep(clean_names, 'L_PP', 'PL');
clean_names = strrep(clean_names, 'L_ZP', 'ZL');
clean_names = strtrim(clean_names);

is_left = contains(label, 'odhad');

electrode_col = {};
side_col = {};
start_col = [];
end_col = [];

unique_clean_names = unique(clean_names);

for u = 1:length(unique_clean_names)
    name = unique_clean_names{u};

    idx_group = find(strcmp(clean_names, name));

    for i = idx_group(:)'
        if isempty(active_times{i})
            continue;
        end

        n_segments = size(active_times{i}, 1);
        names = repmat({name}, n_segments, 1);
        side = repmat({''}, n_segments, 1);

        if is_left(i)
            side(:) = {'Levá'};
        else
            side(:) = {'Pravá'};
        end

        electrode_col = [electrode_col; names];
        side_col = [side_col; side];
        start_col = [start_col; active_times{i}(:, 1)];
        end_col = [end_col; active_times{i}(:, 2)];
    end
end

active_segments_table = table(side_col, electrode_col, start_col, end_col, ...
    'VariableNames', {'Strana', 'Elektroda', 'Start_sec', 'End_sec'});

active_segments_table = sortrows(active_segments_table, 'Start_sec');
writetable(active_segments_table, 'aktivni_useky_sjednocene.xlsx');