function [ Data_matrix, tf ] = data_load( cfg, ListStim, varargin )
% [ Data_matrix, tf ] = data_load(cfg, ListStim) 
% Load stimuli data and generates appropriate matrix Data_matrix for the
% calculation of classification image
% Data_matrix is a N_TrialsLoad X DimCI matrix (Data_matrix(i,f,t) or
% Data_matrix(i,t) or Data_matrix(i,f), where i is the stimulus
% presentation number as specified in ordre_aleatoire ; tf is a structure
% array containing the time and/or frequency index
%
% varargin : 
% matrix_crea(cfg, ListStim, 'dimonly', 'yes') does not calculate
% Data_matrix
% matrix_crea(cfg, ListStim, 'WithSNR', 'yes') takes SNR into account
% matrix_crea(cfg, ListStim, 'zscore', 'no') does not perform a z-scoring
% on the resulting Data_matrix
% matrix_crea(cfg, ListStim, 'ROIs', ROI) returns the mean power in the
% regions defined in ROI
% matrix_crea(cfg, ListStim, 'dB', 'yes') returns the result in dB

dimonly = 0;
WithSNR = 0;
WithSignal = 0;
numerase = 0;
zscore = 1;
ROI = [];
dB = 0;

while ~isempty(varargin)
    if strcmp(varargin{1}, 'dimonly')
        dimonly = isyes(varargin{2});
    elseif strcmp(varargin{1}, 'WithSNR')
        WithSNR = isyes(varargin{2});
    elseif strcmp(varargin{1}, 'WithSignal')
        WithSignal = isyes(varargin{2});
    elseif strcmp(varargin{1}, 'zscore')
        zscore = isyes(varargin{2});
    elseif strcmp(varargin{1}, 'ROIs')
        ROI = varargin{2};
    elseif strcmp(varargin{1}, 'dB')
        dB = isyes(varargin{2});
    else
        error(['error : unrecognised argument : ' varargin{1}]);
    end
    varargin=varargin(3:end);
end

% Extract time and/or frequency index
WavFile = strcat (cfg.FolderInBruit, '/', ListStim(1).name);
[bruit,fs,~] = wavread(WavFile);
bruit=mean(bruit,2);
if strcmp(cfg.DimCI,'tf')
    [~,f,t,p_bruit] = spectrogram(bruit,cfg.Nwindow,cfg.Nwindow*cfg.overlap,cfg.NFFT,fs);
    cfg.freq_analysis_index = find(f>=cfg.freq_analysis(1) & f<=cfg.freq_analysis(2));
    f = f(cfg.freq_analysis_index);
    cfg.time_analysis_index = find(t>=cfg.time_analysis(1) & t<=cfg.time_analysis(2));
    t = t(cfg.time_analysis_index);
    p_bruit = p_bruit(cfg.freq_analysis_index,cfg.time_analysis_index);
    size_CI=size(p_bruit);
    if isempty(ROI)
        Data_matrix = zeros(cfg.N_TrialsLoad, size_CI(1), size_CI(2));
    else
        Data_matrix = zeros(cfg.N_TrialsLoad, size(ROI,3));
    end
    tf.f=f;
    tf.t=t;
elseif strcmp(cfg.DimCI,'lyon')
    Tcoch = 1/(fs/cfg.decimation);
    t=(Tcoch:Tcoch:(length(bruit)/fs))-Tcoch/2;
    [~, f] = DesignLyonFilters(fs,cfg.earQ,cfg.stepfactor);
    f=f(end:-1:1);
    cfg.freq_analysis_index = find(f>=cfg.freq_analysis(1) & f<=cfg.freq_analysis(2));
    f = f(cfg.freq_analysis_index);
    cfg.time_analysis_index = find(t>=cfg.time_analysis(1) & t<=cfg.time_analysis(2));
    t = t(cfg.time_analysis_index);
    size_CI=[length(f),length(t)];
    if isempty(ROI)
        Data_matrix = zeros(cfg.N_TrialsLoad, size_CI(1), size_CI(2));
    else
        Data_matrix = zeros(cfg.N_TrialsLoad, size(ROI,3));
    end
    tf.f=f;
    tf.t=t;
end

% Creating data matrix
if dimonly == 0
    fprintf('\n%% Creating data matrix %%\n');
    if strcmp(cfg.DimCI,'tf')
        for i=1:min(cfg.N_TrialsLoad,cfg.N)
            idx_stim = cfg.ordre_aleatoire(cfg.IdxTrialsLoad(i));
            fprintf(repmat('\b',1,numerase));
            msg=['Loading stim no ' num2str(i) ' of ' num2str(cfg.N_TrialsLoad)];
            if ~isempty(ROI)
                msg = [msg ' (ROIs)'];
            end
            fprintf(msg);
            numerase=numel(msg);
            
            WavFile = strcat(cd, '/', cfg.FolderInBruit, '/', ListStim(idx_stim).name);
            [bruit,fs,~] = wavread(WavFile);
            bruit=mean(bruit,2);
            
            if WithSNR
                Pbruit = mean(bruit.^2);
                A=sqrt((Pbruit)*10^(ListStim(idx_stim).RSB/10));
                bruit = bruit/A;
            end
            
            [~,f,t,p_bruit] = spectrogram(bruit,cfg.Nwindow,cfg.Nwindow*cfg.overlap,cfg.NFFT,fs);
            cfg.freq_analysis_index = find(f>=cfg.freq_analysis(1) & f<=cfg.freq_analysis(2));
            cfg.time_analysis_index = find(t>=cfg.time_analysis(1) & t<=cfg.time_analysis(2));
            p_bruit=p_bruit(cfg.freq_analysis_index,cfg.time_analysis_index);
            
            if ~isempty(ROI)
                for i_ROI=1:size(ROI,3)
                    p_ROI(i_ROI)=sum(sum(p_bruit.*ROI(:,:,i_ROI)));
                end
                Data_matrix(i, :, : ) = p_ROI;
            else
                Data_matrix(i, :, : ) = p_bruit;
            end
        end
    elseif strcmp(cfg.DimCI,'lyon')
        for i=1:min(cfg.N_TrialsLoad,cfg.N)
            idx_stim = cfg.ordre_aleatoire(cfg.IdxTrialsLoad(i));
            fprintf(repmat('\b',1,numerase));
            msg=['Loading stim no ' num2str(i) ' of ' num2str(cfg.N_TrialsLoad)];
            if ~isempty(ROI)
                msg = [msg ' (ROIs)'];
            end
            fprintf(msg);
            numerase=numel(msg);
            WavFile = strcat(cd, '/', cfg.FolderInBruit, '/', ListStim(idx_stim).name);
            [bruit,~,~] = wavread(WavFile);
            bruit=mean(bruit,2);
            if WithSignal
                WavFile = strcat(cd, '/', cfg.FolderInSignal, '/', ListStim(idx_stim).signal, '.wav');
                [signal,~,~] = wavread(WavFile);
                signal=mean(signal,2);
                [ stimulus,~] = Addition_RSB(signal, bruit, ListStim(idx_stim).RSB);
                bruit=stimulus;
            end
            if WithSNR
                Pbruit = mean(bruit.^2);
                A=sqrt((Pbruit)*10^(ListStim(idx_stim).RSB/10));
                bruit = bruit/A;
            end
            bruit_coch = LyonPassiveEar(bruit,fs,cfg.decimation,cfg.earQ,cfg.stepfactor,1,1);
            bruit_coch=bruit_coch(end:-1:1,:);
            bruit_coch=bruit_coch(cfg.freq_analysis_index,cfg.time_analysis_index);
            
            if ~isempty(ROI)
                for i_ROI=1:size(ROI,3)
                    p_ROI(i_ROI)=sum(sum(bruit_coch.*ROI(:,:,i_ROI)));
                end
                Data_matrix(i, :, : ) = p_ROI;
            else
                Data_matrix(i, :, : ) = bruit_coch;
            end
        end
    end
    fprintf('\n');
    
    if dB
        Data_matrix=10*log10(Data_matrix);
    end
    
    if zscore
        Data_matrix=(Data_matrix-mean(Data_matrix(:)))/std(Data_matrix(:));
    end
else
    Data_matrix =[];
end
    
end

