% pop_fda_6 - Perform ERP analysis using PACE: Principal Analysis by 
%             Conditional Expectation. 
%             This function supports both Time-Domain and Frequency-Domain 
%             analyses. It can optionally compare reconstructed signals 
%             against true ERP values.
%
% Usage:
%   pop_fda_6(EEG);   % Opens a GUI prompting the user to select analysis options
%
% Inputs:
%   EEG - EEGLAB dataset containing epoched data with ICA weights computed.
%
% Outputs:
%   In version 0.1, the function does not return any output except a command 
%   string appended to the EEGLAB history.
%
% Version History:
%   Version 0.1 (2025)
%       - Initial release of PACE-based ERP analysis.
%       - GUI interface for selecting PACE options.
%       - Integration of PACE v2.17 MATLAB implementation.
%       - Supports both time-domain and frequency-domain FPCA.
%       - Exports results by frequency bands in separate folders.
%       - Produces Excel files for selected model-selection criteria:
%             * Fraction of Variance Explained (FVE)
%             * Akaike Information Criterion (AIC)
%             * Bayesian Information Criterion (BIC)
%
% Author:
%   Mohammad Fayaz. FDA plug-in GUI code for EEGLAB.
%   E-mail: Mohammad.Fayaz.89@gmail.com
%
%   This function uses (with permission) implementations from:
%
%       PACE Version 2.17  
%       Websites: 
%           https://www.stat.ucdavis.edu/PACE/
%           https://github.com/functionaldata/PACE_matlab
%
%       Main References:
%         [1] Yao, F., Müller, H.G., Clifford, A.J., Dueker, S.R., Follett, J., 
%             Lin, Y., Buchholz, B., Vogel, J.S. (2003).  
%             Shrinkage estimation for functional principal component scores, 
%             with application to the population kinetics of plasma folate.  
%             Biometrics, 59, 676–685.
%
%         [2] Yao, F., Müller, H.G., Wang, J.L. (2005).
%             Functional data analysis for sparse longitudinal data. 
%             Journal of the American Statistical Association, 100(470), 577–590.
%
% Example:
%   % Invoke PACE-based analysis (GUI mode):
%   pop_fda_6(EEG);
%
% Copyright:
%   © 2025 Mohammad Fayaz
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.



function com = pop_fda_6(EEG, varargin)
%%%% Tools -> Regreesion -> FoF
fprintf('Functional Principal Component Analysis with PACE \n');
fprintf('PACE: Principal Analysis by Conditional Expectation Version 2.17 \n');
fprintf('PACE Algorithm is available online at Website  https://www.stat.ucdavis.edu/PACE/  \n');
fprintf('Main References:  \n');
fprintf('[1] Yao, F., Müller, H.G., Clifford, A.J., Dueker, S.R., Follett, J., Lin, Y., Buchholz, B., Vogel, J.S. (2003). Shrinkage estimation for functional principal component scores, with application to the population kinetics of plasma folate. Biometrics 59, 676-685 \n')
fprintf('[2] Yao F, Müller HG, Wang JL. Functional data analysis for sparse longitudinal data. Journal of the American statistical association. 2005 Jun 1;100(470):577-90 \n')


com = '';
if nargin < 1
    help pop_sourcereconstruction;
    return
end

if nargin < 2
    %%%% GUI 
    g = [1 0.5 0.5];
    geometry = { g g g g g g g g g g g g g g g g g g};
    uilist = { ...
      { 'Style', 'text', 'string', 'Channel Selection', 'fontweight', 'bold' } {} {} ...
      { 'Style', 'text', 'string', '   Channel number', 'FontAngle', 'italic'  } ...
      { 'Style', 'edit', 'string', '1' 'tag' 'SelectedChanIndex'} {} ...
      { 'Style', 'text', 'string', '   Time Domain (ms)', 'FontAngle', 'italic'  } ...
      { 'Style', 'edit', 'string', '-0.2:0.8' 'tag' 'TimeDomain'} {} ...
      { 'Style', 'text', 'string', 'Frequency Range (Hz)', 'fontweight', 'bold' } {} {} ...
      { 'Style', 'text', 'string', '   Arbitrary Range (Hz)', 'FontAngle', 'italic'  } ...
      { 'Style', 'checkbox', 'string', 'Yes', 'value', 0, 'tag', 'ArbitraryRange_Yes' } ...
      { 'Style', 'edit', 'string', '0.5:60' 'tag' 'ArbitraryRange_Value'} ...
      { 'Style', 'text', 'string', '   Delta Band (Hz)', 'FontAngle', 'italic'  } ...
      { 'Style', 'checkbox', 'string', 'Yes', 'value', 1, 'tag', 'DeltaBand_Yes' } ...
      { 'Style', 'edit', 'string', '0.5:4' 'tag' 'DeltaBand_Value'} ...
      { 'Style', 'text', 'string', '   Theta Band (Hz)', 'FontAngle', 'italic'  } ...
      { 'Style', 'checkbox', 'string', 'Yes', 'value', 1, 'tag', 'ThetaBand_Yes' } ...
      { 'Style', 'edit', 'string', '4:8' 'tag' 'ThetaBand_Value'} ...
      { 'Style', 'text', 'string', '   Alpha Band (Hz)', 'FontAngle', 'italic'  } ...
      { 'Style', 'checkbox', 'string', 'Yes', 'value', 1, 'tag', 'AlphaBand_Yes' } ...
      { 'Style', 'edit', 'string', '8:13' 'tag' 'AlphaBand_Value'} ...
      { 'Style', 'text', 'string', '   Beta Band (Hz)', 'FontAngle', 'italic'  } ...
      { 'Style', 'checkbox', 'string', 'Yes', 'value', 1, 'tag', 'BetaBand_Yes' } ...
      { 'Style', 'edit', 'string', '13:30' 'tag' 'BetaBand_Value'} ...
      { 'Style', 'text', 'string', '   Gamma Band (Hz)', 'FontAngle', 'italic'  } ...
      { 'Style', 'checkbox', 'string', 'Yes', 'value', 1, 'tag', 'GammaBand_Yes' } ...
      { 'Style', 'edit', 'string', '30:50' 'tag' 'GammaBand_Value'} ...
      { 'Style', 'text', 'string', 'Functional PCA', 'fontweight', 'bold' } {} {} ...
      { 'Style', 'text', 'string', '   Criteria', 'FontAngle', 'italic'  } {} {} ...
      { 'Style', 'text', 'string', '      FVE (%)', 'FontAngle', 'italic'  }  ...
      { 'Style', 'checkbox', 'string', 'Yes', 'value', 1, 'tag', 'FVE_Yes' } ...
      { 'Style', 'edit', 'string', '90' 'tag' 'FVE_Value'} ...
      { 'Style', 'text', 'string', '      AIC', 'FontAngle', 'italic'  }  ...
      { 'Style', 'checkbox', 'string', 'Yes', 'value', 0, 'tag', 'AIC_Yes' } {}...
      { 'Style', 'text', 'string', '      BIC', 'FontAngle', 'italic'  }  ...
      { 'Style', 'checkbox', 'string', 'Yes', 'value', 0, 'tag', 'BIC_Yes' } {}...
      { 'Style', 'text', 'string', '   Plots', 'FontAngle', 'italic'  } ...
      { 'Style', 'checkbox', 'string', 'Yes', 'value', 1, 'tag', 'Plots_Yes' } {}...
      { 'Style', 'text', 'string', 'Save the Results', 'fontweight', 'bold' } {} {} ...
      { 'Style', 'text', 'string', '   File Address (.xlsx)', 'FontAngle', 'italic'  } ...
      { 'Style', 'checkbox', 'string', 'Yes', 'value', 0, 'tag', 'Save_Yes' } ...
      { 'Style', 'edit', 'string', 'C:\' 'tag' 'Save_Value'} ...

};

       [ tmp1 tmp2 strhalt structout ] = inputgui( geometry, uilist, ...
           'pophelp(''pop_newtimef'');', 'EPOCH Analysis | FPCA with PACE for EEGLAB - Version 0.1');


    % Check if user canceled
    if isempty(structout)
        disp('User canceled the dialog. Exiting function.');
        return;  % stop the function safely
end  

%%%%%%%%% Bootstrap --> Residual Bootstrap or model based bootstrapp --->
%%%%%%%%% Improve it by Block Bootstrap

      %%% Import GUI items into the variables.
      %%%%% General Settings

            nChannel = EEG.nbchan;           % Number of Channenls
            nTrials = EEG.trials;            % Number of Trials
            ChanLocs = EEG.chanlocs;         % Channels Location
            Events = EEG.event;              % Events
            Epochs = EEG.epoch;              % Epochs
            Times = EEG.times;               % EEG Times
            EData = EEG.data;                % EEG Dataset
            ICAWeights = EEG.icawinv;        % ICA weights
            ICAData = EEG.icaact;            % ICA Source activation

            %%% Input by User (General)
            %%% Channel (Electrode) Selection
            %%% Channel  
            SelectedChanIndex_1 = str2num(structout.SelectedChanIndex);
            SelectedChan_1 = ChanLocs(:,SelectedChanIndex_1) ;
            NSelectedChan_1 = length(SelectedChanIndex_1);

            %%% Time Interval Selection 
            %%%% Channel 1
            TIS_1 = str2double(strsplit(structout.TimeDomain, ':')) ;
            TimeIntervalSele_1 = TIS_1; %%% Time Interval Selection (defualt: -200, 800 ms)
            TimeLO_1 = min(TimeIntervalSele_1) * 1.0e+03 ;
            TimeUP_1 = max(TimeIntervalSele_1) * 1.0e+03 ;
            TimeIntervalSeleInput_1      = Times(Times(1,:) >= TimeLO_1 & Times(1,:) <= TimeUP_1);
            TimeIntervalSeleIndexInput_1 = Times(1,:) >= TimeLO_1 & Times(1,:) <= TimeUP_1;
            
            %%% Analysis Data Set with selected time intervals
            ADS_EData_1 = EData(SelectedChanIndex_1,TimeIntervalSeleIndexInput_1,:);

            %%% Calcultae Power and Frequency 
            Fs = EEG.srate;
            % --- Inputs ---
            % ADS_EData_1: 1 × nTimePoints × nTrials
            % Fs: sampling frequency in Hz
            
            [~, nTimePoints, nTrials] = size(ADS_EData_1);
            
            % FFT parameters
            nFFT = 2^nextpow2(nTimePoints);       % Zero-padding
            f = Fs/2*linspace(0,1,nFFT/2+1)';    % Column vector of frequencies
            
            % Initialize dataset with first column as frequency
            data_with_freq = zeros(length(f), nTrials+1);
            data_with_freq(:,1) = f;             % First column = frequencies
            
            % Compute power spectra for each trial
            for tr = 1:nTrials
                signal = squeeze(ADS_EData_1(1,:,tr)); % Extract trial data
                Y = fft(signal, nFFT);
                P2 = abs(Y/nTimePoints).^2;           % Power spectrum
                P1 = P2(1:nFFT/2+1);
                P1(2:end-1) = 2*P1(2:end-1);          % Single-sided
                data_with_freq(:, tr+1) = P1;         % Save trial power
            end

         
            %%%% Frequency Settings
            F_AR_Yes   = structout.ArbitraryRange_Yes;
            F_AR_Val   = structout.ArbitraryRange_Value;
            F_DE_Yes   = structout.DeltaBand_Yes;
            F_DE_Val   = structout.DeltaBand_Value;
            F_TH_Yes   = structout.ThetaBand_Yes;
            F_TH_Val   = structout.ThetaBand_Value;
            F_AL_Yes   = structout.AlphaBand_Yes;
            F_AL_Val   = structout.AlphaBand_Value;
            F_BE_Yes   = structout.BetaBand_Yes;
            F_BE_Val   = structout.BetaBand_Value;
            F_GA_Yes   = structout.GammaBand_Yes;
            F_GA_Val   = structout.GammaBand_Value;
            F_FVE_Yes  = structout.FVE_Yes;
            F_FVE_Val  = structout.FVE_Value;
            F_AIC_Yes  = structout.AIC_Yes;
            F_BIC_Yes  = structout.BIC_Yes;
            F_PLT_Yes  = structout.Plots_Yes;
            Save_Yes   = structout.Save_Yes;
            Save_Value = structout.Save_Value;

            %%% Analysis 
            %%%% Make the Structure of Data

            % === Build Bands Structure ===
            Bands = struct();
            
            % --- Arbitrary Range ---
            if F_AR_Yes
                F_AR_Val = str2num(F_AR_Val);
                Bands.Arbitrary.Yes = true;
                Bands.Arbitrary.Name = 'Arbitrary';
                Bands.Arbitrary.LO = min(F_AR_Val);
                Bands.Arbitrary.UP = max(F_AR_Val);
                Bands.Arbitrary.Domain = Bands.Arbitrary.LO:Bands.Arbitrary.UP;
            else
                Bands.Arbitrary.Yes = false;
            end
            
            % --- Delta Band ---
            if F_DE_Yes
                F_DE_Val = str2num(F_DE_Val);
                Bands.Delta.Yes = true;
                Bands.Delta.Name = 'Delta';
                Bands.Delta.LO = min(F_DE_Val);
                Bands.Delta.UP = max(F_DE_Val);
                Bands.Delta.Domain = Bands.Delta.LO:Bands.Delta.UP;
            else
                Bands.Delta.Yes = false;
            end
            
            % --- Theta Band ---
            if F_TH_Yes
                F_TH_Val = str2num(F_TH_Val);
                Bands.Theta.Yes = true;
                Bands.Theta.Name = 'Theta';
                Bands.Theta.LO = min(F_TH_Val);
                Bands.Theta.UP = max(F_TH_Val);
                Bands.Theta.Domain = Bands.Theta.LO:Bands.Theta.UP;
            else
                Bands.Theta.Yes = false;
            end
            
            % --- Alpha Band ---
            if F_AL_Yes
                F_AL_Val = str2num(F_AL_Val);
                Bands.Alpha.Yes = true;
                Bands.Alpha.Name = 'Alpha';
                Bands.Alpha.LO = min(F_AL_Val);
                Bands.Alpha.UP = max(F_AL_Val);
                Bands.Alpha.Domain = Bands.Alpha.LO:Bands.Alpha.UP;
            else
                Bands.Alpha.Yes = false;
            end
            
            % --- Beta Band ---
            if F_BE_Yes
                F_BE_Val = str2num(F_BE_Val);
                Bands.Beta.Yes = true;
                Bands.Beta.Name = 'Beta';
                Bands.Beta.LO = min(F_BE_Val);
                Bands.Beta.UP = max(F_BE_Val);
                Bands.Beta.Domain = Bands.Beta.LO:Bands.Beta.UP;
            else
                Bands.Beta.Yes = false;
            end
            
            % --- Gamma Band ---
            if F_GA_Yes
                F_GA_Val = str2num(F_GA_Val);
                Bands.Gamma.Yes = true;
                Bands.Gamma.Name = 'Gamma';
                Bands.Gamma.LO = min(F_GA_Val);
                Bands.Gamma.UP = max(F_GA_Val);
                Bands.Gamma.Domain = Bands.Gamma.LO:Bands.Gamma.UP;
            else
                Bands.Gamma.Yes = false;
            end


            %%% FPCA based on PACE

            % --- Get list of selected bands ---
            bandFields = fieldnames(Bands);
            selectedBands = {};
            for i = 1:length(bandFields)
                if Bands.(bandFields{i}).Yes
                    selectedBands{end+1} = bandFields{i}; %#ok<SAGROW>
                end
            end          
            
            

            % --- Loop over each selected band for FPCA ---
            for b = 1:length(selectedBands)
                band = selectedBands{b};
                
                % Frequency range for this band
                fRange = [Bands.(band).LO, Bands.(band).UP];
            
                % --- Select indices for this band based on first column of data_with_freq ---
                freq_col = data_with_freq(:,1);           % frequency column
                idx_band = freq_col >= fRange(1) & freq_col <= fRange(2);
            
                % --- Extract band-limited frequencies and power ---
                f_band = freq_col(idx_band);              % selected frequencies
                power_band = data_with_freq(idx_band, 2:end);  % power for all trials
            
                % --- Convert to cell arrays for FPCA ---
                nTrials = size(power_band, 2);
                y_cell = cell(1, nTrials);
                t_cell = cell(1, nTrials);
                for i = 1:nTrials
                    y_cell{i} = power_band(:, i)';      % row vector of power
                    t_cell{i} = f_band(:)';             % row vector of frequency
                end 

                % --- Collect selected methods ---
                methods = {};
                if F_FVE_Yes
                    methods{end+1} = 'FVE';
                    FVE_threshold_val = str2double(F_FVE_Val)/100;  % user-provided threshold
                end
                if F_AIC_Yes
                    methods{end+1} = 'AIC';
                end
                if F_BIC_Yes
                    methods{end+1} = 'BIC';
                end
                
                % --- Create a single figure for this band ---
                % figure('Name', [band ' Band FPCA Results'], 'NumberTitle', 'off');
                
                nMethods = length(methods);
                for m = 1:nMethods
                    sel_method = methods{m};
                
                    % --- Set FPCA options ---
                    if strcmp(sel_method, 'FVE')
                        p = setOptions( ...
                            'yname', ['Power-' band '-' sel_method], ...
                            'regular', 2, ...
                            'selection_k', sel_method, ...
                            'FVE_threshold', FVE_threshold_val, ...
                            'screePlot', 0, ...
                            'designPlot', 0, ...
                            'corrPlot', 0, ...
                            'verbose', 'on' ...
                        );
                    else
                        % AIC or BIC
                        p = setOptions( ...
                            'yname', ['Power-' band '-' sel_method], ...
                            'regular', 2, ...
                            'selection_k', sel_method, ...
                            'screePlot', 0, ...
                            'designPlot', 0, ...
                            'corrPlot', 0, ...
                            'verbose', 'on' ...
                        );
                    end
                
                    % --- Run FPCA ---
                    yy    = FPCA(y_cell, t_cell, p);

               %     param = setOptions_bfda('smethod', 'bhm', ...
               %                             'cgrid', 1, ...
               %                              'mat', 1, ...
               %                              'M', 10000, ...
               %                              'Burnin', 2000, ...
               %                              'w', 1, ...
               %                              'ws', 1);


                %    [out_cgrid, param] = BFDA(y_cell, t_cell, param);


                    
                 %   param_uc = setOptions_bfda('smethod', 'bhm', ...
                 %                              'cgrid', 0, ...
                 %                               'mat', 1, ...
                 %                             'M',10000, ...
                 %                              'Burnin', 2000, ...
                 %                              'pace', 1, ...
                 %                             'ws', 0.1);

                 % [out_ucgrid, param_uc] = BFDA(y_cell, t_cell, param_uc);



                   
                    %extract some of the results for the plots below:
                    out1 = getVal(yy,'out1');      %vector of time points for mu, phi and ypred
                    mu = getVal(yy,'mu');          %estimated mean function
                    out21 = getVal(yy,'out21');    %vector of time points for xcov
                    xcovfit = getVal(yy,'xcovfit');%fitted covariance surface
                    xcov = getVal(yy,'xcov');      %estimated smooth covariance evaluated at out21
                    phi = getVal(yy,'phi');        %estimated eigenfunctions
                    no_opt = getVal(yy,'no_opt');  %number of optimal FPC scores
                    xcorr = getVal(yy,'xcorr');    %fitted correlation surface
                    xi = getVal(yy,'xi_est'); % estimated FPC scores
                    %========================================================================================================
                    ypred = FPCAeval(yy,[],out1); 

                    % --- Plots
                    if F_PLT_Yes 


        % --- Integrated figure with subplots ---
        figure('Name', [band ' Band FPCA Results - ' sel_method], 'NumberTitle', 'off', 'Position', [100 50 1200 800]);

        % --- Panel 1: Scree Plot + FVE ---
        ax1 = subplot(3,2,1);
        My_createScreePlot(yy, ax1);

        % --- Panel 2: Mean Function ---
        ax2 = subplot(3,2,2);
        plot(cell2mat(t_cell),cell2mat(y_cell),'k.'); hold on;
        plot(out1, mu, 'r', 'LineWidth', 2);
        xlabel('Frequency (Hz)'); ylabel('Power');
        title('Mean Function');
        legend('Data','Fitted','Location','Best'); grid on;

        % --- Panel 3: First k Eigenfunctions ---
        ax3 = subplot(3,2,3);
        % Get individual FVE
        FVE_cum = getVal(yy,'FVE');           % cumulative FVE
        FVE_ind = [FVE_cum(1), diff(FVE_cum)]; % individual FVE for each PC

        numEF = min(3, size(phi,2));
        hold on;
        for k1 = 1:numEF
            plot(out1, phi(:,k1), 'LineWidth', 2);
        end
        xlabel('Frequency (Hz)'); ylabel('Eigenfunction Value');
        title(sprintf('First %d Eigenfunctions', numEF));
        legend(arrayfun(@(x) sprintf('PC%d (FVE=%.1f%%)', x, FVE_ind(x)*100), 1:numEF, 'UniformOutput', false));
        grid on;

        % --- Panel 4: K-mode variation for 1st PC ---
        ax4 = subplot(3,2,4);
        My_KModeVariationPlot(yy, 1, ax4);

        % --- Panel 5: K-mode variation for 2nd PC (if exists) ---
        if no_opt >= 2
            ax5 = subplot(3,2,5);
            My_KModeVariationPlot(yy, 2, ax5);
        end

        % --- Panel 6: Reconstructed spectra for first PCs ---
        ax6 = subplot(3,2,6);
        hold on;
        for i = 1:nTrials
            y_recon = mu(:) + phi(:,1:no_opt)*xi(i,1:no_opt)';
            plot(out1, y_recon, 'LineWidth', 1.5);
        end
        xlabel('Frequency (Hz)'); ylabel('Reconstructed Power');
        title(sprintf('Reconstructed Spectra (first %d PCs)', no_opt));
        grid on;

                    end

        if Save_Yes 
           % --- Making Output Excel File
           % ==== Example inputs ====
            outputDir = Save_Value ;   % <-- Change this to your base path
            band = selectedBands{b};              % e.g. 'Alpha'
            methodName = methods{m};              % e.g. 'FPCA'
            
            % ==== Step 1: Create external folder ====
            bandDir = fullfile(outputDir, band);  % Folder per band
            if ~exist(bandDir, 'dir')
                mkdir(bandDir);
            end
            
            % ==== Step 2: Build filename ====
            filename = fullfile(bandDir, sprintf('%s_%s.xlsx', band, methodName));

            % ==== Step 3: ypred (trials × variable #observations) ====
            ypred_mat = cell2mat(ypred);   % Convert cell array into numeric matrix
            nObs = size(ypred_mat, 2);     % Number of observations per trial
            TrialID = (1:size(ypred_mat,1))';
            
            % Create dynamic column names
            obsNames = arrayfun(@(k) sprintf('Predicted%d', k), 1:nObs, 'UniformOutput', false);
            
            % Build table
            ypred_table = array2table([TrialID ypred_mat], ...
                'VariableNames', ['Trial', obsNames]);
            writetable(ypred_table, filename, 'Sheet', 'YPred');
            
            % ==== Step 4: xi (trials × PCs) ====
            nPCs_xi = size(xi, 2);
            xiNames = arrayfun(@(k) sprintf('PC%d_Score', k), 1:nPCs_xi, 'UniformOutput', false);
            
            xi_table = array2table(xi, 'VariableNames', xiNames);
            xi_table.Trial = (1:size(xi,1))'; 
            xi_table = movevars(xi_table, 'Trial', 'before', 1);
            writetable(xi_table, filename, 'Sheet', 'PCScores');
            
            % ==== Step 5: phi (Eigenfunctions) ====
            numPCs = size(phi,2);
            phiNames = arrayfun(@(k) sprintf('PC%d', k), 1:numPCs, 'UniformOutput', false);
            
            phi_table = array2table(phi, 'VariableNames', phiNames);
            phi_table.Frequency = (1:size(phi,1))'; 
            phi_table = movevars(phi_table, 'Frequency', 'before', 1);
            writetable(phi_table, filename, 'Sheet', 'Eigenfunctions');

            % ==== Final message ====
            fprintf('✅ Results successfully exported to Excel file:\n   %s\n', filename);

          end

         end
            end
end
        