data = load("16_5_2025_VRMedical.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);


label = {'PP3', 'PP2', 'PP1', 'ZP3', 'ZP2', 'ZP1A', 'ZP1B', 'PL2'}
length(emgBipol(:,2))

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

figure;
plot(emgBipol(:,4))
figure;
plot(emgBipol(:,1))
 start_of_walk = 93340
end_of_walk = 330000

walk = emgBipol(start_of_walk:end_of_walk, :);

for i = 1:8
    figure;
    plot(walk(:,i));
    ylabel(label{i})
    xlabel('vzorky')
end

%


Fs = 1000;
sig = walk(:,5);
L = length(sig);
Y = fft(sig);
P = abs(Y/L);
f = linspace(0, Fs/2, floor(L/2)+1);

plot(f, P(1:floor(L/2)+1));
xlabel('Frekvence [Hz]');
ylabel('Amplituda');
title('Spektrum signálu ZP2');
grid on;

%


Fs = 1000;         % vzorkovací frekvence
f0 = 160;          % frekvence vibrací (k potlačení)
bw = 10;            % šířka pásma zádrže (např. ±5 Hz)

Fs = 1000; % vzorkovací frekvence

low = 70 / (Fs/2);
high = 130 / (Fs/2);

% Bandpass filtr 4. řádu
[b, a] = butter(4, [low, high], 'bandpass');

filt_sig = filtfilt(b, a, sig);

plot(filt_sig)
ylabel(label(5))
xlabel('vzorky')
Y1 = fft(sig);
P1 = abs(Y1 / length(sig));

Y2 = fft(filt_sig);
P2 = abs(Y2 / length(filt_sig));

f = Fs * (0:floor(length(sig)/2)) / length(sig);

figure;
plot(f, P2(1:length(f)));
xlabel('Frekvence [Hz]'); ylabel('Amplituda');

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
low = 70 / (Fs/2);
high = 130 / (Fs/2);

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

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

    [b, a] = butter(4, [low, high], 'bandpass');
    
    % Aplikace
    filt_sig = filtfilt(b, a, sig);

    % Obálka + vyhlazení
    env = abs(hilbert(filt_sig));

    [b, a] = butter(4, 3/(fs/2), 'low'); 

    env_smooth = filtfilt(b, a, env);

    % Uložení obálky
    envelopes_all{i} = env_smooth;

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

%


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

    % Vykreslení průměrné obálky
    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é segmenty', '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;

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

    % Průměrná obálka
    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 – 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;