%% GASTROGRAM_DEMO.M
% Skript pro analýzu elektrogastrogramu (EGG)
% Cíle: zobrazení v časové oblasti, spektrální analýza, klouzavý průměr, decimace
% Signál: 2-kanálový záznam, 1. část před jídlem, 2. část po jídle

clear, close all

%% Načtení dat
load('cele.txt'); % 2 kanály
fs = 5;           % vzorkovací frekvence v Hz
x = cele;

% Časový úsek: 0–20 min (3000 + 3000 vzorků)
start1 = 1; stop1 = 3000;
start2 = 3001; stop2 = 6000;
time = (0:1/fs:(stop2-start1)/fs)/60; % čas v minutách

%% 1. Časová oblast – původní signál
figure(1)
for k = 1:2
    subplot(2,1,k)
    plot(time, x(start1:stop2,k))
    xlabel('Čas [min]')
    ylabel(['Kanál ' num2str(k)])
    axis([0 20 -0.6 0.6])
end
sgtitle('EGG signál (před jídlem a po jídle), fs = 5 Hz')
text(6, -0.5, 'Před jídlem       –––––––––––––––       Po jídle')

%% 2. Frekvenční oblast – spektrogramy
figure(2)
for k = 1:2
    subplot(2,1,k)
    specgram(x(:,k), 2048, fs*60, 256, 250) % jednotky [cpm]
    xlabel('Čas [min]'); ylabel('frekvence [cpm]')
end
sgtitle('Spektrogramy EGG signálu (2 kanály)')

%% 3. Klouzavý průměr – 30bodový filtr
L = 30; % délka okna
y = filter(ones(1,L)/L, 1, x); % klouzavý průměr

figure(3)
for k = 1:2
    subplot(2,1,k)
    plot(time, y(start1:stop2,k))
    xlabel('Čas [min]')
    ylabel(['Kanál ' num2str(k)])
    axis([0 20 -0.2 0.2])
end
sgtitle('EGG signál – klouzavý průměr (30 vzorků)')

%% 4. Decimace – snížení fs z 5 Hz na 1 Hz
M = 5;
z(:,1) = decimate(x(:,1), M);
z(:,2) = decimate(x(:,2), M);
fs_dec = fs / M;
time_dec = (0:1/fs_dec:(length(z)-1))/fs_dec / 60; % čas v minutách

figure(4)
for k = 1:2
    subplot(2,1,k)
    plot(time_dec, z(:,k))
    xlabel('Čas [min]')
    ylabel(['Kanál ' num2str(k)])
    axis([0 20 -0.25 0.25])
end
sgtitle('Decimovaný EGG signál (fs = 1 Hz)')

%% 5. Parametrická spektra (LPC) – 1. kanál před/po jídle
rad = 50; % řád LPC modelu
figure(5)
segments = [31 300; 331 600]; % cca 5 min

for k = 1:2
    for i = 1:2
        idx = (k-1)*2 + i;
        s = z(segments(i,1):segments(i,2), k);
        s = s - mean(s);
        a = lpc(s, rad);
        [H,~] = freqz(1, a, 300);
        subplot(2,2,idx)
        plot(0:0.1:29.9, abs(H)), grid on
        axis([0 10 0 12.5])
        xlabel('Frekvence [cpm]')
        ylabel(['Kanál ' num2str(k)])
        if i == 1
            stav = 'Před jídlem';
        else
            stav = 'Po jídle';
        end

        f_axis = linspace(0, fs_dec/2, length(H))*60; % převod na cpm
        subplot(2,2,idx)
        plot(f_axis, abs(H)); hold on; grid on
        axis([0 10 0 12.5])
        xlabel('f (cpm)'), ylabel('S(f)')

        % Detekce špičky v okolí 3 cpm
        f_range = (f_axis >= 2 & f_axis <= 4);
        [~, ipeak] = max(abs(H(f_range)));
        peak_freq = f_axis(f_range);
        f_peak = peak_freq(ipeak);
        plot(f_peak, abs(H(f_axis == f_peak)), 'ro', 'MarkerSize', 8)
        legend(sprintf('%.1f cpm', f_peak), 'Location', 'northeast')

        title([stav ' – Kanál ' num2str(k)])
    end
end
sgtitle('Parametrická LPC spektra (fs = 1 Hz)')

%% 6. Periodogram – Welchův průměr
figure(6)
window = boxcar(100);
noverlap = 90;
for k = 1:2
    for i = 1:2
        idx = (k-1)*2 + i;
        s = z(segments(i,1):segments(i,2), k);
        subplot(2,2,idx)
        [Pxx, f] = pwelch(s, window, noverlap, 2048, fs_dec);
        f_cpm = f*60;
        subplot(2,2,idx)
        plot(f_cpm, 10*log10(Pxx)); hold on; grid on
        axis([0 10 -40 -5])
        xlabel('f (cpm)'), ylabel('S(f) [dB]')

        % Detekce špičky v okolí 3 cpm
        f_range = (f_cpm >= 2 & f_cpm <= 4);
        [~, ipeak] = max(Pxx(f_range));
        peak_freq = f_cpm(f_range);
        f_peak = peak_freq(ipeak);
        plot(f_peak, 10*log10(Pxx(f_cpm == f_peak)), 'ro', 'MarkerSize', 8)
        legend(sprintf('%.1f cpm', f_peak), 'Location', 'northeast')

        axis([0 10 -40 -5])
        xlabel('f [cpm]')
        ylabel(['Kanál ' num2str(k)])
        if i == 1
            stav = 'Před jídlem';
        else
            stav = 'Po jídle';
        end
        title([stav ' – Kanál ' num2str(k)])      
    end
end
sgtitle('Periodogramy EGG signálu (Welch, boxcar)')