% Lab02.m - MATLAB skript pro detekci stimulačních časů a motorických evokovaných odpovědí

clc; clear; close all;

%% Načtení dat
file_path = 'ulnaris01.txt';
if exist(file_path, 'file') ~= 2
    error('Soubor %s nebyl nalezen.', file_path);
end

data = load(file_path);
if size(data,2) < 2
    error('Nesprávný formát souboru. Očekávány alespoň dva sloupce.');
end

stim = data(:,1);       % Stimulační pulzy
emg5000 = data(:,2);    % EMG (0,5-5000 Hz)
fs = 50e3;              % Vzorkovací frekvence
t = (0:length(stim)-1)/fs;

%% Detekce počátků stimulů
stimulus_threshold = 30;                     % Prahová hodnota (45)
refractory_period = round(10e-3 * fs);       % 10 ms refrakterní fáze

potential_indices = find(diff(stim) > 0 & stim(1:end-1) < stimulus_threshold & stim(2:end) >= stimulus_threshold);

stimulus_indices = []; last_index = -refractory_period;
for idx = potential_indices'
    if idx > last_index + refractory_period
        stimulus_indices = [stimulus_indices; idx];
        last_index = idx;
    end
end

%% Detekce odpovědí v EMG

% Prahová detekce
threshold = 0.2;
response_indices = [];
for i = 1:length(stimulus_indices)
    search_range = stimulus_indices(i):min(stimulus_indices(i) + round(50e-3 * fs), length(emg5000));
    response_idx = search_range(find(emg5000(search_range) > threshold, 1, 'first'));
    if ~isempty(response_idx)
        response_indices = [response_indices; response_idx];
    end
end

% Detekce odpovědí pomocí výkonové obálky
% window_size = round(5e-3 * fs);                           % 5ms klouzavé okno
% power_envelope = movmean(abs(emg5000).^2, window_size);
% % ekvivalentní zápis pomocí filtfilt
% % b = ones(1, window_size) / window_size; % FIR filtr s rovnoměrnými koeficienty
% % power_envelope = filtfilt(b, 1, abs(emg5000).^2);
% threshold_env = 0.01 * max(power_envelope);               % Prahová hodnota
% response_indices = [];
% for i = 1:length(stimulus_indices)
%     search_range = stimulus_indices(i):min(stimulus_indices(i) + round(50e-3 * fs), length(power_envelope));
%     response_idx = search_range(find(power_envelope(search_range) > threshold_env, 1, 'first'));
%     if ~isempty(response_idx)
%         response_indices = [response_indices; response_idx];
%     end
% end

% % Trojúhelníková metoda detekce odpovědí
% response_indices = [];
% for i = 1:length(stimulus_indices)
%     segment = emg5000(stimulus_indices(i):min(stimulus_indices(i) + round(50e-3 * fs), length(emg5000)));
%     [amp, lmax] = findpeaks(abs(segment)); % Lokální maxima
%     lmax = lmax(amp > 0.05 * max(abs(segment))); % Maxima nad 5 % maxima MEP
%     if isempty(lmax)
%         continue; % Přeskakujeme, pokud není platné maximum
%     end
%     lmax = lmax(1); % První maximum
% 
%     % Trojúhelníková metoda
%     ax = 1; ay = abs(segment(1));       % Vrchol A
%     cx = lmax; cy = abs(segment(lmax)); % Vrchol C
%     S = [];
%     for k = 1:lmax
%         bx = k; by = abs(segment(k));
%         S = [S, 0.5 * abs((cx - ax) * (by - ay) - (cy - ay) * (bx - ax))];
%     end
%     [~, bx_max] = max(S);
%     response_indices = [response_indices; stimulus_indices(i) + bx_max];
% end


%% Výpočet latencí
if length(response_indices) ~= length(stimulus_indices)
    warning('Počet detekovaných odpovědí a stimulací se neshoduje.');
end

latencies = [(response_indices - stimulus_indices)]' / fs;
mean_latency = mean(latencies);

electrode_distance = 0.3;                           % Zadejte skutečnou vzdálenost mezi elektrodami v metrech
mean_velocity = electrode_distance / mean_latency;

%% Vykreslení výsledků
figure;
subplot(2,1,1);
plot(t, stim); hold on;
scatter(t(stimulus_indices), ones(size(stimulus_indices)) * stimulus_threshold, 'ro', 'filled');
title('Stimulační signál'); xlabel('Čas [s]'); ylabel('Amplituda'); grid on;
legend('Stimulační signál', 'Detekované stimuly');

subplot(2,1,2);
plot(t, emg5000); hold on;
scatter(t(response_indices), emg5000(response_indices), 'g', 'filled');
title('Motorické evokované odpovědi'); xlabel('Čas [s]'); ylabel('Amplituda [mV]'); grid on;
legend('EMG signál', 'Detekované odpovědi');

%% Výpis výsledků
% fprintf('Časy počátků stimulů: \n'); disp(t(stimulus_indices)');
% fprintf('Časy počátků odpovědí: \n'); disp(t(response_indices)');
fprintf('Latence [s]: \n'); disp(latencies');

fprintf('\n==== Výsledky detekce ====\n');
fprintf('Počet detekovaných stimulací: %d\n', length(stimulus_indices));
fprintf('Počet detekovaných odpovědí: %d\n', length(response_indices));
fprintf('Průměrná latence: %.5f s\n', mean_latency);
fprintf('Průměrná rychlost vedení nervu: %.2f m/s\n', mean_velocity);
fprintf('===========================\n\n');