% Glotogram_demo_final.m
% Skript – EGG signál, extrémy, práh, CQ

clear; close all; clc;

% --- Načtení dat ---
load sig01.mat  % musí obsahovat: mic, egg, fs
t = (0:length(egg)-1)/fs;

% --- Filtrace EGG ---
[b,a] = butter(4, 20/(fs/2), 'high');
egg_f = filtfilt(b,a,egg);

% --- Detekce extrémů ---
minDist = round(0.002 * fs);  % minimální vzdálenost 2 ms
[maxvals, maxlocs] = findpeaks(egg_f, 'MinPeakDistance', minDist, 'MinPeakHeight', 0);
[minvals, minlocs] = findpeaks(-egg_f, 'MinPeakDistance', minDist, 'MinPeakHeight', 0);
minvals = -minvals;
t_max = t(maxlocs);
t_min = t(minlocs);

% --- Obálky a práh ---
upp = interp1(t_max, maxvals, t, 'pchip', 'extrap');
low = interp1(t_min, minvals, t, 'pchip', 'extrap');
th = low + 0.35 * (upp - low);

% --- CQ prahování ---
egg_f = egg_f(:); th = th(:);
above = egg_f > th;
cross_up = find(diff(above) == 1);
cross_down = find(diff(above) == -1);
t1_all = t(cross_up);
t2_all = t(cross_down);
pairs = [];
i = 1; j = 1;
while i <= length(t1_all) && j <= length(t2_all)
    if t1_all(i) < t2_all(j)
        pairs(end+1, :) = [t1_all(i), t2_all(j)];
        i = i + 1; j = j + 1;
    else
        j = j + 1;
    end
end
t1 = pairs(:,1);
t2 = pairs(:,2);
tau = t2 - t1;
T = [diff(t1); NaN];
CQ = tau ./ T;
CQ(CQ < 0 | CQ > 1) = NaN;

% --- CQ derivace ---
egg_f = filtfilt(b, a, egg);  % znovu filtr
degg = diff(egg_f); degg(end+1) = degg(end);
degg_smooth = movmean(degg, 5);
[maxD, maxLoc] = findpeaks(degg_smooth, 'MinPeakDistance', minDist, 'MinPeakHeight', 0.001);
[minD, minLoc] = findpeaks(-degg_smooth, 'MinPeakDistance', minDist, 'MinPeakHeight', 0.001);
minD = -minD;
t_maxD = t(maxLoc);
t_minD = t(minLoc);
pairsD = [];
i = 1; j = 1;
while i <= length(t_maxD) && j <= length(t_minD)
    if t_maxD(i) < t_minD(j)
        pairsD(end+1, :) = [t_maxD(i), t_minD(j)];
        i = i + 1; j = j + 1;
    else
        j = j + 1;
    end
end
t1D = pairsD(:,1);
t2D = pairsD(:,2);
tauD = t2D - t1D;
T_D = [diff(t1D); NaN];
CQ_der = tauD ./ T_D;
CQ_der(CQ_der < 0 | CQ_der > 1) = NaN;

% --- Detail pro zobrazení ---
ix_zoom = (t >= 1 & t <= 1.1);
ix_CQ = (t1 >= 1 & t1 <= 1.1);
ixD = (t1D >= 1 & t1D <= 1.1);
t_zoom = t(ix_zoom);
egg_zoom = egg(ix_zoom);
egg_f_zoom = egg_f(ix_zoom);
degg_zoom = degg_smooth(ix_zoom);
upp_zoom = upp(ix_zoom);
low_zoom = low(ix_zoom);
th_zoom = th(ix_zoom);
CQ_zoom = CQ(ix_CQ);

%% --- Figure 1: Celkový přehled ---
figure;
subplot(3,1,1); plot(t, mic, 'b'); ylabel('Mic [mV]'); title('Mikrofon'); grid on;
subplot(3,1,2); plot(t, egg, 'k'); ylabel('EGG'); title('EGG signál'); grid on;
subplot(3,1,3); plot(t1, CQ, 'b'); hold on; plot(t1D, CQ_der, 'm');
legend('Prahování','Derivace'); ylabel('CQ [-]');
xlabel('t [s]'); title('Kontaktní kvocient CQ (celý signál)'); ylim([0 1]); grid on;

%% --- Figure 2: Detail mikrofon a EGG ---
figure;
subplot(2,1,1); plot(t_zoom, mic(ix_zoom), 'b'); ylabel('Mic [mV]');
title('Detail mikrofon'); grid on; axis tight;
subplot(2,1,2); plot(t_zoom, egg_zoom, 'k'); ylabel('EGG'); title('Detail EGG'); grid on; axis tight;

%% --- Figure 3: Prah a CQ ---
figure;
subplot(3,1,1); plot(t_zoom, egg_f_zoom, 'k'); hold on;
plot(t_zoom, upp_zoom, 'r'); plot(t_zoom, low_zoom, 'b'); plot(t_zoom, th_zoom, 'm');
legend('EGG_f', 'Horní obálka', 'Dolní obálka', '35 % práh');
xlabel('t [s]'); ylabel('Amplituda'); title('EGG_f s obálkami a prahem (výřez 100 ms)'); grid on; axis tight;

subplot(3,1,2);
plot(t1(ix_CQ), CQ(ix_CQ), 'b'); ylabel('CQ [-]');
xlabel('t [s]'); title('Kontaktní kvocient – prahování'); ylim([0 1]); grid on; axis([1 1.1 0 1]);

subplot(3,1,3);
binary_state = zeros(size(t_zoom));
for k = 1:length(t1)
    if t2(k) < t_zoom(1) || t1(k) > t_zoom(end)
        continue;
    end
    idx = t_zoom >= t1(k) & t_zoom <= t2(k);
    binary_state(idx) = 1;
end
plot(t_zoom, binary_state, 'k', 'LineWidth', 1.5);
ylim([-0.1 1.1]);
xlabel('t [s]'); ylabel('Stav');
title('Stav hlasivek: 1 = kontakt, 0 = otevřeno'); grid on; axis tight;

%% --- Figure 4: Derivace a CQ ---
figure;
subplot(4,1,1);
plot(t_zoom, egg_zoom, 'k'); hold on;
plot(t1D(ixD), interp1(t, egg_f, t1D(ixD)), 'r.', 'MarkerSize', 10);
plot(t2D(ixD), interp1(t, egg_f, t2D(ixD)), 'm.', 'MarkerSize', 10);
ylabel('EGG'); title('Detail EGG s detekcí'); grid on; axis tight;

subplot(4,1,2);
plot(t_zoom, degg_zoom, 'b'); hold on;
plot(t_maxD(t_maxD >= 1 & t_maxD <= 1.1), degg_smooth(maxLoc(t_maxD >= 1 & t_maxD <= 1.1)), 'ro');
plot(t_minD(t_minD >= 1 & t_minD <= 1.1), degg_smooth(minLoc(t_minD >= 1 & t_minD <= 1.1)), 'mo');
ylabel('dEGG'); title('Derivace s extrémy'); grid on; axis tight;

subplot(4,1,3);
plot(t1D(ixD), CQ_der(ixD), 'b'); ylabel('CQ [-]');
title('CQ z derivace'); ylim([0 1]); grid on; axis([1 1.1 0 1]);

subplot(4,1,4);
binState = zeros(size(t));
for i = 1:length(t1D)
    if t2D(i) > t1D(i)
        binState(t >= t1D(i) & t <= t2D(i)) = 1;
    end
end
plot(t_zoom, binState(ix_zoom), 'k', 'LineWidth', 1.5);
xlabel('t [s]'); ylabel('Stav');
title('Stav hlasivek: 1 = kontakt, 0 = otevřeno'); ylim([1 2.2]);
yticks([1 2]); yticklabels({'otevřeno','zavřeno'}); grid on; axis tight;