%% Analisi FBG – Media per gruppo + plot a schermo 
clear; clc; close all;

dataFolder = 'C:\Users\Franc\Desktop\MatlabProve\14 Novembre';

savePlotsFolder = 'C:\Users\Franc\Desktop\MatlabProve\PlotsFinaliConMassimi';
if ~exist(savePlotsFolder, 'dir')
    mkdir(savePlotsFolder);
end


%% Elenco file
files = dir(fullfile(dataFolder, '*.mat'));
if isempty(files)
    error('Nessun file .mat trovato in %s', dataFolder);
end

%% Definizione gruppi 
groupNames = {'Prova-14-11-2025-A', 'Prova-14-11-2025-B', ...
              'Prova-14-11-2025-C', 'Prova2-14-11-2025-A', 'Prova2-14-11-2025-B', ...
              'Prova2-14-11-2025-C','Prova3-14-11-2025-A', 'Prova3-14-11-2025-B', ...
              'Prova3-14-11-2025-C'};

%% Parametri FFT
minFreq = 5;
maxFreq = 400;
hp_fc   = 1; % Hz

fprintf('Trovati %d file da ragruppare.\n', numel(files));

%% Preallocazione strutture
avgSpettro = struct();

for g = 1:numel(groupNames)
    thisGroup = groupNames{g};
    fprintf('\n=== Gruppo: %s ===\n', thisGroup);

    % Trova i file che appartengono a questo gruppo
    matches = contains({files.name}, thisGroup, 'IgnoreCase', true);
    groupFiles = files(matches);
    if isempty(groupFiles)
        warning('Nessun file per %s', thisGroup);
        continue;
    end

    %% Variabili cumulative per media (frequenza + tempo)
    allF = [];
    cumAmp  = cell(1,3);

    % --- AGGIUNTA: per media nel tempo ---
    timeRef = [];
    cumTime = cell(1,3);
    % -------------------------------------

    for fidx = 1:numel(groupFiles)
        inFile = fullfile(groupFiles(fidx).folder, groupFiles(fidx).name);
        S = load(inFile);

        if ~isfield(S,'time')
            warning('File %s corrotto.', groupFiles(fidx).name);
            continue;
        end

        t = S.time(:);
        data = [S.sensor1(:), S.sensor2(:), S.sensor3(:)];
        [nSamples, nSens] = size(data);

        %% Pre-processing
        dt = mean(diff(t));
        fs = 1/dt;
        hp = designfilt('highpassiir','FilterOrder',4, ...
            'HalfPowerFrequency',hp_fc,'SampleRate',fs);
        data_f = zeros(size(data));
        for k = 1:nSens
            x = detrend(data(:,k),'linear');
            data_f(:,k) = filtfilt(hp,x);
        end

       


        %% FFT con finestra di Hann
        N = size(data_f,1);
        N2 = floor(N/2);
        f = (0:N2-1)'*(fs/N);
        mask = f >= minFreq & f <= maxFreq;
        f_low = f(mask);

        w = hann(N);
        cg = mean(w);
        Amp = zeros(sum(mask), nSens);
        for k = 1:nSens
            X = fft(data_f(:,k).*w);
            Xpos = X(1:N2);
            Amp(:,k) = abs(Xpos(mask)) * (2/(N*cg));
        end

        %% Allineamento e somma per media in frequenza
        if isempty(allF)
            allF = f_low;
            for k = 1:nSens
                cumAmp{k} = Amp(:,k);
            end
        else
            for k = 1:nSens
                cumAmp{k}(:,end+1) = interp1(f_low, Amp(:,k), allF, 'linear', 'extrap');
            end
        end
    end

    %% Calcola media tra file del gruppo (FREQUENZA)
    for k = 1:3
        avgSpettro(g).Amp(:,k) = mean(cumAmp{k}, 2, 'omitnan');
    end
    avgSpettro(g).f = allF;
    avgSpettro(g).name = thisGroup;

    %% Calcola media tra file del gruppo (TEMPO) – AGGIUNTA
    if ~isempty(timeRef)
        avgSpettro(g).time = timeRef;
        avgSpettro(g).timeResp = zeros(numel(timeRef),3);
        for k = 1:3
            avgSpettro(g).timeResp(:,k) = mean(cumTime{k}, 2, 'omitnan');
        end
    end

    fprintf('  → %d file mediati.\n', numel(groupFiles));

    %% === PLOT COMBINATO DEI TRE SENSORI IN FREQUENZA (visibile) ===
    figure('Name',sprintf('Media FFT – %s', thisGroup),'NumberTitle','off');
    plot(allF, avgSpettro(g).Amp(:,1), 'r','LineWidth',1.2); hold on;
    plot(allF, avgSpettro(g).Amp(:,2), 'g','LineWidth',1.2);
    plot(allF, avgSpettro(g).Amp(:,3), 'b','LineWidth',1.2);
    grid on; xlim([0 maxFreq]);
    xlabel('Frequenza [Hz]'); ylabel('Ampiezza (rel.)');
    title(sprintf('Media FFT – %s', thisGroup), 'Interpreter','none');
    legend({'sensor1','sensor2','sensor3'}, 'Location','northeast');

    %% === PLOT COMBINATO DEI TRE SENSORI NEL TEMPO (visibile) ===
    if isfield(avgSpettro(g),'time')
        figure('Name',sprintf('Media risposta nel tempo – %s', thisGroup),'NumberTitle','off');
        % Se vuoi partire da t=0:
        t_plot = avgSpettro(g).time - avgSpettro(g).time(1);
        plot(t_plot, avgSpettro(g).timeResp(:,1), 'r','LineWidth',1.2); hold on;
        plot(t_plot, avgSpettro(g).timeResp(:,2), 'g','LineWidth',1.2);
        plot(t_plot, avgSpettro(g).timeResp(:,3), 'b','LineWidth',1.2);
        grid on;
        xlabel('Tempo [s]'); ylabel('Ampiezza filtrata (rel.)');
        title(sprintf('Media risposta nel tempo – %s', thisGroup), 'Interpreter','none');
        legend({'sensor1','sensor2','sensor3'}, 'Location','northeast');
    end
end

%% === PLOT CUMULATIVO (tutti i gruppi mediati) – FREQUENZA ===
colors = lines(numel(groupNames));
figure('Name','Confronto cumulativo – Tutte le prove','NumberTitle','off');
for g = 1:numel(groupNames)
    if ~isfield(avgSpettro(g),'f'), continue; end
    meanAmp = mean(avgSpettro(g).Amp, 2, 'omitnan');
    plot(avgSpettro(g).f, meanAmp, ...
        'LineWidth',1.3,'DisplayName',avgSpettro(g).name,'Color',colors(g,:)); hold on;
end
grid on; xlim([0 maxFreq]);
xlabel('Frequenza [Hz]'); ylabel('Ampiezza media (rel.)');
title('Confronto cumulativo – Medie per prova');
legend('Location','northeast');



%% === PLOT CUMULATIVO PER OGNI PROVA (A/B/C) – FREQUENZA ===
% Assumendo l'ordine:
% 1-3  -> Prova-14-11-2025- (A,B,C)
% 4-6  -> Prova2-14-11-2025- (A,B,C)
% 7-9  -> Prova3-14-11-2025- (A,B,C)

idxProva1 = [1 2 3];
idxProva2 = [4 5 6];
idxProva3 = [7 8 9];

% Prova 1: A, B, C
figure('Name','Confronto cumulativo – Prova 1','NumberTitle','off');
hold on;
for idx = idxProva1
    if ~isfield(avgSpettro(idx),'f'), continue; end
    meanAmp = mean(avgSpettro(idx).Amp, 2, 'omitnan');
    plot(avgSpettro(idx).f, meanAmp, 'LineWidth',1.3, ...
        'DisplayName', avgSpettro(idx).name);
end
grid on; xlim([0 maxFreq]);
xlabel('Frequenza [Hz]'); ylabel('Ampiezza media (rel.)');
title('Confronto cumulativo – Prova 1 (A-B-C)');
legend('Location','northeast');

% Prova 2: A, B, C
figure('Name','Confronto cumulativo – Prova 2','NumberTitle','off');
hold on;
for idx = idxProva2
    if ~isfield(avgSpettro(idx),'f'), continue; end
    meanAmp = mean(avgSpettro(idx).Amp, 2, 'omitnan');
    plot(avgSpettro(idx).f, meanAmp, 'LineWidth',1.3, ...
        'DisplayName', avgSpettro(idx).name);
end
grid on; xlim([0 maxFreq]);
xlabel('Frequenza [Hz]'); ylabel('Ampiezza media (rel.)');
title('Confronto cumulativo – Prova 2 (A-B-C)');
legend('Location','northeast');

% Prova 3: A, B, C
figure('Name','Confronto cumulativo – Prova 3','NumberTitle','off');
hold on;
for idx = idxProva3
    if ~isfield(avgSpettro(idx),'f'), continue; end
    meanAmp = mean(avgSpettro(idx).Amp, 2, 'omitnan');
    plot(avgSpettro(idx).f, meanAmp, 'LineWidth',1.3, ...
        'DisplayName', avgSpettro(idx).name);
end
grid on; xlim([0 maxFreq]);
xlabel('Frequenza [Hz]'); ylabel('Ampiezza media (rel.)');
title('Confronto cumulativo – Prova 3 (A-B-C)');
legend('Location','northeast');

% %% === PLOT NEL DOMINIO DEL TEMPO (taglio a 9 s se possibile) ===
% tMax = 9;  % [s] tempo massimo da plottare
% 
% fprintf('\n=== Plot nel dominio del tempo (taglio a %.1f s dove possibile) ===\n', tMax);
% 
% for g = 1:numel(groupNames)
%     thisGroup = groupNames{g};
%     fprintf('\n--- Gruppo: %s ---\n', thisGroup);
% 
%     % Trova i file che appartengono a questo gruppo
%     matches    = contains({files.name}, thisGroup, 'IgnoreCase', true);
%     groupFiles = files(matches);
%     if isempty(groupFiles)
%         warning('Nessun file per %s (per plot nel tempo).', thisGroup);
%         continue;
%     end
% 
%     for fidx = 1:numel(groupFiles)
% 
%         inFile = fullfile(groupFiles(fidx).folder, groupFiles(fidx).name);
% 
%         % Caricamento sicuro
%         try
%             S = load(inFile);
%         catch ME
%             warning('Impossibile caricare il file %s: %s', groupFiles(fidx).name, ME.message);
%             continue;
%         end
% 
%         % Controllo campi necessari
%         if ~isfield(S,'time') || ~isfield(S,'sensor1') || ...
%            ~isfield(S,'sensor2') || ~isfield(S,'sensor3')
%             warning('File %s: mancano time/sensor1/sensor2/sensor3. Salto.', groupFiles(fidx).name);
%             continue;
%         end
% 
%         % Estrazione segnali grezzi
%         t  = S.time(:);
%         s1 = S.sensor1(:);
%         s2 = S.sensor2(:);
%         s3 = S.sensor3(:);
% 
%         % Allineo alla stessa lunghezza minima, per sicurezza
%         minLen = min([numel(t), numel(s1), numel(s2), numel(s3)]);
%         t  = t(1:minLen);
%         s1 = s1(1:minLen);
%         s2 = s2(1:minLen);
%         s3 = s3(1:minLen);
% 
%         data = [s1, s2, s3];
% 
%         % Taglio a 9 s se il segnale è più lungo
%         if t(end) > tMax
%             maskTime = t <= tMax;
%         else
%             maskTime = true(size(t));
%         end
% 
%         tPlot    = t(maskTime);
%         dataPlot = data(maskTime,:);
% 
%         % Controllo minimo numero di campioni
%         if numel(tPlot) < 5
%             warning('File %s: meno di 5 campioni dopo il taglio. Salto.', groupFiles(fidx).name);
%             continue;
%         end
% 
%         % Pre-processing (come nella FFT: detrend + high-pass)
%         dt = mean(diff(tPlot));
%         fs = 1/dt;
% 
%         try
%             hp = designfilt('highpassiir','FilterOrder',4, ...
%                 'HalfPowerFrequency',hp_fc,'SampleRate',fs);
% 
%             data_f = zeros(size(dataPlot));
%             for k = 1:3
%                 x = detrend(dataPlot(:,k),'linear');
%                 data_f(:,k) = filtfilt(hp,x);
%             end
%         catch ME
%             warning('File %s: errore nel filtraggio: %s. Salto.', ...
%                 groupFiles(fidx).name, ME.message);
%             continue;
%         end
% 
%         % Plot nel tempo (3 sensori sovrapposti)
%         figure('Name',sprintf('Segnali nel tempo - %s - file %d', ...
%                               thisGroup, fidx), ...
%                'NumberTitle','off');
%         plot(tPlot, data_f(:,1), 'r', 'LineWidth', 1.1); hold on;
%         plot(tPlot, data_f(:,2), 'g', 'LineWidth', 1.1);
%         plot(tPlot, data_f(:,3), 'b', 'LineWidth', 1.1);
%         grid on;
%         xlabel('Tempo [s]');
%         ylabel('Ampiezza (filtrata, rel.)');
%         titStr = sprintf('Segnali nel tempo (max %.1f s) - %s - %s', ...
%                          tMax, thisGroup, groupFiles(fidx).name);
%         title(titStr, 'Interpreter','none');
%         legend({'sensor1','sensor2','sensor3'}, 'Location','best');
% 
%     end
% end
