Return to Question

edited title
Link
edited May 8 at 9:46
toolic
16.4k
6
29
221
Where does My MATLAB Gait Analysis Script mistake and how can I correct it?

Source Link
asked May 8 at 9:42
elena conidi
Where does My MATLAB Gait Analysis Script mistake and how can I correct it?

Could someone kindly help me review and fix any errors or inefficiencies in this MATLAB script? I’m trying to make sure it’s clean, robust, and optimized.
Brief description of what the code does:
This MATLAB script processes biomechanical data from different subject groups (e.g., Healthy, PD, Stroke, TBI) during gait tasks. It loads and combines motion capture data, filters for valid entries, computes spatiotemporal gait metrics (e.g., stride length, speed, duration), calculates variability measures like coefficient of variation (CV), phase coordination index (PCI), and gait asymmetry (GA), and conducts statistical analysis (ANOVA and Tukey HSD) to compare gait characteristics across groups. It also includes pre-processing for sensor signals and attempts to correlate gait events with signal features (e.g., heel strikes via accelerometer peaks).
clear

path = 'C:\Users\franc\Desktop\BIOMECCANICA';

pop = dir(path); pop = {pop(3:end).name};

[indx,~] = listdlg('ListString',pop);
pop_true = pop(indx);


file_to_load = {'ROOT' 'INDIP' 'LWRIST' 'RWRIST' 'LFOOT' 'RFOOT' 'LBACK' 'STERN' 'FHEAD'};
[indx,~] = listdlg('ListString',file_to_load);
file_true = file_to_load(indx);
TAB = table;
TAB1 = table;

%Carico tutta la tabella TAB
for po = 1:length(pop_true)

    files = dir([path '\' char(pop_true(po))]); files = {files(3:end).name};
    files_true = files(contains(files,file_true));

    for fi = 1:length(files_true)

        temp = load([path '\' char(pop_true(po)) '\' char(files_true(fi))]);
        fld = char(fieldnames(temp)); ind = strfind(fld,'_');
        if contains(fld,'INDIP') == 1 || contains(fld,'ROOT') == 1
            temp = temp.(fld);
        else
            temp = table(temp.(fld));
            temp.Properties.VariableNames = {fld(ind+1:end)};
        end
        TAB1 = [TAB1 temp];

    end

    TAB = [TAB; TAB1];

    TAB1 = table;
end

TAB = removevars(TAB,'b2');
varnames = TAB.Properties.VariableNames;
varnames(contains(varnames,'b1')) = {'INDIP'};

TAB.Properties.VariableNames = varnames;

TAB =  movevars(TAB,'GROUP','Before',1); TAB = movevars(TAB,{'SUBJECT' 'TASK' 'TRIAL' 'ORIGINAL PATH' 'AGE' 'MASS' 'HEIGHT' 'INDIP'},'After','GROUP');

for ii = 1:size(TAB,1)
    TAB.HEIGHT{ii} = cell2mat(TAB.HEIGHT(ii))*0.53;
end


c = contains(TAB.("ORIGINAL PATH"),'FAST');
TAB = TAB(c,:);

timing = repmat({'T0'},size(TAB,1),1);

TAB.TIMING = timing;
TAB = movevars(TAB,'TIMING','Before','TASK');

%% CLEANING
for ii = 1:size(TAB,1)

    mwb = TAB.INDIP(ii).MicroWB;

    ind(ii) =   isempty(mwb);
end

ind = logical(ind);
TAB(ind,:) = [];
ind = [];

for ii = 1:size(TAB,1)

    mwb = TAB.FHEAD(ii).acc;

    ind(ii) =   isempty(mwb);
end

ind = logical(ind);
TAB(ind,:) = [];
ind = [];


for ii = 1:size(TAB,1)

    mwb = TAB.LBACK(ii).acc;
    ind(ii) =   isempty(mwb);
end

ind = logical(ind);
TAB(ind,:) = [];
ind = [];


for ii = 1:size(TAB,1)

    mwb = TAB.STERN(ii).acc;

    ind(ii) =   isempty(mwb);
end

ind = logical(ind);
TAB(ind,:) = [];
ind = [];


% %% GET SPATIOTEMPORAL
for ii = 1:size(TAB,1)

    mwb = TAB.INDIP(ii).MicroWB;

    TAB.STRIDE_FQ(ii) = {nanmean(cell2mat({mwb.StrideFrequency}))};


    spa = [];

    for jj = 1:size(mwb,2)
        temp = mwb(jj).Stride_Length;
        spa = [spa temp];
    end

    TAB.STRIDE_LENGTH(ii) = {(spa)};


    spa = [];

    for jj = 1:size(mwb,2)
        temp = mwb(jj).Stride_Speed;
        spa = [spa temp];
    end


    TAB.STRIDE_SPEED(ii) = {(spa)};


    spa = [];

    for jj = 1:size(mwb,2)
        temp = mwb(jj).Stride_Duration;
        spa = [spa temp];
    end

    TAB.STRIDE_DURATION(ii) = {(spa)};

    spa = [];

    for jj = 1:size(mwb,2)
        temp = mwb(jj).Stride_Length;
        spa = [spa temp];
    end

    TAB.STRIDE_LENGTH(ii) = {(spa)};

    spa = [];

    for jj = 1:size(mwb,2)
        temp = mwb(jj).Stance_Duration;
        spa = [spa temp];
    end

    TAB.STANCE_DURATION(ii) = {(spa)};

    spa = [];

    for jj = 1:size(mwb,2)
        temp = mwb(jj).Stance_Length;
        spa = [spa temp];
    end

    TAB.STANCE_LENGTH(ii) = {(spa)};

    spa = [];

    for jj = 1:size(mwb,2)
        temp = mwb(jj).Stance_Speed;
        spa = [spa temp];
    end

    TAB.STANCE_SPEED(ii) = {(spa)};


    spa = [];

    for jj = 1:size(mwb,2)
        temp = mwb(jj).Swing_Duration;
        spa = [spa temp];
    end

    TAB.SWING_DURATION(ii) = {(spa)};

    spa = [];

    for jj = 1:size(mwb,2)
        temp = mwb(jj).Swing_Length;
        spa = [spa temp];
    end

    TAB.SWING_LENGTH(ii) = {(spa)};

    spa = [];

    for jj = 1:size(mwb,2)
        temp = mwb(jj).Swing_Speed;
        spa = [spa temp];
    end

    TAB.SWING_SPEED(ii) = {(spa)};

    spa = [];

    for jj = 1:size(mwb,2)
        temp = mwb(jj).SingleSupport_Duration;
        spa = [spa temp];
    end

    TAB.SS_DURATION(ii) = {(spa)};

    spa = [];

    for jj = 1:size(mwb,2)
        temp = mwb(jj).DoubleSupport_Duration;
        spa = [spa temp];
    end

    TAB.DS_DURATION(ii) = {(spa)};

    spa = [];

    for jj = 1:size(mwb,2)
        temp = mwb(jj).Stride_InitialContacts;
        spa = [spa; temp];
    end

        [r,c] = find(isnan(spa)); r = unique(r); spa(r,:) = [];
    TAB.CONTANCTS(ii) = {(spa)};

    spa = [];

    for jj = 1:size(mwb,2)
        temp = mwb(jj).FinalContact_LeftRight;
        spa = [spa temp];
    end

    %     spa(:,r) = [];
    TAB.SIDE(ii) = {(spa)};


end

Dividiamo in TAB10 e TAB8
% Selezioniamo le righe in cui TASK è '10MWT' (cammino rettilineo)
TAB10 = TAB(strcmp(TAB.TASK, '10MWT'), :)

% Selezioniamo le righe in cui TASK è '8MWT' (cammino a 8)
TAB8 = TAB(strcmp(TAB.TASK, 'Fo8WT'), :)

Creo una tabella ridotta (TABreduced), rispetto al task '10MTW' (TAB10), contenente le prime 20 righe per ogni malattia
TABreduced = table; % Inizializzo una tabella vuota per il risultato

% Filtro i gruppi per ciascun tipo di malattia solo per TASK = '10MWT'+
groups = {'Healthy_Control', 'PD', 'STROKE', 'TBI'};
start_rows = [1, 254, 413, 602]; % Indici di inizio per ogni gruppo
end_rows = [253, 412, 601, 680]; % Indici di fine per ogni gruppo

for i = 1:length(groups)
    % Seleziono le righe per il gruppo corrente da TAB10
    group_data = TAB10(strcmp(TAB10.GROUP, groups{i}), :);

    % Prendo solo le prime 15 righe 
    num_rows = min(15, size(group_data, 1));   %prendo al massimo 15 righe del gruppo, ma se il gruppo ha meno di 15 righe, allora prenderai tutte le righe disponibili. Quindi, se il gruppo ha meno di 15 righe, il codice prende tutte le righe, altrimenti prende solo le prime 15 righe.
    group_data_reduced = group_data(1:num_rows, :);    %Seleziona le righe dalla 1 alla num_rows (ad esempio, se ci sono 10 righe, prende le righe da 1 a 10).

    % Aggiungo i dati ridotti a TABreduced
    TABreduced = [TABreduced; group_data_reduced];

    end


Calcolo, per tutte le malattie, il CV.
CV = 100SD∕mean where mean is the mean step length and SD is the standard deviation over the entire step length for each subject. The higher the CV, the higher is the variability of the step length. Un valore basso del CV indica una maggiore regolarità (i passi sono più simili tra loro), mentre un valore alto indica una maggiore variabilità tra gli stride.
stride=l'intervallo di movimento compreso tra due appoggi consecutivi dello stesso piede.

CV_individuale = zeros(height(TABreduced),1); 

for ii = 1:height(TABreduced)
    stride = TABreduced.STRIDE_LENGTH{ii}; 
    CV_ind = 100 * std(stride,'omitmissing') / mean(stride,'omitmissing');
    CV_individuale(ii) = CV_ind;
end

% Aggiungo la colonna alla tabella
TABreduced.CV_individuale = CV_individuale;

CV della stride_length medio per gruppo
CV_Healthy = mean(TABreduced.CV_individuale(1:15))
CV_PD = mean(TABreduced.CV_individuale(16:30))
CV_ST = mean(TABreduced.CV_individuale(31:45))
CV_TBI= mean(TABreduced.CV_individuale(46:60))

PCI individuale per tutti i pazienti

PCI_individuale = zeros(height(TABreduced),1); %Creo nuova colonna

for ii = 1:height(TABreduced)


ic_times = TABreduced.INDIP(ii).ContinuousWalkingPeriod.InitialContact_Event;
ic_side  = TABreduced.INDIP(ii).ContinuousWalkingPeriod.InitialContact_LeftRight;

ic_times = ic_times(:);
ic_side = ic_side(:);

right_idx = [];
left_idx = [];

for c = 1:length(ic_side)
% Separazione in destra e sinistra
    if ic_side(c) == "Right"
    right_idx = [right_idx c];
else
    left_idx = [left_idx c];
    end
end
% Assumo A = destro, B = sinistro
% Calcolo fasi φi tra un heel strike sinistro tra due destri (o viceversa)
phi = [];

for i = 1:length(right_idx)-1
    A1 = ic_times(right_idx(i));
    A2 = ic_times(right_idx(i+1));
    
    % Cerco heel strike sinistro che cade tra A1 e A2
    mid_L_idx = find(ic_times(left_idx) > A1 & ic_times(left_idx) < A2);
    
    if ~isempty(mid_L_idx)
        B = ic_times(left_idx(mid_L_idx(1)));
        phi_i = ((B - A1) / (A2 - A1)) * 360;
        phi(end+1) = phi_i;
    end
% PCI = CV_phi + ME_phi
mean_phi = mean(phi);
std_phi = std(phi);
CV_phi = (std_phi / mean_phi) * 100;
ME_phi = mean(abs(phi - 180));

PCI = CV_phi + ME_phi;

PCI_individuale(ii) = PCI;

    end
end

% Aggiungo la colonna alla tabella
TABreduced.PCI_individuale = PCI_individuale;


% Crea il vettore dei gruppi (20 per ciascuno)
%analisi statistica
n = 15;
gruppi = [repmat({'Healthy_Control'}, n, 1); ...
          repmat({'PD'}, n, 1); ...
          repmat({'STROKE'}, n, 1); ...
          repmat({'TBI'}, n, 1)];
gruppi = categorical(gruppi);

% Esempio: PCI = dati a caso (solo per test)
PCI = TABreduced.PCI_individuale;  

% ANOVA a una via
[p, tbl, stats] = anova1(PCI, gruppi);
% Test post-hoc Tukey HSD
figure;
[c,~,~,gnames] = multcompare(stats);
%p-value per confermare differenze viste nel multcompare

tbl = array2table(c, 'VariableNames', ...
    ["GroupA", "GroupB", "LowerLimit", "A_B", "UpperLimit", "P_value"]);

tbl.GroupA = gnames(tbl.GroupA);
tbl.GroupB = gnames(tbl.GroupB)

GAIT ASIMMETRY (GA)
GA_individuale = zeros(height(TABreduced),1); % create new column

for ii = 1:height(TABreduced)
    swing_times = TABreduced.INDIP(ii).ContinuousWalkingPeriod.Swing_Duration;
    swing_side  = TABreduced.INDIP(ii).ContinuousWalkingPeriod.InitialContact_LeftRight;

    % Ensure these are columns
    swing_times = swing_times(:);
    swing_side = swing_side(:);

    right_idx_swing = [];
    left_idx_swing = [];

    % Separate right and left swings
    for c = 1:length(swing_side)
        if swing_side(c) == "Right"
            right_idx_swing = [right_idx_swing c];
        else
            left_idx_swing = [left_idx_swing c];
        end
    end

    % Assumptions: C = right, D = left
    GA = []; % Initialize GA for this individual

 if length(right_idx_swing) && ~isempty(left_idx_swing)
        for i = 1:(length(right_idx_swing)-1)
            C1 = swing_times(right_idx_swing(i));    
            C2 = swing_times(right_idx_swing(i+1));

    mid_L_idx_swing = find(swing_times(left_idx_swing) > C1 & swing_times(left_idx_swing) < C2);
        if ~isempty(mid_L_idx_swing) && all(left_idx_swing(mid_L_idx_swing) <= length(swing_times))
    D = swing_times(left_idx_swing(mid_L_idx_swing(1))); 
        end
            % Calculate GA for this pair (use log to calculate asymmetry)
            GA = 100 * log(mean(C1) / mean(D));

        end
    end
end

GA_individuale(ii) = GA; 
TABreduced.GA_individuale = GA_individuale;
    

HEALTHY
%% calcolo Healthy solo per TABreduced, prendendo i dati da 1 a 20 (tutti gli Healthy per il cammino rettilineo nella tabella ridotta)

for ii = 1:15
    lb_acc_Healthy10 = TABreduced.LBACK(ii).acc; lb_acc_Healthy10 = fillmissing(lb_acc_Healthy10,'spline'); lb_acc_Healthy10 = bwfilt(lb_acc_Healthy10,2,128,10,'low');
    lb_gyr_Healthy10 = TABreduced.LBACK(ii).gyro; lb_gyr_Healthy10 = fillmissing(lb_gyr_Healthy10,'spline'); lb_gyr_Healthy10 = bwfilt(lb_gyr_Healthy10,2,128,6,'low');
    % 
    % st_acc = TAB10.STERN(ii).acc; st_acc = fillmissing(st_acc,'spline'); st_acc = bwfilt(st_acc,2,128,10,'low');
    % st_gyr = TAB10.STERN(ii).gyro; st_gyr = fillmissing(st_gyr,'spline'); st_gyr = bwfilt(st_gyr,2,128,6,'low');
    % 
    % fh_acc = TAB10.FHEAD(ii).acc; fh_acc = fillmissing(fh_acc,'spline'); fh_acc = bwfilt(fh_acc,2,128,10,'low');
    % fh_gyr = TAB10.FHEAD(ii).gyro; fh_gyr = fillmissing(fh_gyr,'spline'); fh_gyr = bwfilt(fh_gyr,2,128,6,'low');
    % 

    contacts = TABreduced.CONTANCTS{ii};
    side = TABreduced.SIDE{ii};

    [r,c] = (find(isnan(contacts))); r = unique(r); contacts(r,:) = []; %TAB.CONTANCTS{ii} = contacts;
    side(:,r) = []; %TAB.SIDE{ii} = side;
       

    for co = 1:length(contacts)
% if

    end
end

DST (Double support Time) ->è un indice temporale che permette di discriminare la gravità del parkinson, non è un indice di variabilità motoria. pazienti con patologie neurologiche tipo parkinson hanno un DST piu lungo perche per paura di cadere rimangono piu tempo con i piedi a terra a differenza degli Healthy che sono piu sicuri durante la camminata
IL CoV di DST POTREBBE ESSERE CALCOLATO
% DefiniscO un threshold per rilevare i picchi nei segnali di accelerazione
threshold = 0.5;

% TrovO i picchi nel segnale di accelerazione
acc_z = lb_acc_Healthy10(:,3); % Prendo solo la componente Z
[peaks, locs] = findpeaks(abs((acc_z)), 'MinPeakHeight', threshold, 'MinPeakDistance', 200); 

% Aggiungo un ciclo per associare i picchi con i contatti
heel_strikes = {}; % Inizializzo una lista per i momenti di heel strike

% Associo i picchi ai momenti di contatto
for i = 1:length(ic_times)
    contact_time = ic_times(i); % Tempo del contatto
    side = ic_side(i);          % Sinistro o destro

    % Trova i picchi che avvengono vicino al tempo del contatto (0.8s)
    %dal grafico abbiamo visto la distanza tra una linea tratteggiata, che rappresenta il contact time, el'altra
 idx = find(locs/128 >= contact_time - 1 & locs/128 <= contact_time + 1);
t = (1:length(acc_z)) / 128; % asse dei tempi
plot(t, acc_z); hold on;
plot(locs/128, acc_z(locs), 'rv'); % picchi
xline(contact_time, 'g--'); % tempo del contatto
legend('acc_z','picchi','contact time');
    if ~isempty(idx)
        for j = 1:length(idx)
             
            heel_strikes{end+1,1} = contact_time;
            heel_strikes{end,  2} = side; % stringa 'Right' o 'Left'
            heel_strikes{end,  3} = locs(idx(j))/128; % tempo del picco
            heel_strikes{end,  4} = peaks(idx(j));
        end
    end
end
%nel workspace, in heel_strikes troviamo: prima colonna è il contact time
%(quando il piede tocca terra), la seconda (il piede), la terza (istante in
%cui si rileva un picco, che corrisponde a quando uno dei due piedi tocca
%terra), la quarta rappresenta il valore del picco dell'accelerazione del
%piede che tocca terra, vedi foto gruppo


rs = []; % Right side
ls = []; % Left side
% Scorri tutte le righe di heel_strikes
for i = 1:size(heel_strikes, 1)
    side = heel_strikes{i, 2};       % 'Right' o 'Left'
    contact_time = heel_strikes{i, 1}; % tempo di contatto

    % Aggiungi il contact_time al vettore corrispondente
    if strcmp(side, 'Right')
        rs(end+1, 1) = contact_time;
    elseif strcmp(side, 'Left')
        ls(end+1, 1) = contact_time;
    end
end
% STEP 2: CALCOLO DST PER OGNI CICLO
numCycles = min(length(rs), length(ls));
DST = zeros(numCycles-1,1);
for k = 2:numCycles

    % sovrapposizione: da max dei due contatti precedenti a min dei due successivi
    DST(k-1) = min(rs(k), ls(k)) - max(rs(k-1), ls(k-1));
end

% STEP 3: AGGIUNTA ALLA TABELLA
TABreduced.DST{ii}     = DST;                       % vettore DST per trial
TABreduced.MeanDST(ii) = mean(DST,'omitnan');      % media DST
TABreduced.CVDST(ii)   = 100 * std(DST,'omitnan') / TABreduced.MeanDST(ii);  % coeff. di variazione



PD
%% calcolo PD solo per TAB10

for ii = 16:30
    lb_acc_PD10 = TAB10.LBACK(ii).acc; lb_acc_PD10 = fillmissing(lb_acc_PD10,'spline'); lb_acc_PD10 = bwfilt(lb_acc_PD10,2,128,10,'low');
    lb_gyr_PD10 = TAB10.LBACK(ii).gyro; lb_gyr_PD10 = fillmissing(lb_gyr_PD10,'spline'); lb_gyr_PD10 = bwfilt(lb_gyr_PD10,2,128,6,'low');

    % st_acc = TAB10.STERN(ii).acc; st_acc = fillmissing(st_acc,'spline'); st_acc = bwfilt(st_acc,2,128,10,'low');
    % st_gyr = TAB10.STERN(ii).gyro; st_gyr = fillmissing(st_gyr,'spline'); st_gyr = bwfilt(st_gyr,2,128,6,'low');
    % 
    % fh_acc = TAB10.FHEAD(ii).acc; fh_acc = fillmissing(fh_acc,'spline'); fh_acc = bwfilt(fh_acc,2,128,10,'low');
    % fh_gyr = TAB10.FHEAD(ii).gyro; fh_gyr = fillmissing(fh_gyr,'spline'); fh_gyr = bwfilt(fh_gyr,2,128,6,'low');


    contacts = TAB10.CONTANCTS{ii};
    side = TAB10.SIDE{ii};

    [r,c] = (find(isnan(contacts))); r = unique(r); contacts(r,:) = []; %TAB.CONTANCTS{ii} = contacts;
    side(:,r) = []; %TAB.SIDE{ii} = side;


    for co = 1:length(contacts)

%if
    end 
end

 ST
for ii = 31:45
    lb_acc_ST10 = TAB10.LBACK(ii).acc; lb_acc_ST10 = fillmissing(lb_acc_ST10,'spline'); lb_acc_ST10 = bwfilt(lb_acc_ST10,2,128,10,'low');
    lb_gyr_ST10 = TAB10.LBACK(ii).gyro; lb_gyr_ST10 = fillmissing(lb_gyr_ST10,'spline'); lb_gyr_ST10 = bwfilt(lb_gyr_ST10,2,128,6,'low');

    % st_acc = TAB10.STERN(ii).acc; st_acc = fillmissing(st_acc,'spline'); st_acc = bwfilt(st_acc,2,128,10,'low');
    % st_gyr = TAB10.STERN(ii).gyro; st_gyr = fillmissing(st_gyr,'spline'); st_gyr = bwfilt(st_gyr,2,128,6,'low');
    % 
    % fh_acc = TAB10.FHEAD(ii).acc; fh_acc = fillmissing(fh_acc,'spline'); fh_acc = bwfilt(fh_acc,2,128,10,'low');
    % fh_gyr = TAB10.FHEAD(ii).gyro; fh_gyr = fillmissing(fh_gyr,'spline'); fh_gyr = bwfilt(fh_gyr,2,128,6,'low');


    contacts = TAB10.CONTANCTS{ii};
    side = TAB10.SIDE{ii};

    [r,c] = (find(isnan(contacts))); r = unique(r); contacts(r,:) = []; %TAB.CONTANCTS{ii} = contacts;
    side(:,r) = []; %TAB.SIDE{ii} = side;


    for co = 1:length(contacts)

%if
       
        end
end

% TBI
    %% For TBI (Traumatic Brain Injury)
    for ii = 46:60   
    lb_acc_TBI10 = TAB10.LBACK(ii).acc; lb_acc_TBI10 = fillmissing(lb_acc_TBI10,'spline'); lb_acc_TBI10 = bwfilt(lb_acc_TBI10,2,128,10,'low');
    lb_gyr_TBI10 = TAB10.LBACK(ii).gyro; lb_gyr_TBI10 = fillmissing(lb_gyr_TBI10,'spline'); lb_gyr_TBI10 = bwfilt(lb_gyr_TBI10,2,128,6,'low');

    % st_acc = TAB10.STERN(ii).acc; st_acc = fillmissing(st_acc,'spline'); st_acc = bwfilt(st_acc,2,128,10,'low');
    % st_gyr = TAB10.STERN(ii).gyro; st_gyr = fillmissing(st_gyr,'spline'); st_gyr = bwfilt(st_gyr,2,128,6,'low');
    % 
    % fh_acc = TAB10.FHEAD(ii).acc; fh_acc = fillmissing(fh_acc,'spline'); fh_acc = bwfilt(fh_acc,2,128,10,'low');
    % fh_gyr = TAB10.FHEAD(ii).gyro; fh_gyr = fillmissing(fh_gyr,'spline'); fh_gyr = bwfilt(fh_gyr,2,128,6,'low');


    contacts = TAB10.CONTANCTS{ii};
    side = TAB10.SIDE{ii};

    [r,c] = (find(isnan(contacts))); r = unique(r); contacts(r,:) = []; %TAB.CONTANCTS{ii} = contacts;
    side(:,r) = []; %TAB.SIDE{ii} = side;


    for co = 1:length(contacts)

%if
    end
    end



   
Estrae il segnale Z da ogni soggetto.
Riempie i dati mancanti e applica un filtro passa basso.
Rimuove il trend dal segnale.
Costruisce lo spazio delle fasi (embedding) con:
Dimensione m = 5
Ritardo τ = 10 (≈ 80 ms se fs = 128 Hz)
Trova per ogni punto il vicino più vicino non temporale.
Calcola la divergenza logaritmica nel tempo per 20 passi.
Usa una regressione lineare per stimare la pendenza = Lyapunov exponent.
    %% Calcolo dell’Esponente di Lyapunov per ogni soggetto in TABreduced
LyE_individuale = zeros(height(TABreduced), 1); % inizializzazione

for ii = 1:height(TABreduced)
    % Dati accelerometrici: scegli asse verticale (Z)
    acc_data = TABreduced.LBACK(ii).acc;
    
    if isempty(acc_data)
        LyE_individuale(ii) = NaN;
        continue;
    end

    acc_data = fillmissing(acc_data,'spline');
    acc_data = bwfilt(acc_data,2,128,10,'low'); % filtro passa basso

    signal = acc_data(:,3); % asse Z
    signal = detrend(signal); % rimuove trend lineare

    % Parametri per embedding
    m = 5; % dimensione dell'embedding
    tau = 10; % ritardo in punti (≈80ms con 128Hz)
    fs = 128; % frequenza di campionamento

    % Ricostruzione dello spazio delle fasi
    N = length(signal) - (m-1)*tau;
    X = zeros(N, m);
    for k = 1:m
        X(:,k) = signal(1+(k-1)*tau : N+(k-1)*tau);
    end

    % Calcolo delle distanze iniziali
    d0 = zeros(N,1);
    for k = 1:N
        % Trova punto vicino escludendo vicini temporali
        dist = vecnorm(X - X(k,:), 2, 2);
        dist(abs((1:N)' - k) < 30) = Inf; % escludi 30 campioni vicini
        [d0(k), j(k)] = min(dist);
    end
% Calcolo divergenza media nel tempo
    max_iter = 20; % numero di passi temporali
    divergence = zeros(max_iter,1);

    for n = 1:max_iter
        valid = (j+n <= N) & (1:N-n)';
        dist_n = vecnorm(X(valid + n,:) - X(j(valid) + n,:), 2, 2);
        divergence(n) = mean(log(dist_n ./ d0(valid)), 'omitnan');
    end

    % Regressione lineare per stimare la Lyapunov exponent
    t = (1:max_iter) / fs;
    p = polyfit(t, divergence(1:max_iter)', 1);
    LyE_individuale(ii) = p(1); % coefficiente angolare = esponente di Lyapunov
end

% Aggiungo alla tabella
TABreduced.LyE = LyE_individuale;

% Visualizzazione
figure;
boxplot(LyE_individuale, gruppi);
ylabel('Lyapunov Exponent');
title('LyE per Gruppo (Healthy vs Patologici)');

% ANOVA
[p_lye, tbl_lye, stats_lye] = anova1(LyE_individuale, gruppi);
[c_lye,gnames_lye] = multcompare(stats_lye);
csv matlab