%% ========================================================================= %% PLAYER LOAD EQUATION COMPARISON — FIGURE GENERATION SCRIPT %% ========================================================================= % % Callaway, A.J., Ellis, S.A., & Williams, J. % "Not All Player Load Equations Are Equivalent: Formula Choice Alters % Load Magnitude and Interpretation in Elite Football" % % Produces all four manuscript figures: % % Figure 1 — Winter residual analysis for Butterworth cutoff selection % (reads raw game CSV files directly) % % Figure 2 — Bland-Altman agreement plots: PlayerLoad (EQ3) vs AcelT (EQ1) % and PlayerLoad (EQ3) vs PL_RE (EQ5), Kalman condition % % Figure 3 — Within-match rank-order preservation across equation families, % Kalman condition (four panels: A-D) % % Figure 4 — Effect of signal processing condition on PlayerLoad (EQ3): % distribution comparison and raw vs Kalman scatter % % Also produces (not a manuscript figure): % Figure_PSD — Power spectral density of positional signals % % REQUIREMENTS: % - MATLAB R2020b or later % - Signal Processing Toolbox (for filtfilt, butter, pwelch) % - Step 2 output file: Step2_PlayerLoad_Summary.mat % - Raw game CSV files: game_XXXX_xy.csv (for Figure 1 only) % % USAGE: % 1. Set all paths in the CONFIGURATION section below % 2. Run the full script (or run sections individually using Ctrl+Enter) % % ========================================================================= clear; clc; close all; %% ========================================================================= %% SECTION 0: CONFIGURATION — UPDATE ALL PATHS BEFORE RUNNING %% ========================================================================= % --- Input files ---------------------------------------------------------- % Step 2 summary .mat file (required for Figures 2, 3, 4) CONFIG.step2_mat = '/path/to/Step2_PlayerLoad_Summary.mat'; % Folder containing raw game CSV files: game_XXXX_xy.csv (required for Figure 1) CONFIG.csv_folder = '/path/to/game_csv_files/'; % Optional: folder containing Step 1 .mat files (required for PSD figure only) % Leave empty to skip PSD figure: CONFIG.step1_folder = ''; CONFIG.step1_folder = '/path/to/step1_mat_files/'; % --- Output --------------------------------------------------------------- % Folder where all figures will be saved (.tiff and .svg) CONFIG.fig_folder = '/path/to/output/figures/'; % --- Signal parameters ---------------------------------------------------- CONFIG.fps = 29.9697; % nominal sampling rate (Hz) CONFIG.dt = 1 / CONFIG.fps; CONFIG.fig_dpi = 300; % --- Figure dimensions (cm) ----------------------------------------------- CONFIG.fig_w_single = 8.5; % single column CONFIG.fig_w_double = 17.5; % double column CONFIG.fig_h_single = 8.0; CONFIG.fig_h_double = 9.0; % --- Winter residual analysis parameters (Figure 1) ---------------------- CONFIG.fc_min = 0.5; % minimum candidate cutoff (Hz) CONFIG.fc_max = 10.0; % maximum candidate cutoff (Hz) CONFIG.fc_step = 0.25; % step size (Hz) CONFIG.butter_order = 2; % filter order (zero-phase via filtfilt = effective 4th order) CONFIG.n_games = 20; % number of game CSVs to sample from CONFIG.n_players = 100; % target number of player-period sequences CONFIG.min_frames = 150; % minimum sequence length (~5 s) CONFIG.max_frames = 900; % maximum sequence length (~30 s) CONFIG.rng_seed = 42; % --- Colour scheme (colourblind-friendly, consistent across all figures) -- COL.instant = [0.122, 0.471, 0.706]; % blue — instantaneous magnitude family COL.jerk = [0.200, 0.627, 0.173]; % green — jerk-based family COL.energy = [0.890, 0.102, 0.110]; % red — energy-weighted family COL.neutral = [0.500, 0.500, 0.500]; % grey COL.optimal = [0.600, 0.000, 0.600]; % purple — Winter optimal cutoff COL.hz2 = [0.122, 0.471, 0.706]; % blue — 2 Hz line COL.hz4 = [0.890, 0.102, 0.110]; % red — 4 Hz line COL.fit = [0.400, 0.400, 0.400]; % grey — fitted lines % --- Equation metadata ---------------------------------------------------- EQ_NAMES = {'AcelT','PL_{RT}','PlayerLoad','ImpulseLoad','PL_{RE}', ... 'TotalLoad','PL_x','PL_y','PL_z','PL2D','PLTM'}; EQ_MAT_NAMES = {'EQ1_ActiGraph','EQ2_PL_RT_Cum','EQ3_PlayerLoad', ... 'EQ4_ImpLoad','EQ5_PL_RE','EQ6_TotalLoad', ... 'EQ7_PL_x','EQ8_PL_y','EQ9_PL_z', ... 'EQ10_PL2D','EQ11_PLTM'}; EQ_FAMILY = [1, 2, 2, 1, 3, 2, 2, 2, 2, 2, 2]; % 1 = instantaneous magnitude, 2 = jerk-based, 3 = energy-weighted % --- Positional groups (used for colour coding in Figure 2) --------------- POS_GROUPS = {'GK','CB','FB','CM','W','ST'}; POS_MAP = containers.Map( ... {'GK','RCB','LCB','MCB','RB','LB','RWB','LWB','CM','DM','LW','RW','LM','RM','AM','CF'}, ... {'GK','CB', 'CB', 'CB', 'FB','FB','FB', 'FB', 'CM','CM','W', 'W', 'W', 'W', 'W', 'ST'}); POS_COLS = containers.Map( ... POS_GROUPS, ... {[0.89, 0.10, 0.11], ... % GK red [0.12, 0.47, 0.71], ... % CB blue [0.20, 0.63, 0.17], ... % FB green [1.00, 0.50, 0.00], ... % CM orange [0.60, 0.31, 0.64], ... % W purple [0.65, 0.34, 0.14]}); % ST brown % ========================================================================= fprintf('========== PLAYER LOAD FIGURE GENERATION ==========\n\n'); if ~exist(CONFIG.fig_folder, 'dir') mkdir(CONFIG.fig_folder); fprintf('Created output folder: %s\n\n', CONFIG.fig_folder); end %% ========================================================================= %% SECTION 1: FIGURE 1 — WINTER RESIDUAL ANALYSIS %% ========================================================================= % Implements Winter (2009) residual analysis to determine optimal Butterworth % cutoff frequency for raw XY positional data. Reads directly from game CSVs. % % Output: Figure1_WinterResidual.tiff / .svg % ========================================================================= fprintf('--- FIGURE 1: Winter Residual Analysis ---\n\n'); rng(CONFIG.rng_seed); fc_vals = CONFIG.fc_min : CONFIG.fc_step : CONFIG.fc_max; n_fc = numel(fc_vals); % Precompute Butterworth filter coefficients for each candidate cutoff butter_b = cell(n_fc, 1); butter_a = cell(n_fc, 1); for ifc = 1:n_fc Wn = min(fc_vals(ifc) / (CONFIG.fps / 2), 0.9999); [butter_b{ifc}, butter_a{ifc}] = butter(CONFIG.butter_order, Wn, 'low'); end % --- Collect position sequences from game CSV files --- fprintf(' Collecting x_raw sequences from CSV files...\n'); csv_files = dir(fullfile(CONFIG.csv_folder, 'game_*_xy.csv')); n_csv = numel(csv_files); if n_csv == 0 warning('No game CSV files found in %s — Figure 1 skipped.', CONFIG.csv_folder); else n_games_use = min(CONFIG.n_games, n_csv); game_idx = randperm(n_csv, n_games_use); sequences = {}; sequences_ks = {}; for ig = 1:n_games_use if numel(sequences) >= CONFIG.n_players; break; end csvFile = fullfile(CONFIG.csv_folder, csv_files(game_idx(ig)).name); try T = readtable(csvFile, 'TextType', 'string'); catch ME fprintf(' LOAD ERROR (%s): %s\n', csv_files(game_idx(ig)).name, ME.message); continue; end T = T(T.team ~= "ball", :); % Resolve duplicate frames grpKey = strcat(string(T.game_id),'_',string(T.frame),'_', ... string(T.team),'_',string(T.jersey)); [~, ia, ic] = unique(grpKey, 'stable'); grpCounts = accumarray(ic, 1); dupGroups = find(grpCounts > 1); for dg = 1:numel(dupGroups) rows = find(ic == dupGroups(dg)); T.x_raw(rows(1)) = mean(T.x_raw(rows), 'omitnan'); end T = T(ia, :); playerKeys = strcat(string(T.team),'_',string(T.jersey),'_',string(T.period)); uniqueKeys = uniqueKeys_shuffled(playerKeys); for ik = 1:numel(uniqueKeys) if numel(sequences) >= CONFIG.n_players; break; end seg = sortrows(T(playerKeys == uniqueKeys(ik), :), 'frame'); xraw = double(seg.x_raw); xsmth = double(seg.x_smooth); valid = isfinite(xraw) & isfinite(xsmth); xraw = xraw(valid); xsmth = xsmth(valid); if numel(xraw) < CONFIG.min_frames; continue; end if numel(xraw) > CONFIG.max_frames s = randi(numel(xraw) - CONFIG.max_frames); xraw = xraw(s : s + CONFIG.max_frames - 1); xsmth = xsmth(s : s + CONFIG.max_frames - 1); end sequences{end+1} = xraw; %#ok sequences_ks{end+1} = xsmth; %#ok end fprintf(' [%d/%d] %s — sequences so far: %d\n', ig, n_games_use, ... csv_files(game_idx(ig)).name, numel(sequences)); end n_used = numel(sequences); fprintf('\n Sequences collected: %d\n\n', n_used); if n_used == 0 warning('No valid sequences — Figure 1 skipped.'); else % --- Compute residuals --- fprintf(' Computing residuals...\n'); residuals = NaN(n_used, n_fc); for iseq = 1:n_used sig = sequences{iseq}; for ifc = 1:n_fc try sig_filt = filtfilt(butter_b{ifc}, butter_a{ifc}, sig); residuals(iseq, ifc) = sqrt(mean((sig - sig_filt).^2)); catch residuals(iseq, ifc) = NaN; end end if mod(iseq, 10) == 0 || iseq == n_used fprintf(' %d/%d sequences\n', iseq, n_used); end end res_mean = mean(residuals, 1, 'omitnan'); res_sd = std(residuals, 0, 1, 'omitnan'); res_sem = res_sd / sqrt(n_used); % --- Elbow detection: two-segment linear fit --- best_sse = Inf; best_split = round(n_fc/2); best_p1 = [0 0]; best_p2 = [0 0]; for isplit = 3:(n_fc-3) x1 = fc_vals(1:isplit)'; y1 = res_mean(1:isplit)'; x2 = fc_vals(isplit:end)'; y2 = res_mean(isplit:end)'; p1 = polyfit(x1, y1, 1); p2 = polyfit(x2, y2, 1); sse = sum((y1-polyval(p1,x1)).^2) + sum((y2-polyval(p2,x2)).^2); if sse < best_sse best_sse = sse; best_split = isplit; best_p1 = p1; best_p2 = p2; end end if abs(best_p1(1) - best_p2(1)) > 1e-10 fc_optimal = (best_p2(2)-best_p1(2)) / (best_p1(1)-best_p2(1)); fc_optimal = max(CONFIG.fc_min, min(CONFIG.fc_max, fc_optimal)); else fc_optimal = fc_vals(best_split); warning('Lines nearly parallel — using split point as optimal cutoff.'); end % --- Kalman equivalence --- kalman_res = NaN(n_used, 1); for iseq = 1:n_used xr = sequences{iseq}; xk = sequences_ks{iseq}; n = min(numel(xr), numel(xk)); if n >= CONFIG.min_frames kalman_res(iseq) = sqrt(mean((xr(1:n)-xk(1:n)).^2)); end end kalman_res_mean = mean(kalman_res, 'omitnan'); kalman_res_sd = std(kalman_res, 'omitnan'); fc_kalman_interp = interp1(res_mean(end:-1:1), fc_vals(end:-1:1), ... kalman_res_mean, 'linear', 'extrap'); fc_kalman_interp = max(CONFIG.fc_min, min(CONFIG.fc_max, fc_kalman_interp)); fprintf('\n Optimal cutoff: %.2f Hz\n', fc_optimal); fprintf(' Kalman RMS residual: %.4f m (equiv. ~%.2f Hz)\n\n', ... kalman_res_mean, fc_kalman_interp); % --- Plot Figure 1 --- fig1 = figure('Units','centimeters', ... 'Position',[2 2 CONFIG.fig_w_double CONFIG.fig_h_double*1.1], ... 'Color','w', 'PaperPositionMode','auto'); % Panel A: mean ± SD with elbow fit ax1a = subplot(1,2,1); hold on; y_max = max(res_mean + res_sd) * 1.15; fill([fc_vals, fliplr(fc_vals)], [res_mean+res_sd, fliplr(res_mean-res_sd)], ... COL.jerk, 'FaceAlpha', 0.15, 'EdgeColor', 'none', 'HandleVisibility','off'); plot(fc_vals, res_mean, '-', 'Color', COL.jerk, 'LineWidth', 2.0, ... 'DisplayName', sprintf('Mean \\pm SD (n = %d)', n_used)); x_left = linspace(fc_vals(1), fc_vals(best_split), 60); x_right = linspace(fc_vals(best_split), fc_vals(end), 60); plot(x_left, polyval(best_p1, x_left), '--', 'Color', COL.fit, ... 'LineWidth', 1.2, 'HandleVisibility','off'); plot(x_right, polyval(best_p2, x_right), '--', 'Color', COL.fit, ... 'LineWidth', 1.2, 'HandleVisibility','off'); plot([2 2], [0 y_max], '--', 'Color', COL.hz2, 'LineWidth', 1.5, 'HandleVisibility','off'); plot([4 4], [0 y_max], '--', 'Color', COL.hz4, 'LineWidth', 1.5, 'HandleVisibility','off'); plot([fc_optimal fc_optimal], [0 y_max], '-', 'Color', COL.optimal, 'LineWidth', 2.0, 'HandleVisibility','off'); plot([CONFIG.fc_min CONFIG.fc_max], [kalman_res_mean kalman_res_mean], ... ':', 'Color', [0.3 0.3 0.3], 'LineWidth', 1.4, 'HandleVisibility','off'); text(2.08, y_max*0.97, '2 Hz', 'Color', COL.hz2, 'FontSize', 8, 'FontWeight','bold', 'VerticalAlignment','top'); text(4.08, y_max*0.97, '4 Hz', 'Color', COL.hz4, 'FontSize', 8, 'FontWeight','bold', 'VerticalAlignment','top'); text(fc_optimal+0.15, y_max*0.80, sprintf('Optimal\n%.2f Hz', fc_optimal), ... 'Color', COL.optimal, 'FontSize', 8, 'FontWeight','bold', 'VerticalAlignment','top'); text(CONFIG.fc_max*0.98, kalman_res_mean*1.08, sprintf('Kalman \\approx %.2f Hz', fc_kalman_interp), ... 'Color', [0.2 0.2 0.2], 'FontSize', 8, 'FontWeight','bold', 'HorizontalAlignment','right'); xlabel('Cutoff frequency (Hz)', 'FontSize', 10); ylabel('Residual (RMS, m)', 'FontSize', 10); legend('Location','northwest', 'FontSize', 8); xlim([CONFIG.fc_min CONFIG.fc_max]); ylim([0 y_max]); box off; grid on; set(gca, 'FontSize', 9); text(-0.14, 1.05, 'A', 'Units','normalized', 'FontSize', 11, 'FontWeight','bold'); % Panel B: individual traces + mean ± SEM ax1b = subplot(1,2,2); hold on; y_max2 = max(res_mean) * 1.4; for iseq = 1:n_used plot(fc_vals, residuals(iseq,:), '-', 'Color', [COL.jerk 0.10], ... 'LineWidth', 0.6, 'HandleVisibility','off'); end fill([fc_vals, fliplr(fc_vals)], [res_mean+res_sem, fliplr(res_mean-res_sem)], ... COL.jerk, 'FaceAlpha', 0.40, 'EdgeColor','none', 'HandleVisibility','off'); plot(fc_vals, res_mean, '-', 'Color', COL.jerk, 'LineWidth', 2.0, ... 'DisplayName', 'Mean \\pm SEM'); plot([2 2], [0 y_max2], '--', 'Color', COL.hz2, 'LineWidth', 1.5, 'HandleVisibility','off'); plot([4 4], [0 y_max2], '--', 'Color', COL.hz4, 'LineWidth', 1.5, 'HandleVisibility','off'); plot([fc_optimal fc_optimal], [0 y_max2], '-', 'Color', COL.optimal, 'LineWidth', 2.0, 'HandleVisibility','off'); text(2.08, y_max2*0.97, '2 Hz', 'Color', COL.hz2, 'FontSize', 8, 'FontWeight','bold', 'VerticalAlignment','top'); text(4.08, y_max2*0.97, '4 Hz', 'Color', COL.hz4, 'FontSize', 8, 'FontWeight','bold', 'VerticalAlignment','top'); xlabel('Cutoff frequency (Hz)', 'FontSize', 10); ylabel('Residual (RMS, m)', 'FontSize', 10); legend('Location','northwest', 'FontSize', 8); xlim([CONFIG.fc_min CONFIG.fc_max]); ylim([0 y_max2]); box off; grid on; set(gca, 'FontSize', 9); text(-0.14, 1.05, 'B', 'Units','normalized', 'FontSize', 11, 'FontWeight','bold'); save_figure(fig1, fullfile(CONFIG.fig_folder, 'Figure1_WinterResidual'), CONFIG.fig_dpi); fprintf(' Saved Figure 1\n\n'); end end %% ========================================================================= %% SECTION 2: LOAD STEP 2 DATA (required for Figures 2, 3, 4) %% ========================================================================= fprintf('--- Loading Step 2 summary data ---\n'); S = load(CONFIG.step2_mat, 'T'); T = S.T; fprintf(' %d rows loaded\n\n', height(T)); % Pool periods: sum P1+P2 per player-match per condition fprintf(' Pooling periods (P1+P2)...\n'); grpVars = {'game_id','team','jersey','position_raw','smoothing_condition'}; Tp = groupsummary(T, grpVars, 'sum', EQ_MAT_NAMES); for i = 1:numel(EQ_MAT_NAMES) Tp.Properties.VariableNames{['sum_' EQ_MAT_NAMES{i}]} = EQ_MAT_NAMES{i}; end fprintf(' %d player-match rows after pooling\n\n', height(Tp)); % Kalman-smoothed subset (primary analysis condition) Tk = Tp(Tp.smoothing_condition == "kalman", :); fprintf(' Kalman rows: %d\n\n', height(Tk)); % Build positional group colour vectors for Figure 2 raw_pos = string(Tk.position_raw); grp_col = repmat([0.7 0.7 0.7], height(Tk), 1); grp_lbl = repmat("OTHER", height(Tk), 1); pos_keys_cell = keys(POS_MAP); for ik = 1:numel(pos_keys_cell) pk = pos_keys_cell{ik}; grp = POS_MAP(pk); msk = raw_pos == string(pk); if any(msk) grp_col(msk, :) = repmat(POS_COLS(grp), sum(msk), 1); grp_lbl(msk) = string(grp); end end %% ========================================================================= %% SECTION 3: FIGURE 2 — BLAND-ALTMAN AGREEMENT PLOTS %% ========================================================================= % Two-panel Bland-Altman, Kalman condition: % Panel A: PlayerLoad (EQ3) vs AcelT (EQ1) — jerk vs instantaneous magnitude % Panel B: PlayerLoad (EQ3) vs PL_RE (EQ5) — jerk vs energy-weighted % % Legend shows positional groups only (bias/LoA/regression lines hidden). % Output: Figure2_BlandAltman.tiff / .svg % ========================================================================= fprintf('--- FIGURE 2: Bland-Altman plots ---\n'); pairs_fig2 = { ... {'EQ3_PlayerLoad', 'EQ1_ActiGraph', 'PlayerLoad vs AcelT'}, ... {'EQ3_PlayerLoad', 'EQ5_PL_RE', 'PlayerLoad vs PL_{RE}'}}; fig2 = figure('Units','centimeters', ... 'Position',[2 2 CONFIG.fig_w_double*1.2 CONFIG.fig_h_double*1.5], ... 'Color','w', 'PaperPositionMode','auto'); panel_lab2 = {'A','B'}; for ip = 1:2 eq_x_name = pairs_fig2{ip}{1}; eq_y_name = pairs_fig2{ip}{2}; x_vals = Tk.(eq_x_name); y_vals = Tk.(eq_y_name); valid = isfinite(x_vals) & isfinite(y_vals); x_v = x_vals(valid); y_v = y_vals(valid); lbl_v = grp_lbl(valid); diff_v = x_v - y_v; mean_v = (x_v + y_v) / 2; bias = mean(diff_v); sd_d = std(diff_v); loa_lo = bias - 1.96*sd_d; loa_hi = bias + 1.96*sd_d; r2_val = corr(x_v, y_v)^2; p_prop = polyfit(mean_v, diff_v, 1); ax2 = subplot(1, 2, ip); hold on; % Scatter points coloured by positional group for ig = 1:numel(POS_GROUPS) msk = lbl_v == POS_GROUPS{ig}; if any(msk) scatter(mean_v(msk), diff_v(msk), 12, POS_COLS(POS_GROUPS{ig}), 'filled', ... 'MarkerFaceAlpha', 0.45, 'MarkerEdgeAlpha', 0, ... 'DisplayName', POS_GROUPS{ig}); end end xl = [min(mean_v)*0.98, max(mean_v)*1.02]; % Bias and LoA lines — hidden from legend plot(xl, [bias bias], '-', 'Color', [0.15 0.15 0.15], 'LineWidth', 1.8, 'HandleVisibility','off'); plot(xl, [loa_lo loa_lo], '--', 'Color', [0.4 0.4 0.4], 'LineWidth', 1.2, 'HandleVisibility','off'); plot(xl, [loa_hi loa_hi], '--', 'Color', [0.4 0.4 0.4], 'LineWidth', 1.2, 'HandleVisibility','off'); % Proportional bias regression line — hidden from legend x_fit = linspace(xl(1), xl(2), 200); plot(x_fit, polyval(p_prop, x_fit), '-.', 'Color', [0.8 0.3 0.0], 'LineWidth', 1.4, ... 'HandleVisibility','off'); % Value annotations y_range = loa_hi - loa_lo; x_ann = xl(1) + 0.02*(xl(2)-xl(1)); text(xl(2)*0.98, bias + 0.04*y_range, sprintf('Bias = %.1f AU', bias), 'HorizontalAlignment','right', 'FontSize', 8, 'FontWeight','bold'); text(xl(2)*0.98, loa_hi + 0.04*y_range, sprintf('+1.96 SD = %.1f AU', loa_hi), 'HorizontalAlignment','right', 'FontSize', 7.5, 'Color', [0.4 0.4 0.4]); text(xl(2)*0.98, loa_lo - 0.06*y_range, sprintf('\x22121.96 SD = %.1f AU', loa_lo), 'HorizontalAlignment','right', 'FontSize', 7.5, 'Color', [0.4 0.4 0.4]); text(x_ann, loa_hi - 0.04*y_range, sprintf('R^2 = %.3f', r2_val), 'FontSize', 8.5, 'Color', [0.2 0.2 0.2]); text(x_ann, loa_hi - 0.10*y_range, sprintf('N = %d', sum(valid)), 'FontSize', 8.5, 'Color', [0.2 0.2 0.2]); xlim(xl); % Suppress x10^n exponent to prevent overlap with xlabel ax2.XAxis.Exponent = 0; xtickformat('%,.0f'); xlabel('Mean of PlayerLoad and comparator (AU)', 'FontSize', 9); ylabel('PlayerLoad \minus comparator (AU)', 'FontSize', 10); legend('Location','best', 'FontSize', 7.5, 'NumColumns', 2); box off; grid on; set(gca, 'FontSize', 9); % Expand bottom inset so xlabel is not clipped on export set(ax2, 'LooseInset', max(get(ax2,'TightInset'), [0.15 0.18 0.05 0.05])); % Panel label top-left text(-0.10, 1.05, panel_lab2{ip}, 'Units','normalized', ... 'FontSize', 12, 'FontWeight','bold', 'VerticalAlignment','bottom'); end save_figure(fig2, fullfile(CONFIG.fig_folder, 'Figure2_BlandAltman'), CONFIG.fig_dpi); fprintf(' Saved Figure 2\n\n'); %% ========================================================================= %% SECTION 4: FIGURE 3 — RANK-ORDER PRESERVATION %% ========================================================================= % Four-panel scatter of within-match percentile ranks, Kalman condition: % Panel A: EQ3 vs EQ2 (within-family — jerk) % Panel B: EQ3 vs EQ6 (within-family — jerk) % Panel C: EQ3 vs EQ1 (cross-family — jerk vs instantaneous) % Panel D: EQ3 vs EQ5 (cross-family — jerk vs energy-weighted) % % Output: Figure3_RankOrder.tiff / .svg % ========================================================================= fprintf('--- FIGURE 3: Rank-order preservation ---\n'); % Compute within-match percentile ranks for each equation (Kalman) fprintf(' Computing within-match percentile ranks...\n'); game_ids = Tk.game_id; uniq_games = unique(game_ids); n_obs = height(Tk); rank_mat = NaN(n_obs, 11); for ig = 1:numel(uniq_games) gmsk = game_ids == uniq_games(ig); idx = find(gmsk); nPl = numel(idx); if nPl < 2; continue; end for ie = 1:11 vals = Tk.(EQ_MAT_NAMES{ie})(idx); if sum(isfinite(vals)) < 2; continue; end [~, rank_order] = sort(vals); pct_rank = NaN(nPl, 1); pct_rank(rank_order) = (1:nPl)' / nPl * 100; rank_mat(idx, ie) = pct_rank; end end pairs_fig3 = { ... 3, 2, 'A', 'PlayerLoad vs PL_{RT}'; % within-family 3, 6, 'B', 'PlayerLoad vs TotalLoad'; % within-family 3, 1, 'C', 'PlayerLoad vs AcelT'; % cross-family 3, 5, 'D', 'PlayerLoad vs PL_{RE}'}; % cross-family panel_colors_fig3 = {COL.jerk, COL.jerk, COL.instant, COL.energy}; fig3 = figure('Units','centimeters', ... 'Position',[2 2 CONFIG.fig_w_double*1.2 CONFIG.fig_h_double*1.5], ... 'Color','w', 'PaperPositionMode','auto'); for ip = 1:4 ei = pairs_fig3{ip, 1}; ej = pairs_fig3{ip, 2}; p_lab = pairs_fig3{ip, 3}; rx = rank_mat(:, ei); ry = rank_mat(:, ej); valid = isfinite(rx) & isfinite(ry); rx_v = rx(valid); ry_v = ry(valid); rho = corr(rx_v, ry_v, 'Type', 'Spearman'); rmse = sqrt(mean((rx_v - ry_v).^2)); top10_mask = rx_v >= 90; if sum(top10_mask) > 0 pct_low = 100 * sum(ry_v(top10_mask) < 50) / sum(top10_mask); else pct_low = NaN; end axf = subplot(2, 2, ip); hold on; scatter(rx_v, ry_v, 6, panel_colors_fig3{ip}, 'filled', ... 'MarkerFaceAlpha', 0.25, 'MarkerEdgeAlpha', 0); % Identity line plot([0 100], [0 100], 'k-', 'LineWidth', 1.5); % ±10 percentile shaded band x_fill = [0 100 100 0]; fill(x_fill, min(max([10 110 100 0], 0), 100), [0.85 0.85 0.85], 'FaceAlpha', 0.25, 'EdgeColor','none'); fill(x_fill, min(max([-10 90 100 10], 0), 100), [0.85 0.85 0.85], 'FaceAlpha', 0.25, 'EdgeColor','none'); xlim([0 100]); ylim([0 100]); xlabel('PlayerLoad within-match percentile rank', 'FontSize', 9); ylabel(sprintf('%s within-match percentile rank', EQ_NAMES{ej}), 'FontSize', 9); % Statistics annotation box ann_str = sprintf('\\rho = %.3f\nRMSE_{rank} = %.1f%%', rho, rmse); if ~isnan(pct_low) ann_str = sprintf('%s\nTop-10%% misclassified: %.0f%%', ann_str, pct_low); end text(3, 93, ann_str, 'FontSize', 8.5, 'VerticalAlignment','top', ... 'BackgroundColor', [0.97 0.97 0.97], 'EdgeColor', [0.8 0.8 0.8], 'Margin', 3); box off; grid on; set(gca, 'FontSize', 9); axis square; % Panel label top-left text(-0.10, 1.05, p_lab, 'Units','normalized', ... 'FontSize', 12, 'FontWeight','bold', 'VerticalAlignment','bottom'); end save_figure(fig3, fullfile(CONFIG.fig_folder, 'Figure3_RankOrder'), CONFIG.fig_dpi); fprintf(' Saved Figure 3\n\n'); %% ========================================================================= %% SECTION 5: FIGURE 4 — SMOOTHING CONDITION COMPARISON %% ========================================================================= % Two panels: % Left: box-whisker distribution of EQ3 across all four conditions % Right: scatter of EQ3 raw vs Kalman on log-log scale % % Output: Figure4_Smoothing.tiff / .svg % ========================================================================= fprintf('--- FIGURE 4: Smoothing condition comparison ---\n'); conditions = {'raw','kalman','butter_2hz','butter_4hz'}; cond_labels = {'Raw','Kalman','Butterworth 2 Hz','Butterworth 4 Hz'}; cond_cols = {[0.7 0.7 0.7], COL.jerk, COL.instant, COL.energy}; fig4 = figure('Units','centimeters', ... 'Position',[2 2 CONFIG.fig_w_double*1.2 CONFIG.fig_h_double*1.4], ... 'Color','w', 'PaperPositionMode','auto'); % Panel A: EQ3 distribution by smoothing condition subplot(1, 2, 1); hold on; for ic = 1:4 mask = Tp.smoothing_condition == conditions{ic}; bp = boxplot(Tp.EQ3_PlayerLoad(mask), 'Positions', ic, 'Widths', 0.55, ... 'Colors', cond_cols{ic}, 'Symbol', '.'); set(bp, 'LineWidth', 1.2); end xlim([0.3 4.7]); xticks(1:4); xticklabels(cond_labels); xtickangle(20); ylabel('PlayerLoad (AU)', 'FontSize', 10); box off; grid on; set(gca, 'FontSize', 9); text(-0.14, 1.05, 'A', 'Units','normalized', 'FontSize', 11, 'FontWeight','bold'); % Panel B: raw vs Kalman scatter, log-log subplot(1, 2, 2); hold on; t_raw = Tp(Tp.smoothing_condition == "raw", {'game_id','team','jersey','EQ3_PlayerLoad'}); t_kal = Tp(Tp.smoothing_condition == "kalman", {'game_id','team','jersey','EQ3_PlayerLoad'}); t_raw.Properties.VariableNames{'EQ3_PlayerLoad'} = 'EQ3_raw'; t_kal.Properties.VariableNames{'EQ3_PlayerLoad'} = 'EQ3_kal'; merged = innerjoin(t_raw, t_kal, 'Keys', {'game_id','team','jersey'}); valid_m = isfinite(merged.EQ3_raw) & isfinite(merged.EQ3_kal) & ... merged.EQ3_raw > 0 & merged.EQ3_kal > 0; scatter(merged.EQ3_kal(valid_m), merged.EQ3_raw(valid_m), 10, ... COL.jerk, 'filled', 'MarkerFaceAlpha', 0.3, 'MarkerEdgeAlpha', 0); lim_lo = min([merged.EQ3_kal(valid_m); merged.EQ3_raw(valid_m)]); lim_hi = max([merged.EQ3_kal(valid_m); merged.EQ3_raw(valid_m)]); plot([lim_lo lim_hi], [lim_lo lim_hi], 'k-', 'LineWidth', 1.4); ratio_mean = mean(merged.EQ3_raw(valid_m) ./ merged.EQ3_kal(valid_m)); ratio_sd = std(merged.EQ3_raw(valid_m) ./ merged.EQ3_kal(valid_m)); text(0.05, 0.92, sprintf('Mean ratio raw/Kalman = %.1f \\pm %.1f', ratio_mean, ratio_sd), ... 'Units','normalized', 'FontSize', 9, ... 'BackgroundColor', [0.97 0.97 0.97], 'EdgeColor', [0.8 0.8 0.8], 'Margin', 3); set(gca, 'XScale','log', 'YScale','log'); xlabel('PlayerLoad Kalman-smoothed (AU)', 'FontSize', 10); ylabel('PlayerLoad Raw (AU)', 'FontSize', 10); box off; grid on; set(gca, 'FontSize', 9); text(-0.14, 1.05, 'B', 'Units','normalized', 'FontSize', 11, 'FontWeight','bold'); save_figure(fig4, fullfile(CONFIG.fig_folder, 'Figure4_Smoothing'), CONFIG.fig_dpi); fprintf(' Saved Figure 4\n\n'); %% ========================================================================= %% SECTION 6: PSD FIGURE — NOT A MANUSCRIPT FIGURE %% ========================================================================= % Power spectral density of positional signals (raw vs Kalman). % Retained for supplementary or methods documentation use. % Requires Step 1 .mat files. Set CONFIG.step1_folder = '' to skip. % % Output: Figure_PSD_notInManuscript.tiff / .svg % ========================================================================= fprintf('--- PSD figure (not in manuscript) ---\n'); if isempty(CONFIG.step1_folder) fprintf(' Skipped (CONFIG.step1_folder is empty)\n\n'); else step1_files = dir(fullfile(CONFIG.step1_folder, 'step1_*.mat')); if isempty(step1_files) fprintf(' Skipped (no step1_*.mat files found)\n\n'); else S1 = load(fullfile(CONFIG.step1_folder, step1_files(1).name), 'data'); d1 = S1.data; mask_cm = string(d1.position_raw) == "CM" & string(d1.smoothing_condition) == "kalman"; if sum(mask_cm) < 300 mask_cm = string(d1.smoothing_condition) == "kalman"; end idx_cm = find(mask_cm, 3000, 'first'); ax_seg = double(d1.ax(idx_cm)); mask_rw = string(d1.smoothing_condition) == "raw"; idx_rw = find(mask_rw, numel(idx_cm), 'first'); ax_raw = double(d1.ax(idx_rw)); ax_seg = ax_seg(isfinite(ax_seg)); ax_raw = ax_raw(isfinite(ax_raw)); n_use = min(numel(ax_seg), numel(ax_raw)); ax_seg = ax_seg(1:n_use); ax_raw = ax_raw(1:n_use); nw = min(512, floor(n_use/4)); nfft = 2^nextpow2(nw*2); [pxx_k, f_k] = pwelch(ax_seg, hann(nw), round(nw/2), nfft, CONFIG.fps); [pxx_r, f_r] = pwelch(ax_raw, hann(nw), round(nw/2), nfft, CONFIG.fps); fig_psd = figure('Units','centimeters', ... 'Position',[2 2 CONFIG.fig_w_double CONFIG.fig_h_double*1.1], ... 'Color','w', 'PaperPositionMode','auto'); subplot(1,2,1); hold on; y_lim_l = [min(10*log10(pxx_r)), max(10*log10(pxx_r))]; plot(f_r, 10*log10(pxx_r), 'Color', [0.6 0.6 0.6], 'LineWidth', 1.2, 'DisplayName','Raw'); plot(f_k, 10*log10(pxx_k), 'Color', COL.jerk, 'LineWidth', 1.6, 'DisplayName','Kalman-smoothed'); plot([2 2], y_lim_l, '--', 'Color', COL.instant, 'LineWidth', 1.5, 'HandleVisibility','off'); plot([4 4], y_lim_l, '--', 'Color', COL.energy, 'LineWidth', 1.5, 'HandleVisibility','off'); text(2, y_lim_l(2)*0.97, ' 2 Hz', 'Color', COL.instant, 'FontSize', 8, 'FontWeight','bold', 'VerticalAlignment','top'); text(4, y_lim_l(2)*0.97, ' 4 Hz', 'Color', COL.energy, 'FontSize', 8, 'FontWeight','bold', 'VerticalAlignment','top'); xlabel('Frequency (Hz)', 'FontSize', 10); ylabel('Power spectral density (dB Hz^{-1})', 'FontSize', 10); legend('Location','northeast', 'FontSize', 8); xlim([0 min(15, CONFIG.fps/2)]); box off; grid on; set(gca, 'FontSize', 9); subplot(1,2,2); hold on; cum_k = cumsum(pxx_k)./sum(pxx_k)*100; cum_r = cumsum(pxx_r)./sum(pxx_r)*100; plot(f_r, cum_r, 'Color', [0.6 0.6 0.6], 'LineWidth', 1.2, 'DisplayName','Raw'); plot(f_k, cum_k, 'Color', COL.jerk, 'LineWidth', 1.6, 'DisplayName','Kalman-smoothed'); plot([2 2], [0 105], '--', 'Color', COL.instant, 'LineWidth', 1.5, 'HandleVisibility','off'); plot([4 4], [0 105], '--', 'Color', COL.energy, 'LineWidth', 1.5, 'HandleVisibility','off'); plot([0 min(15,CONFIG.fps/2)], [80 80], ':', 'Color', [0.3 0.3 0.3], 'LineWidth', 1.2, 'HandleVisibility','off'); pct2 = interp1(f_k, cum_k, 2); pct4 = interp1(f_k, cum_k, 4); text(2.15, pct2-4, sprintf('%.0f%% below 2 Hz', pct2), 'FontSize', 8, 'Color', COL.instant); text(4.15, pct4-4, sprintf('%.0f%% below 4 Hz', pct4), 'FontSize', 8, 'Color', COL.energy); xlabel('Frequency (Hz)', 'FontSize', 10); ylabel('Cumulative spectral power (%)', 'FontSize', 10); legend('Location','southeast', 'FontSize', 8); xlim([0 min(15, CONFIG.fps/2)]); ylim([0 105]); box off; grid on; set(gca, 'FontSize', 9); save_figure(fig_psd, fullfile(CONFIG.fig_folder, 'Figure_PSD_notInManuscript'), CONFIG.fig_dpi); fprintf(' Saved PSD figure\n\n'); end end %% ========================================================================= %% DONE %% ========================================================================= fprintf('========== ALL FIGURES COMPLETE ==========\n'); fprintf('Output folder: %s\n\n', CONFIG.fig_folder); figfiles = dir(fullfile(CONFIG.fig_folder, 'Figure*.tiff')); for i = 1:numel(figfiles) fprintf(' %s\n', figfiles(i).name); end %% ========================================================================= %% LOCAL FUNCTIONS %% ========================================================================= function save_figure(fig, filepath_noext, dpi) % Save figure as both high-resolution TIFF and SVG. set(fig, 'PaperPositionMode', 'auto'); try exportgraphics(fig, [filepath_noext '.tiff'], 'Resolution', dpi, 'ContentType', 'image'); catch print(fig, [filepath_noext '.tiff'], '-dtiff', sprintf('-r%d', dpi)); end try print(fig, [filepath_noext '.svg'], '-dsvg'); catch warning('SVG export failed for %s', filepath_noext); end end function ukeys = uniqueKeys_shuffled(playerKeys) % Return unique keys in shuffled order (helper for Figure 1 sequence sampling). ukeys = unique(playerKeys, 'stable'); ukeys = ukeys(randperm(numel(ukeys))); end