text:
% PURPOSE: Converts decoded "bin descriptor(s)" into a numeric structure, in order to be used by binlister.
%
% FORMAT:
%
% [BIN, isparsednumerically] = bdf2struct(BIN)
%
% Inputs:
%
%   BIN                    - BIN structure containing "bin descriptor(s)"
%
% Output
%
% BIN                      - BIN structure
% isparsednumerically      - boolean, 1=everything ok; 0= there is error.
%
%
% See binlister.m, decodebdf.m
%
%
% *** This function is part of ERPLAB Toolbox ***
% Author: Javier Lopez-Calderon
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2009

%b8d3721ed219e65100184c6b95db209bb8d3721ed219e65100184c6b95db209b
%
% ERPLAB Toolbox
% Copyright ? 2007 The Regents of the University of California
% Created by Javier Lopez-Calderon and Steven Luck
% Center for Mind and Brain, University of California, Davis,
% javlopez@ucdavis.edu, sjluck@ucdavis.edu
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

function [BIN, isparsednumerically] = bdf2struct(BIN, allowMixedSigns)

% allowMixedSigns ... allows the event codes within one '{}' to have
%                     different signs

if nargin < 2 || isempty(allowMixedSigns)
    allowMixedSigns = false;
end

if nargin < 1
        help bdf2struct
        return
end

%
% For debugging...
%
condition  = 0;  % 1 = displays processing comments;  0 = don't do that!

%
% Since the data extraction routine is the same for prehome, athome, and
% posthome, we will use the variable logEventSelector for looping.
%
isparsednumerically = 1;    % Numeric parsing is approved by default
logEventSelector{1} = 'prehome';  %(log event selector = LES)
logEventSelector{2} = 'athome';
logEventSelector{3} = 'posthome';
nBin = length(BIN);
fprintf('\n');

for iBin = 1:nBin       % Bin's loop

        isworkingFine = 1;   %boolean condition if everything goes well
        jLogES = 1;

        while  (jLogES <= length(logEventSelector)) && (isworkingFine) % LES's loop

        % Take each field and divide it by curly braces = sequencers
                [segmentedBrackets]    =   regexp(BIN(iBin).expression(jLogES), '.+?\d+.*?}', 'match');
                % Detects special codes as * (all) and ~* (none)
                [spCodes specialCodesToken]    =   regexp(BIN(iBin).expression(jLogES), '(\~)*(\*+)', 'match','tokens');  % 21 febrero 2008!

                if ~isempty(segmentedBrackets{:}) && isempty(spCodes{:})  % The sequencer exists and is numeric
                        nSequencer = length(segmentedBrackets{:});             % Number of sequencers
                        auxt1pre = [];
                        auxt2pre = [];

                        for kSequencer = 1:nSequencer                          % Sequencer's loop
                                %
                                % extract each sequencer from one to nSequencer, one by one
                                %
                                sequencerString = char(segmentedBrackets{1}(kSequencer));
                                sequencerString = regexprep(sequencerString, ':rt<".*?">', ''); % hides reaction time (if any)

                                %
                                % Check if sequencer are simple or they
                                % have more complex conditions like time,
                                % flags, etc.

                                storeCode   = [];
                                storeSign   = [];
                                catch_IC    = [];
                                catch_Sign  = []; % catchs interval code ans sign

                                %
                                % First, searchs event code given as an interval e.g {1-10;30-40}
                                % Careful: There might be several of these
                                % defined intervals
                                %
                                [intervalString intervalToken] = regexp(sequencerString, '[{;,>](\~*)(\d+)\-(\d+)',...
                                        'match','tokens'); % 4 February: added on">"

                                if ~isempty(intervalToken)   % Only if intervals were found
                                        coluInTok = size(intervalToken,2);
                                        pointerInterval   = 1;
                                        dispcond( condition, 'There is a specified interval x !')

                                        % generate all number given by interval
                                        for k = 1:coluInTok
                                                intervalCodeStart = str2double(cell2mat(intervalToken{k}(2)));
                                                intervalCodeStop  = str2double(cell2mat(intervalToken{k}(3)));

                                                if intervalCodeStart == intervalCodeStop
                                                        catch_IC(pointerInterval) = intervalCodeStart;

                                                        if strcmp(char(intervalToken{k}(1)),'~')    % There is negation
                                                                catch_Sign(pointerInterval) = 0;
                                                        else                                        % There is no negation
                                                                catch_Sign(pointerInterval) = 1;
                                                        end

                                                        pointerInterval = pointerInterval + 1;
                                                        dispcond( condition, 'but the values are equals...')

                                                elseif intervalCodeStart < intervalCodeStop

                                                        intervalArray = intervalCodeStart : intervalCodeStop;
                                                        catch_IC( pointerInterval:length(intervalArray) + pointerInterval - 1 ) = intervalArray;

                                                        if strcmp(char(intervalToken{k}(1)),'~')  % There is negation
                                                                catch_Sign( pointerInterval:length(intervalArray) + pointerInterval - 1 ) = zeros(1,length(intervalArray));
                                                        else                             % THere is no negation
                                                                catch_Sign( pointerInterval:length(intervalArray) + pointerInterval - 1 ) = ones(1,length(intervalArray));
                                                        end

                                                        pointerInterval = pointerInterval + length(intervalArray);
                                                        dispcond( condition, 'building a list from the specified interval !')

                                                else
                                                        fprintf('Fatal Error: bad ordered interval at BIN:%g\n', iBin);
                                                        isparsednumerically = 0;  % Numeric parsin was not approved
                                                        return
                                                end
                                                storeCode = catch_IC;      %28 feb 2008
                                                storeSign = catch_Sign;              %28 feb 2008
                                        end
                                        % sort generated event code list
                                        [auxMatrix] = cat(1,storeCode, storeSign)';
                                        [MatSorted] = sortrows(auxMatrix,1);
                                        storeCode   = MatSorted(:,1)';
                                        storeSign   = MatSorted(:,2)';
                                        storeTime   = repmat(-1,length(storeCode),2);
                                        storeFlag   = repmat(uint16(0),1, length(storeCode));
                                        storeFmask  = repmat(uint16(0),1, length(storeCode));
                                        storeWrite  = repmat(uint16(0),1, length(storeCode));
                                        storeWmask  = repmat(uint16(0),1, length(storeCode));
                                        % replacements...
                                        sequencerString = regexprep(sequencerString, '>(\~*)(\d+)\-(\d+)', '>1000000');
                                        sequencerString = regexprep(sequencerString, '(\~*)(\d+)\-(\d+):f', '2000000:f');
                                        sequencerString = regexprep(sequencerString, '(\~*)(\d+)\-(\d+):w', '3000000:w');
                                        sequencerString = regexprep(sequencerString, '(\~*)(\d+)\-(\d+)[,;}]', '');
                                        dispcond(condition,'Erasing Interval Expressions......')
                                else
                                        storeCode  = []; % There is no defined interval eventcode
                                        storeSign  = [];
                                        storeTime  = [];
                                        storeFlag  = [];
                                        storeFmask = [];
                                        storeWrite = [];
                                        storeWmask = [];
                                        dispcond( condition, 'there was not a specified interval ...')
                                end

                                %
                                % Second, detect discrete format for event code e.g {1;2;3;4} or {1,2,3,4}
                                %
                                [ma2 cod] = regexp(sequencerString, '[{;,](\~*)(\d+)(?=[;,}])',...
                                        'match','tokens') ;

                                if ~isempty(cod)
                                        dispcond( condition, 'There is a discrete event code!')
                                        auxSign = zeros(1,length(ma2));
                                        catchSingleCode    = zeros(1,length(ma2));

                                        for m = 1:length(ma2) % For each unique code detected
                                                dispcond( condition, 'extracting discrete event code!')
                                                catchSingleCode(m) = str2double(cell2mat(cod{m}(2)));

                                                if strcmp(char(cod{m}(1)),'~')  % There is negation
                                                        auxSign(m) = 0;
                                                else
                                                        auxSign(m) = 1;         % There is no negation
                                                end
                                        end

                                        storeCode   = cat(2,storeCode, catchSingleCode);
                                        storeSign   = cat(2,storeSign, auxSign);
                                        [auxMatrix] = cat(1,storeCode, storeSign)';
                                        [MatSorted] = sortrows(auxMatrix,1);
                                        dispcond( condition, 'Showing "simple" Values of complex expressions:',...
                                                storeCode, storeSign, storeTime, storeFlag, storeWrite)
                                        storeCode   = MatSorted(:,1)';
                                        storeSign   = MatSorted(:,2)';
                                        storeTime   =  repmat(-1,   length(storeCode),2);
                                        storeFlag   =  repmat(uint16(0),1, length(storeCode));
                                        storeFmask  =  repmat(uint16(0),1, length(storeCode));
                                        storeWmask  =  repmat(uint16(0),1, length(storeCode));
                                        storeWrite  =  repmat(uint16(0),1, length(storeCode));

                                        %**********************************************************
                                else
                                        dispcond( condition, 'there was not a discret complex interval ...')
                                end

                                %
                                % Now it searches for time specification and/or flag
                                %
                                dispcond( condition, 'EXTRACTING T, F, y W')
                                % All possible combinations!!!

                                expp =  cell(1);
                                expp{1}  = 't<(\d+)+-(\d+)+>(\~)*(\d+)+:f<(\w+)+>:w<(\w+)+>';    % 1  full-flag full-write
                                expp{2}  = 't<(\d+)+-(\d+)+>(\~)*(\d+)+:f<(\w+)+>(\TONGO)*';     % 2  full-flag ----------
                                expp{3}  = 't<(\d+)+-(\d+)+>(\~)*(\d+)+(\TONGO)*:w<(\w+)+>';     % 3  --------- full-write
                                expp{4}  = 't<(\d+)+-(\d+)+>(\~)*(\d+)+(\TONGO)*(\TONGO)*';      % 4  --------- ----------
                                expp{5}  = '(\TONGO)*(\TONGO)*(\~)*(\d+):f<(\w+)+>:w<(\w+)+>';   % 5  full-flag full-write
                                expp{6}  = '(\TONGO)*(\TONGO)*(\~)*(\d+):f<(\w+)+>(\TONGO)*';    % 6  full-flag ----------
                                expp{7}  = '(\TONGO)*(\TONGO)*(\~)*(\d+)+(\TONGO)*:w<(\w+)+>';   % 7  --------- full-write

                                expp{8}  = 't<(\d+)+-(\d+)+>(\~)*(\d+)+:fa<(\w+)+>:wa<(\w+)+>';  % 8  low-flag  low-write
                                expp{9}  = 't<(\d+)+-(\d+)+>(\~)*(\d+)+:fa<(\w+)+>:wb<(\w+)+>';  % 9  low-flag  high-write
                                expp{10} = 't<(\d+)+-(\d+)+>(\~)*(\d+)+:fb<(\w+)+>:wa<(\w+)+>';  % 10  high-flag low-write
                                expp{11} = 't<(\d+)+-(\d+)+>(\~)*(\d+)+:fb<(\w+)+>:wb<(\w+)+>';  % 11 high-flag high-write

                                expp{12} = 't<(\d+)+-(\d+)+>(\~)*(\d+)+:fa<(\w+)+>(\TONGO)*';    % 12 low-flag  ----------
                                expp{13} = 't<(\d+)+-(\d+)+>(\~)*(\d+)+:fb<(\w+)+>(\TONGO)*';    % 13 high-flag ----------
                                expp{14} = 't<(\d+)+-(\d+)+>(\~)*(\d+)+(\TONGO)*:wa<(\w+)+>';    % 14 --------- low-write
                                expp{15} = 't<(\d+)+-(\d+)+>(\~)*(\d+)+(\TONGO)*:wb<(\w+)+>';    % 15 --------- high-write

                                expp{16} = '(\TONGO)*(\TONGO)*(\~)*(\d+):fa<(\w+)+>:wa<(\w+)+>'; % 16 low-flag low-write
                                expp{17} = '(\TONGO)*(\TONGO)*(\~)*(\d+):fa<(\w+)+>:wb<(\w+)+>'; % 17 low-flag high-write
                                expp{18} = '(\TONGO)*(\TONGO)*(\~)*(\d+):fb<(\w+)+>:wa<(\w+)+>'; % 18 high-flag low-write
                                expp{19} = '(\TONGO)*(\TONGO)*(\~)*(\d+):fb<(\w+)+>:wb<(\w+)+>'; % 19 high-flag high-write

                                expp{20} = '(\TONGO)*(\TONGO)*(\~)*(\d+):fa<(\w+)+>(\TONGO)*';   % 20 low-flag  ---------
                                expp{21} = '(\TONGO)*(\TONGO)*(\~)*(\d+):fb<(\w+)+>(\TONGO)*';   % 21 high-flag ---------
                                expp{22} = '(\TONGO)*(\TONGO)*(\~)*(\d+)+(\TONGO)*:wa<(\w+)+>';  % 22 --------- low-write
                                expp{23} = '(\TONGO)*(\TONGO)*(\~)*(\d+)+(\TONGO)*:wb<(\w+)+>';  % 23 --------- high-write

                                auxstr   = sequencerString;
                                countexp = zeros(1,length(expp));  % January 28,2008
                                istimedetected = 0;   % boolean
                                isflagdetected = 0;   % boolean
                                isWritedetected= 0;   % boolean

                                for e = 1:length(expp)

                                        [mtr rtimes] = regexp(auxstr, expp{e},'match','tokens');

                                        if ~isempty(mtr)
                                                countexp(e) = length(mtr); % Test the take-all event (*)

                                                for n = 1:countexp(e)

                                                        catchComplexExpression = [rtimes{n}(1) rtimes{n}(2) rtimes{n}(4) rtimes{n}(5)...
                                                                rtimes{n}(6)];

                                                        t1 = str2double(cell2mat(catchComplexExpression(1)));
                                                        t2 = str2double(cell2mat(catchComplexExpression(2)));
                                                        captureCode = str2double(cell2mat(catchComplexExpression(3)));

                                                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                                        % Capture Flag (f, fa, and fb)
                                                        %
                                                        strflag = catchComplexExpression{4};
                                                        lestrf  = length(strflag);

                                                        if ismember_bc2(e,[1 2 5 6]) % full flag (16bits)
                                                                if lestrf > 16
                                                                        fprintf('Fatal Error: BIN %g. 16-bit Flag (f<>) format was violated',iBin)
                                                                        isparsednumerically = 0;  % Numeric parsing was not approved
                                                                        return
                                                                elseif lestrf <= 16 && lestrf >= 1

                                                                        %  Steve's option: always 'x' replaces blanks
                                                                        %  Example:  Flag = 10011  means xxxxxxxxxxx10011
                                                                        %  Example:  Flag = x10011  means xxxxxxxxxxx10011
                                                                        %
                                                                        if strflag(1)=='1' || strflag(1)=='0'
                                                                                strflag = [repmat('x',1, 16-lestrf) strflag];
                                                                        else
                                                                                strflag = [repmat('x',1, 16-lestrf) strflag];
                                                                        end
                                                                end
                                                        elseif ismember_bc2(e,[8 9 12 16 17 20]) % fa = low flag  = LSB byte = artifact flag

                                                                if lestrf > 8
                                                                        fprintf('Fatal Error: BIN %g. 8-bit Flag (fa<>) format was violated',iBin)
                                                                        isparsednumerically = 0;  % Numeric parsing was not approved
                                                                        return
                                                                elseif lestrf <= 8 && lestrf >= 1

                                                                        %  Steve's option: always 'x' replaces blanks
                                                                        %  Example:  Flag = 10011   means xxxxxxxxxxx10011
                                                                        %  Example:  Flag = x10011  means xxxxxxxxxxx10011
                                                                        %
                                                                        if strflag(1)=='1' || strflag(1)=='0'
                                                                                strflag = [repmat('x',1, 16-lestrf) strflag];
                                                                        else
                                                                                strflag = [repmat('x',1, 16-lestrf) strflag];
                                                                        end
                                                                end
                                                        elseif ismember_bc2(e,[10 11 13 18 19 21]) % fb = high flag  = MSB byte = user? flag

                                                                if lestrf > 8
                                                                        fprintf('Fatal Error: BIN %g. 8-bit Flag (fb<>) format was violated',iBin)
                                                                        isparsednumerically = 0;  % Numeric parsing was not approved
                                                                        return
                                                                elseif lestrf <= 8 && lestrf >= 1

                                                                        %  Steve's option: always 'x' replaces blanks
                                                                        %  Example:  Flag = 10011   means xxx10011xxxxxxxx
                                                                        %  Example:  Flag = x10011  means xxx10011xxxxxxxx
                                                                        %
                                                                        if strflag(1)=='1' || strflag(1)=='0'
                                                                                strflag = [repmat('x',1, 8-lestrf) strflag repmat('x',1, 8)];
                                                                        else
                                                                                strflag = [repmat('x',1, 8-lestrf) strflag repmat('x',1, 8)];
                                                                        end
                                                                end
                                                        end

                                                        [matf, tokf, indxf] = regexpi(strflag, '0|1','match','tokens','start');

                                                        if isempty(indxf)
                                                                captureFlag    = uint16(bin2dec(strflag));
                                                                maskFlagstring = '1111111111111111';
                                                                captFlagMask   = uint16(bin2dec(maskFlagstring)); % Flag's Logic Mask
                                                        else
                                                                strflag        =  regexprep(strflag, 'x','0','ignorecase');
                                                                captureFlag    = uint16(bin2dec(strflag));
                                                                maskFlagstring = repmat('0', 1, 16);
                                                                maskFlagstring(indxf) = '1';
                                                                captFlagMask   = uint16(bin2dec(maskFlagstring)); % Flag's Logic Mask
                                                        end

                                                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                                        % Capture Write and masking  10-27-2008
                                                        %
                                                        strWrite = catchComplexExpression{5};
                                                        lestrs  = length(strWrite);

                                                        if ismember_bc2(e,[1 3 5 7]) % full Write (16bits)
                                                                if lestrs > 16
                                                                        fprintf('Fatal Error: BIN %g. 16-bit Write (w<>) format was violated',iBin)
                                                                        isparsednumerically = 0;  % Numeric parsing was not approved
                                                                        return
                                                                elseif lestrs <= 16 && lestrs >= 1

                                                                        %  Steve's option: always 'x' replaces blanks
                                                                        %  Example:  Flag = 10011  means xxxxxxxxxxx10011
                                                                        %  Example:  Flag = x10011  means xxxxxxxxxxx10011
                                                                        %
                                                                        if strWrite(1)=='1' || strWrite(1)=='0'
                                                                                strWrite = [repmat('x',1, 16-lestrs) strWrite];
                                                                        else
                                                                                strWrite = [repmat('x',1, 16-lestrs) strWrite];
                                                                        end
                                                                end
                                                        elseif ismember_bc2(e,[8 10 14 16 18 22]) % wa = low write  = LSB byte = artifact write

                                                                if lestrs > 8
                                                                        fprintf('Fatal Error: BIN %g. 8-bit Write (wa<>) format was violated',iBin)
                                                                        isparsednumerically = 0;  % Numeric parsing was not approved
                                                                        return
                                                                elseif lestrs <= 8 && lestrs >= 1

                                                                        %  Steve's option: always 'x' replaces blanks
                                                                        %  Example:  write = 10011   means xxxxxxxxxxx10011
                                                                        %  Example:  write = x10011  means xxxxxxxxxxx10011
                                                                        %
                                                                        if strWrite(1)=='1' || strWrite(1)=='0'
                                                                                strWrite = [repmat('x',1, 16-lestrs) strWrite];
                                                                        else
                                                                                strWrite = [repmat('x',1, 16-lestrs) strWrite];
                                                                        end
                                                                end
                                                        elseif ismember_bc2(e,[9 11 15 17 19 23]) % wb = high write  = MSB byte = user? write

                                                                if lestrs > 8
                                                                        fprintf('Fatal Error: BIN %g. 8-bit Write (wb<>) format was violated',iBin)
                                                                        isparsednumerically = 0;  % Numeric parsing was not approved
                                                                        return
                                                                elseif lestrs <= 8 && lestrs >= 1

                                                                        %  Steve's option: always 'x' replaces blanks
                                                                        %  Example:  write = 10011   means xxx10011xxxxxxxx
                                                                        %  Example:  write = x10011  means xxx10011xxxxxxxx
                                                                        %
                                                                        if strWrite(1)=='1' || strWrite(1)=='0'
                                                                                strWrite = [repmat('x',1, 8-lestrs) strWrite repmat('x',1, 8)];
                                                                        else
                                                                                strWrite = [repmat('x',1, 8-lestrs) strWrite repmat('x',1, 8)];
                                                                        end
                                                                end
                                                        end

                                                        [mats, toks, indxs] = regexpi(strWrite, '0|1','match','tokens','start');

                                                        if isempty(indxs)
                                                                captureWrite    = uint16(bin2dec(strWrite));
                                                                maskWriteString = '1111111111111111';
                                                                captWriteMask   = uint16(bin2dec(maskWriteString)); % Write's Logic Mask
                                                        else
                                                                strWrite        =  regexprep(strWrite, 'x','0','ignorecase');
                                                                captureWrite    = uint16(bin2dec(strWrite));
                                                                maskWriteString = repmat('0', 1, 16);
                                                                maskWriteString(indxs) = '1';
                                                                captWriteMask   = uint16(bin2dec(maskWriteString)); % Write's Logic Mask
                                                        end
                                                        if ~isempty(t1) && ~isempty(t2)
                                                                if t1 >= t2
                                                                        fprintf('\nFatal Error: BIN %g, Sequencer %g, from "%s", has a erroneous time condition range\n',...
                                                                                iBin, kSequencer, logEventSelector{jLogES})
                                                                        isparsednumerically = 0;  % Numeric parsing was not approved
                                                                        return
                                                                end

                                                                if countexp(e) > 1
                                                                        dispcond(condition,'Hey! you are using more than one definition of time')
                                                                end
                                                                storeTime      = cat(1,storeTime, [t1 t2]);
                                                                t1catch        = t1;
                                                                t2catch        = t2;
                                                                istimedetected = 1;
                                                        else
                                                                if captureCode ~= 1000000 && captureCode ~= 2000000 && captureCode ~= 3000000
                                                                        storeTime = cat(1,storeTime, [-1 -1]);
                                                                end
                                                        end
                                                        if ~isempty(captureCode)
                                                                if captureCode ~= 1000000 && captureCode ~= 2000000 && captureCode ~= 3000000
                                                                        storeCode = cat(2,storeCode, captureCode);
                                                                        if char(rtimes{n}(3))=='~' % There is negation
                                                                                storeSign = cat(2,storeSign, 0);
                                                                        else                       % There is no negation
                                                                                storeSign = cat(2,storeSign, 1);
                                                                        end
                                                                end
                                                        else
                                                                fprintf('Fatal Error: expression without code!!! \n');
                                                                isparsednumerically = 0;  % Numeric parsing was not approved
                                                                return
                                                        end
                                                        if ~isempty(captureFlag)
                                                                if captureFlag > 65535
                                                                        fprintf('Fatal Error: BIN %g. 16-bit Flag format was violated',iBin)
                                                                        isparsednumerically = 0;  % Numeric parsing was not approved
                                                                        return
                                                                else
                                                                        flagcatch      = captureFlag;
                                                                        fmaskcatch     = captFlagMask;
                                                                        isflagdetected = 1;
                                                                        if captureCode ~= 2000000
                                                                                storeFlag   = cat(2,storeFlag, captureFlag );
                                                                                storeFmask  = cat(2,storeFmask, captFlagMask );
                                                                        end
                                                                end
                                                        else
                                                                if captureCode ~= 1000000 && captureCode ~= 2000000 && captureCode ~= 3000000
                                                                        storeFlag   = cat(2,storeFlag, 0);
                                                                        storeFmask  = cat(2,storeFmask, 0 );
                                                                end
                                                        end
                                                        if ~isempty(captureWrite)
                                                                if captureWrite > 65535
                                                                        fprintf('Fatal Error: BIN %g. 16-bit Write format was violated',iBin)
                                                                        isparsednumerically = 0;  % Numeric parsing was not approved
                                                                        return
                                                                else
                                                                        Writecatch     = captureWrite;
                                                                        smaskcatch     = captWriteMask;
                                                                        isWritedetected  = 1;
                                                                        if captureCode ~= 3000000
                                                                                storeWrite   = cat(2,storeWrite, captureWrite);
                                                                                storeWmask   = cat(2,storeWmask, captWriteMask );
                                                                        end
                                                                end
                                                        else

                                                                if captureCode ~= 1000000 && captureCode ~= 2000000 && captureCode ~= 3000000
                                                                        storeWrite   = cat(2,storeWrite, 0);
                                                                        storeWmask   = cat(2,storeWmask, 0 );
                                                                end
                                                        end
                                                end
                                                auxstr = regexprep(auxstr, expp{e}, '');
                                        else
                                                countexp(e) = 0;  %for future application
                                        end
                                end
                                if istimedetected
                                        lco = size(storeCode,2);
                                        storeTime = repmat([t1catch t2catch],lco,1);
                                end
                                if isflagdetected
                                        lco = size(storeCode,2);
                                        storeFlag  = repmat(flagcatch,1,lco);
                                        storeFmask = repmat(fmaskcatch,1,lco);
                                end
                                if isWritedetected
                                        lco = size(storeCode,2);
                                        storeWrite = repmat(Writecatch,1,lco);
                                        storeWmask = repmat(smaskcatch,1,lco);
                                end

                                %
                                % Structure
                                % Sorting by eventcode from lesser to
                                % greater including its time interval, flags, and writes
                                %
                                [ValMatrix]  = cat(1,num2cell(storeCode), num2cell(storeSign), num2cell(storeTime'), num2cell(storeFlag),...
                                        num2cell(storeFmask), num2cell(storeWrite), num2cell(storeWmask))';

                                auxforsort   = ValMatrix(:,1);
                                auxforsort   = [auxforsort{:}];
                                [indx, indx] = sort(auxforsort);
                                [MatrixSorted] = ValMatrix(indx,:);
                                storeCode      = cell2mat(MatrixSorted(:,1)');
                                storeSign      = cell2mat(MatrixSorted(:,2)');
                                storeTime      = cell2mat(MatrixSorted(:,3:4));
                                storeFlag      = cell2mat(MatrixSorted(:,5)');
                                storeFmask     = cell2mat(MatrixSorted(:,6)');
                                storeWrite     = cell2mat(MatrixSorted(:,7)');
                                storeWmask     = cell2mat(MatrixSorted(:,8)');

                                % Bad values filter
                                ecodeneg = find(storeCode <= 0, 1);
                                etimeneg = find(storeTime ~= -1 & storeTime <= 0 , 1);

                                if ~isempty(ecodeneg) || ~isempty(etimeneg)
                                        fprintf('Fatal Error: BIN %g has a bad numeric expression for event code or time range\n',iBin)
                                        isparsednumerically = 0;  % Numeric parsing was not approved
                                        return
                                end
                                %
                                % Detect if there is at least one negation and if so, everything is negated.
                                nWishedCodes  = nnz(storeSign); % number of nonzero elements in storeSign (nWishedCodes).
                                nDifferentSign = length(storeSign)-nWishedCodes;

                                if ~allowMixedSigns
                                    if nDifferentSign ~= length(storeSign) && nDifferentSign ~= 0
                                        storeSign = storeSign*0;
                                    end
                                end
                                if istimedetected
                                        t1now  = storeTime(1,1);
                                        t2now  = storeTime(1,2);
                                        sense =  jLogES-2;

                                        if ~isempty(auxt1pre)
                                                if ((sense*t2now)<(sense*auxt2pre)) || ((sense*t1now)<(sense*auxt1pre))
                                                        fprintf('\nFatal Error: BIN %g, Sequencer %g, from "%s", has a time condition range LESSER than previous one\n',...
                                                                iBin, kSequencer, logEventSelector{jLogES})
                                                        isparsednumerically = 0;  % No paso el parsing numerico
                                                        return
                                                end
                                        end
                                        auxt1pre = t1now;
                                        auxt2pre = t2now;
                                end

                                % Build a datastructure for the bin
                                % descriptor values
                                dispcond( condition, 'Showing Values (of complex expressions):',...
                                        storeCode, storeSign, storeTime, storeFlag, storeWrite)
                                %******************************************************
                                BIN(iBin).(logEventSelector{jLogES})(kSequencer).eventcode    =   single(storeCode);
                                BIN(iBin).(logEventSelector{jLogES})(kSequencer).eventsign    =   single(storeSign);
                                BIN(iBin).(logEventSelector{jLogES})(kSequencer).timecode     =   single(storeTime);
                                BIN(iBin).(logEventSelector{jLogES})(kSequencer).flagcode     =   uint16(storeFlag);
                                BIN(iBin).(logEventSelector{jLogES})(kSequencer).flagmask     =   uint16(storeFmask);
                                BIN(iBin).(logEventSelector{jLogES})(kSequencer).writecode    =   uint16(storeWrite);
                                BIN(iBin).(logEventSelector{jLogES})(kSequencer).writemask    =   uint16(storeWmask);
                        end  % Sequencer's loop

                elseif   isempty(segmentedBrackets{:}) && ~isempty(spCodes{:})    % The sequencer is not numeric

                        if char(specialCodesToken{1,1}{1,1}(1)) == '~'
                                BIN(iBin).(logEventSelector{jLogES}).eventcode =   -13;  % Code is a non-value!
                                BIN(iBin).(logEventSelector{jLogES}).eventsign =   0;
                        else
                                BIN(iBin).(logEventSelector{jLogES}).eventcode =   -7;   % Code is a non-value!
                                BIN(iBin).(logEventSelector{jLogES}).eventsign =   1;
                        end

                        BIN(iBin).(logEventSelector{jLogES-1})            =  []; % delete the prehome
                        BIN(iBin).(logEventSelector{jLogES}).timecode     =  single([-1 -1]);
                        BIN(iBin).(logEventSelector{jLogES}).flagcode     =  uint16(0);
                        BIN(iBin).(logEventSelector{jLogES}).flagmask     =  uint16(65535);
                        BIN(iBin).(logEventSelector{jLogES}).writecode    =  uint16(0);
                        BIN(iBin).(logEventSelector{jLogES}).writemask    =  uint16(65535);
                        BIN(iBin).(logEventSelector{jLogES+1})            =  [];  % Delete the posthome
                        isworkingFine = 0;  % Don't continue with the post home because the prehome was rejected.

                else   % No numeric sequencer and no special sequencer (~* or *)
                        BIN(iBin).(logEventSelector{jLogES}) =  []; % If there is no expression for the field, then it will be empty
                end
                jLogES = jLogES + 1;  % Field Number (pre, at, and post)
        end   % LES's loop
        % name the bin been processed
        BIN(iBin).namebin  =   ['BIN ' num2str(iBin)];
end % Bin's loop

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

% Feb 6, 2019
groundTruthFile = 'C:\NeuroSLAM_Datasets\02_NeuroSLAM_Groudtruth\01_SynPerData_GT.txt';
expMapFile = 'C:\NeuroSLAM_Datasets\03_NeuroSLAM_Experiments_Results\SynPerData\01_exp_map_ml.txt';
% plot_3d_multilayer_experience_map(groundTruthFile, expMapFile, xExpMapScaling, yExpMapScaling, zExpMapScaling, xExpMapTrans, yExpMapTrans, zExpMapTrans, xGtScaling, yGtScaling, zGtScaling)
plot_3d_multilayer_experience_map(groundTruthFile, expMapFile, 0.1, 0.1, 0.1, 1.2, 1.3, 0, 21.5,23, 23);

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function [Fo,Oq,Oqval,DEOPA,goodperiods,OqS,OqvalS,DEOPAS,goodperiodsS,simppeak,dSIG,SdSIG] = ...
    FO(COEF,method,propthresh,resampC,maxF,SIG,FS)

% FO: detection of Fundamental frequency and Open quotient on the basis of the
% derivative of the electroglottographic signal.
% This is the script that calls the various functions for analysis of the
% electroglottographic signal.
%
% The sequence of operations is as follows:
% - the function CRO detects the points where the DEGG signal crosses the threshold
% set for peak detection, and calls the function RIM which determines where peaks
% begin and end
%
% - the function AMPOS reinterpolates the signal within the peak regions (to
% obtain a more accurate measurement of peak height) and calls the function
% PEAKSHAPE to detect the number of summits and the ratio of the amplitude of
% the summits
%
% - the function DETECTMULT treats double peaks of another kind: double peaks
% that are almost completely distinct, and belong in different "peak bunches"
% (regions created by CRO), but which are too close for their Fo to make sense,
% and are thus considered as making up one double/multiple peak.
% There are thus two vectors containing information on double peaks: verif, for
% peaks which are unique at threshold level and fork at a higher level; and
% doublepeaks for peaks which fork below threshold level.
%
% - the function SIMP calculates an approximative Fo value.

% retrieving user choice for smoothing step. Reminder: 0 means no smoothing; 1
% means smoothing 1 point to the left and right, i.e. over 3 data points; 2
% means smoothing 2 points to the left and right, i.e. over 5 data points, etc.
C_SMOO = COEF(2);

% retrieving user choice for amplitude threshold setting method
toggle = COEF(3);

% retrieving user choice for amplitude threshold
AMP = COEF(4);

% Placing the sampling rate in first slot of <COEF> vector
COEF(1) = FS;

% Sampling period: (T for: period, s for: sampling)
Ts = 1 / COEF(1);

% check that the signal is Mono; if it is stereo: de-multiplexing it
[S,P] = size(SIG);
%%%% calculating the derivative of the signal
% <dSIG> is the first derivative of <SIG>. No smoothing is done at this stage, as the
% detection of the closing peaks is best done using the unsmoothed signal. If a
% smoothing step has been specified, it is used later: in the detection of opening peaks.
dSIG = [];
for w = 1 : length(SIG) - 1
   dSIG (w) = SIG (w + 1) - SIG (w);
end

%%%% smoothing the derivative of the signal
if C_SMOO > 0
    disp('Smoothing signal...')
    SdSIG = smoo(dSIG,C_SMOO);
    disp('Smoothed.')
else
    SdSIG = dSIG;
end

% detecting points where the DEGG signal crosses the threshold; this is done
% using the smoothed signal, otherwise the number of crossings may be very high
% due to a saw-like shape of the signal
[rims] = CRO(SdSIG,COEF);

% detection of the amplitude and position of the peaks (on the reinterpolated signal)
[Tgci,valid,validtime,numberofpeaks] = AMPOS(SIG,rims,resampC,Ts,method,propthresh);

% detection of double / multiple closing peaks ("peak clusters")
doublepeaks = detectmult(Tgci,maxF);

% choice of closing peaks, creating a matrix with only one maximum per cycle.
[simppeak,BUND] = simp(Tgci,doublepeaks,method);

%%%%%%%%%%% Calculating Fo on the basis of the unique peaks
Fo = [];
for w = 1:length(simppeak(:,1)) - 1
    Fo(w) = 1 / (simppeak(w + 1,1) - simppeak(w,1));
end

% Retrieving opening peaks where possible on unsmoothed signal
[Oq,Oqval,DEOPA,goodperiods] = OP(simppeak,SIG,dSIG,valid,doublepeaks,COEF,resampC,Ts,method);

%%% Second call to the function OP, with smoothed signal
[OqS,OqvalS,DEOPAS,goodperiodsS] = OP(simppeak,SIG,SdSIG,valid,doublepeaks,COEF,resampC,Ts,method);

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
%VOCEVALSEG Evaluates a set of segmentation results.
% VOCEVALSEG(VOCopts,ID); prints out the per class and overall
% segmentation accuracies. Accuracies are given using the intersection/union
% metric:
%   true positives / (true positives + false positives + false negatives)
%
% [ACCURACIES,AVACC,CONF] = VOCEVALSEG(VOCopts,ID) returns the per class
% percentage ACCURACIES, the average accuracy AVACC and the confusion
% matrix CONF.
%
% [ACCURACIES,AVACC,CONF,RAWCOUNTS] = VOCEVALSEG(VOCopts,ID) also returns
% the unnormalised confusion matrix, which contains raw pixel counts.
function [accuracies,avacc,conf,rawcounts] = tvg_VOCevalseg(VOCopts,id)

% image test set
[gtids,t]=textread(sprintf(VOCopts.seg.imgsetpath,VOCopts.testset),'%s %d');

% number of labels = number of classes plus one for the background
num = VOCopts.nclasses+1;
confcounts = zeros(num);
count=0;
tic;
for i=1:length(gtids)
    % display progress
    if toc>1
        fprintf('test confusion: %d/%d\n',i,length(gtids));
        drawnow;
        tic;
    end

    imname = gtids{i};

    % ground truth label file
    gtfile = sprintf(VOCopts.seg.clsimgpath,imname);
    [gtim,map] = imread(gtfile);
    gtim = double(gtim);

    % results file
    resfile = sprintf(VOCopts.seg.clsrespath,id,VOCopts.testset,imname);
    [resim,map] = imread(resfile);
    resim = double(resim);

    % Check validity of results image
    maxlabel = max(resim(:));
    if (maxlabel>VOCopts.nclasses),
        error('Results image ''%s'' has out of range value %d (the value should be <= %d)',imname,maxlabel,VOCopts.nclasses);
    end

    szgtim = size(gtim); szresim = size(resim);
    if any(szgtim~=szresim)
        error('Results image ''%s'' is the wrong size, was %d x %d, should be %d x %d.',imname,szresim(1),szresim(2),szgtim(1),szgtim(2));
    end

    %pixel locations to include in computation
    locs = gtim<255;

    % joint histogram
    sumim = 1+gtim+resim*num;
    hs = histc(sumim(locs),1:num*num);
    count = count + numel(find(locs));
    confcounts(:) = confcounts(:) + hs(:);
end

% confusion matrix - first index is true label, second is inferred label
%conf = zeros(num);
conf = 100*confcounts./repmat(1E-20+sum(confcounts,2),[1 size(confcounts,2)]);
rawcounts = confcounts;

% Percentage correct labels measure is no longer being used.  Uncomment if
% you wish to see it anyway
%overall_acc = 100*sum(diag(confcounts)) / sum(confcounts(:));
%fprintf('Percentage of pixels correctly labelled overall: %6.3f%%\n',overall_acc);

accuracies = zeros(VOCopts.nclasses,1);
fprintf('Accuracy for each class (intersection/union measure)\n');
for j=1:num

   gtj=sum(confcounts(j,:));
   resj=sum(confcounts(:,j));
   gtjresj=confcounts(j,j);
   % The accuracy is: true positive / (true positive + false positive + false negative)
   % which is equivalent to the following percentage:
   accuracies(j)=100*gtjresj/(gtj+resj-gtjresj);

   clname = 'background';
   if (j>1), clname = VOCopts.classes{j-1};end;
   fprintf('  %14s: %6.3f%%\n',clname,accuracies(j));
end
accuracies = accuracies(1:end);
avacc = mean(accuracies);
fprintf('-------------------------\n');
fprintf('Average accuracy: %6.3f%%\n',avacc);

author: torrvision

repo: crfasrnn

save_path: github-repos/MATLAB/torrvision-crfasrnn


text:
function [] = aps_era_ECMWF_Python(timing,wheatherregstr)
% script for the creation of download python script files for ERA-I model
% downloads monthly data in splitted times
% input str: timing in the format 'YYYYMMDDHH'
%
%     Copyright (C) 2015
%     With permission of Hannes Bathke
%     Email: eedpsb@leeds.ac.uk or davidbekaert.com
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
% By Hannes Bathke -- January 2014
%
% modifications:
% DB    02/2014 Integrate in the APS toolbox
% DB    03/2014 Auto download included

year = timing(1:4);
month = timing(5:6);
day = timing(7:8);
htime = timing(9:10);

%mfield = calendar(str2double(year),str2double(month));
%nday = max(max(mfield));
datestr = [year,'-',month,'-',day];

pythonsc_path = [timing,'.py'];
fid = fopen(pythonsc_path,'w');
fprintf(fid,'#!/usr/bin/env python\n');
fprintf(fid,'\n');
fprintf(fid,'from ecmwfapi import ECMWFDataServer\n');
fprintf(fid,'\n');
fprintf(fid,'# To run this example, you need an API key\n');
fprintf(fid,'# available from https://api.ecmwf.int/v1/key/\n');
fprintf(fid,'\n');
fprintf(fid,'server = ECMWFDataServer()\n');
fprintf(fid,'\n');
fprintf(fid,'server.retrieve({\n');
fprintf(fid,'          ''dataset'' : "interim",\n');
fprintf(fid,'          ''class''   : "ei",\n');
fprintf(fid,'          ''type''    : "an",\n');
fprintf(fid,'          ''stream''  : "oper",\n');
fprintf(fid,'          ''levtype'' : "pl",\n');
fprintf(fid,'          ''levelist'': "1/2/3/5/7/10/20/30/50/70/100/125/150/175/200/225/250/300/350/400/450/500/550/600/650/700/750/775/800/825/850/875/900/925/950/975/1000",\n');
fprintf(fid,'          ''param''   : "129.128/130.128/157.128/246.128",\n');
fprintf(fid,'          ''date''    : "%s",\n',datestr);
fprintf(fid,'          ''time''    : "%s",\n',htime);
fprintf(fid,'          ''step''    : "0",\n');
fprintf(fid,'          ''format''  : "netcdf",\n');
fprintf(fid,'          ''resol''   : "auto",\n');
fprintf(fid,'          ''area''    : "%s",\n',wheatherregstr);
fprintf(fid,'          ''grid''    : "0.75/0.75",\n');
fprintf(fid,['          ''target''  : "ggap' num2str(timing(1:10)) '00.nc",\n']);
fprintf(fid,'                })\n');

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
classdef TestFindCirclesGrid
    %TestFindCirclesGrid

    methods (Static)
        function test_1
            img = imread(fullfile(mexopencv.root(),'test','acircles_pattern.png'));
            img = cv.resize(img, 0.5, 0.5);
            sz = [4 11];
            [centers,found] = cv.findCirclesGrid(img, sz, 'SymmetricGrid',false);
            validateattributes(centers, {'cell'}, {'vector'});
            cellfun(@(v) validateattributes(v, {'numeric'}, ...
                {'vector', 'numel',2, 'real'}), centers);
            validateattributes(found, {'logical'}, {'scalar'});
            if found
                assert(numel(centers) == prod(sz));
            end
        end

        function test_2
            img = imread(fullfile(mexopencv.root(),'test','acircles_pattern.png'));
            sz = [4 11];
            [centers,found] = cv.findCirclesGrid(img, sz, 'SymmetricGrid',false, ...
                'BlobDetector',{'SimpleBlobDetector', 'MaxArea',100000});
            validateattributes(centers, {'cell'}, {'vector'});
            cellfun(@(v) validateattributes(v, {'numeric'}, ...
                {'vector', 'numel',2, 'real'}), centers);
            validateattributes(found, {'logical'}, {'scalar'});
            if found
                assert(numel(centers) == prod(sz));
            end
        end

        function test_error_argnum
            try
                cv.findCirclesGrid();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [ theta, phi ] = xyz_to_tp ( xyz )

%*****************************************************************************80
%
%% XYZ_TO_TP converts (X,Y,Z) to (Theta,Phi) coordinates.
%
%  Discussion:
%
%    Given an XYZ point, regard it as lying on a sphere of radius R,
%    centered at the origin, whose axis is the Z axis.
%
%    We assume that the actual value of R is of no interest, and do
%    not report it.  This is especially appropriate if the point is
%    expected to lie on the unit sphere, for instance.
%
%    THETA measures the "longitude" of the point, between 0 and 2 PI.
%
%    PHI measures the angle from the "north pole", between 0 and PI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real XYZ(3,1), the coordinates of a point in 3D.
%
%    Output, real THETA, PHI, the longitude and declination
%    of the point.
%
  r = norm ( xyz(1:3,1) );

  if ( r == 0.0 )
    theta = 0.0;
    phi = 0.0;
    return
  end

  phi = r8_acos ( xyz(3,1) / r );

  theta = r8_atan ( xyz(2,1), xyz(1,1) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_tutorial_natmeg2014_preprocessing

% WALLTIME 00:20:00
% MEM 5gb
% DEPENDENCY

% this script executes the MATLAB content from
% http://www.fieldtriptoolbox.org/tutorial/natmeg2014/timefrequency
%
% it corresponds to the wiki version of 7 October 2014

clear all
close all

cd(dccnpath('/home/common/matlab/fieldtrip/data/ftp/workshop/natmeg2014'));

cfg = [];
cfg.dataset    = 'oddball1_mc_downsampled.fif';
cfg.continuous = 'yes';
cfg.channel    = 'MEG*1';
cfg.viewmode   = 'vertical';
cfg.blocksize  = 1; % Length of data to display, in seconds
ft_databrowser(cfg);

cfg = [];
cfg.dataset = 'oddball1_mc_downsampled.fif';
cfg.channel = {'MEG*2', 'MEG*3'};
cfg.viewmode = 'vertical';
cfg.blocksize = 1;                             % Length of data to display, in seconds
ft_databrowser(cfg);

cfg = [];
cfg.dataset = 'oddball1_mc_downsampled.fif';
cfg.channel = 'EEG';
cfg.viewmode = 'vertical';
cfg.blocksize = 1;                             % Length of data to display, in seconds
cfg.preproc.demean = 'yes';                    % Demean the data before display
cfg.ylim = [-4e-6 4e-6];
ft_databrowser(cfg);

cfg = [];
cfg.dataset = 'oddball1_mc_downsampled.fif';

cfg.trialdef.prestim        = 1;
cfg.trialdef.poststim       = 1;
cfg.trialdef.std_triggers   = 1;
cfg.trialdef.stim_triggers  = [1 2]; % 1 for standard, 2 for deviant
cfg.trialdef.odd_triggers   = 2;
cfg.trialdef.rsp_triggers   = [256 4096];
cfg.trialfun                = 'trialfun_oddball_stimlocked';
cfg                         = ft_definetrial(cfg);

cfg.continuous              = 'yes';
cfg.hpfilter                = 'no';
cfg.detrend                 = 'no';
cfg.continuous              = 'yes';
cfg.demean                  = 'yes';
cfg.dftfilter               = 'yes';
cfg.dftfreq                 = [50 100];
cfg.channel                 = 'MEG';
data_MEG                    = ft_preprocessing(cfg);

if false
  % skip the interactive section
  % separately for magnetometers
  cfg               = [];
  cfg.metric        = 'zvalue';
  cfg.layout        = 'neuromag306all.lay';
  cfg.channel       = 'MEG*1';
  cfg.keepchannel   = 'yes';  % This keeps those channels that are not displayed in the data
  data_MEG_clean    = ft_rejectvisual(cfg,data_MEG);
  % separately for gradiometers
  cfg.channel = {'MEG*2','MEG*3'};
  data_MEG_clean    = ft_rejectvisual(cfg,data_MEG_clean);
else
  % simply copy it over
  data_MEG_clean = data_MEG;
end

cfg = [];
cfg.lpfilter        = 'yes';
cfg.lpfreq          = 25;
cfg.demean          = 'yes';
cfg.baselinewindow  = [-0.5 0];
data_MEG_filt       = ft_preprocessing(cfg,data_MEG_clean);

cfg = [];
cfg.trials          = find(data_MEG_filt.trialinfo(:,1) == 1);
ERF_standard        = ft_timelockanalysis(cfg,data_MEG_filt);

cfg.trials          = find(data_MEG_filt.trialinfo(:,1) == 2);
ERF_oddball         = ft_timelockanalysis(cfg,data_MEG_filt);

cfg = [];
cfg.operation = 'x1 - x2';
cfg.parameter = 'avg';
ERF_diff = ft_math(cfg, ERF_oddball, ERF_standard);

cfg = [];
cfg.fontsize = 6;
cfg.layout = 'neuromag306mag.lay';
cfg.ylim = [-2.5e-13 2.5e-13];
cfg.xlim = [-0.2 0.6];

figure;
ft_multiplotER(cfg, ERF_standard, ERF_oddball, ERF_diff );
legend({'Standard';'Oddball';'Difference'});

cfg = [];
cfg.fontsize = 6;
cfg.layout   = 'neuromag306mag.lay';
cfg.xlim     = [-0.2 0.6];
cfg.ylim     = [-3e-13 3e-13];
cfg.channel  = 'MEG0211';

figure;
ft_singleplotER(cfg, ERF_standard, ERF_oddball, ERF_diff);
legend({'Standard';'Oddball';'Difference'});

cfg                 = [];
cfg.layout          = 'neuromag306mag.lay'; % name will change
cfg.zlim            = [-3e-13 3e-13];
cfg.xlim            = [0.08 0.15];
cfg.style           = 'straight';
cfg.comment         = 'no';
cfg.marker          = 'off';
cfg.colorbar        = 'southoutside';

figure;
subplot(1,3,1);
ft_topoplotER(cfg, ERF_standard);
title('Standard');
axis tight

subplot(1,3,2);
ft_topoplotER(cfg, ERF_oddball);
title('Oddball');
axis tight

subplot(1,3,3);
ft_topoplotER(cfg, ERF_diff);
title('Difference');
axis tight

% Combine planar
cfg = [];
ERF_standard_cmb    = ft_combineplanar(cfg, ERF_standard);
ERF_oddball_cmb     = ft_combineplanar(cfg, ERF_oddball);
ERF_diff_cmb        = ft_combineplanar(cfg, ERF_diff);

cfg = [];
cfg.fontsize = 6;
cfg.layout   = 'neuromag306cmb.lay';
cfg.ylim     = [0 8e-12];
cfg.xlim     = [-0.2 0.6];

figure;
ft_multiplotER(cfg, ERF_standard_cmb, ERF_oddball_cmb, ERF_diff_cmb);
legend({'Standard?, ?Oddball?, ?Difference'});

cfg = [];
cfg.showlabels = 'yes';
cfg.fontsize   = 6;
cfg.layout     = 'neuromag306cmb.lay';
cfg.xlim       = [-0.2 0.6];
cfg.ylim       = [0 8e-12];
cfg.channel    = 'MEG0222+0223';

figure;
ft_singleplotER(cfg, ERF_standard_cmb, ERF_oddball_cmb, ERF_diff_cmb);
legend({'Standard?, ?Oddball?, ?Difference'});

cfg                 = [];
cfg.layout          = 'neuromag306cmb.lay'; % name will change
cfg.zlim            = 'zeromax';
cfg.xlim            = [0.08 0.15];
cfg.style           = 'straight';
cfg.comment         = 'no';
cfg.marker          = 'off';
cfg.colorbar        = 'southoutside';

figure;
subplot(1,3,1);
ft_topoplotER(cfg, ERF_standard_cmb);
title('Standard');
axis tight

subplot(1,3,2);
ft_topoplotER(cfg, ERF_oddball_cmb);
title('Deviant');
axis tight

subplot(1,3,3);
ft_topoplotER(cfg, ERF_diff_cmb);
title('Difference');
axis tight

cfg = [];
cfg.dataset = 'oddball1_mc_downsampled.fif';

cfg.trialdef.prestim        = 1;
cfg.trialdef.poststim       = 1;
cfg.trialdef.std_triggers   = 1;
cfg.trialdef.stim_triggers  = [1 2];
cfg.trialdef.odd_triggers   = 2;
cfg.trialdef.rsp_triggers   = [256 4096];
cfg.trialfun                = 'trialfun_oddball_stimlocked';
cfg                         = ft_definetrial(cfg);

cfg.continuous              = 'yes';
cfg.hpfilter                = 'no';
cfg.detrend                 = 'no';
cfg.continuous              = 'yes';
cfg.demean                  = 'yes';
cfg.dftfilter               = 'yes';
cfg.dftfreq                 = [50 100];
cfg.channel                 = 'EEG';

cfg.reref                   = 'yes'; % recorded with left mastoid
cfg.refchannel              = 'all';
data_EEG                    = ft_preprocessing(cfg);

if false
  % skip the interactive section
  cfg               = [];
  cfg.metric        = 'zvalue';
  cfg.layout        = 'natmeg_customized_eeg1005.lay';
  data_EEG_clean    = ft_rejectvisual(cfg,data_EEG);
else
  % simply copy the data over
  data_EEG_clean = data_EEG;
end

cfg = [];
cfg.lpfilter        = 'yes';
cfg.lpfreq          = 25;
cfg.demean          = 'yes';
cfg.baselinewindow  = [-0.5 0];
data_EEG_filt       = ft_preprocessing(cfg,data_EEG_clean);

cfg = [];
cfg.trials          = find(data_EEG_filt.trialinfo(:,1) == 1);
ERP_standard        = ft_timelockanalysis(cfg, data_EEG_filt);
cfg.trials          = find(data_EEG_filt.trialinfo(:,1) == 2);
ERP_oddball         = ft_timelockanalysis(cfg, data_EEG_filt);

cfg = [];
cfg.operation = 'x1 - x2';
cfg.parameter = 'avg';
ERP_diff = ft_math(cfg, ERP_oddball, ERP_standard);

cfg          = [];
cfg.fontsize = 6;
cfg.layout   = 'natmeg_customized_eeg1005.lay';
cfg.ylim     = [-3e-6 3e-6];
cfg.xlim     = [-0.2 0.6];

figure;
ft_multiplotER(cfg, ERP_standard, ERP_oddball, ERP_diff);

cfg            = [];
cfg.showlabels = 'yes';
cfg.fontsize   = 6;
cfg.layout     = 'natmeg_customized_eeg1005.lay';
cfg.xlim       = [-0.2 0.6];
cfg.ylim       = [-8e-6 8e-6];
cfg.channel    = 'EEG020';

figure;
ft_singleplotER(cfg, ERP_standard, ERP_oddball, ERP_diff);
legend({'Standard';'Oddball';'Difference'});

% Topo
cfg                 = [];
cfg.layout          = 'natmeg_customized_eeg1005.lay';
cfg.zlim            = [-3e-6 3e-6];
cfg.xlim            = [0.08 0.15];
cfg.style           = 'straight';
cfg.comment         = 'no';
cfg.marker          = 'off';
cfg.colorbar        = 'southoutside';

figure;
subplot(1,3,1);
ft_topoplotER(cfg,ERP_standard);
title('Standard');
axis tight

subplot(1,3,2);
ft_topoplotER(cfg,ERP_oddball);
title('Deviant');
axis tight

subplot(1,3,3);
ft_topoplotER(cfg,ERP_diff);
title('Difference');
axis tight

cfg                 = [];
cfg.method          = 'finite';
cfg.elec            = ERP_standard.elec;

scd_ERP_standard    = ft_scalpcurrentdensity(cfg, ERP_standard);
scd_ERP_oddball     = ft_scalpcurrentdensity(cfg, ERP_oddball);
scd_ERP_diff        = ft_scalpcurrentdensity(cfg, ERP_diff);

cfg                 = [];
cfg.layout          = 'natmeg_customized_eeg1005.lay'; % name will change
cfg.zlim            = 'maxabs';
cfg.xlim            = [0.08 0.15];
cfg.style           = 'straight';
cfg.comment         = 'no';
cfg.marker          = 'off';
cfg.colorbar        = 'southoutside';

figure;
subplot(1,3,1);
ft_topoplotER(cfg,scd_ERP_standard);
title('Standard');
axis tight;

subplot(1,3,2);
ft_topoplotER(cfg,scd_ERP_oddball);
title('Oddball');
axis tight;

subplot(1,3,3);
ft_topoplotER(cfg,scd_ERP_diff);
title('Difference');
axis tight;

cfg      = [];
data_all = ft_appenddata(cfg, data_MEG, data_EEG);

if false
  % skip the interactive section
  cfg = [];
  cfg.channel    = 'EEG';
  cfg.metric     = 'zvalue';
  cfg.keepchannel= 'yes';
  cfg.layout     = 'neuromag306all.lay';
  data_all_clean = ft_rejectvisual(cfg, data_all);

  cfg.channel    = 'MEGMAG';
  data_all_clean = ft_rejectvisual(cfg, data_all_clean);

  cfg.channel    = 'MEGGRAD';
  data_all_clean = ft_rejectvisual(cfg, data_all_clean);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function g = log(f, pref)
%LOG   Natural logarithm of a CHEBFUN.
%   LOG(F) returns the natural logarithm of F. If F has any roots in its domain,
%   then the representation is likely to be inaccurate.
%
% See also LOG1P, LOG2, LOG10, EXP, REALLOG.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

if ( nargin < 2 )
    pref = chebfunpref();
end

% Add breaks at the roots of f:
f = addBreaksAtRoots(f);

% Call COMPOSE():
g = compose(f, @(x) log(x), [], pref);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function gamma_nk = calc_gamma_E_step(N_nk, w_k)
%CALC_GAMMA Calculate gamma in E step while scaling both the denominator
% and numerator by the largest quantity in N.
% Input:
%   N_nk - logrithm of Gaussian evaluations (N by K)
%   w_k - prior probabilities (1 by K)
% Output:
%   gamma_nk - soft membership in E-step (N by K)
K1 = length(w_k);
N = size(N_nk, 1);
max_N_nk = max(N_nk, [], 2);
N_nk1 = exp(N_nk - repmat(max_N_nk, [1,K1])); % avoid too large negative

gamma_nk = (ones(N,1)*w_k) .* N_nk1;
gamma_nk = gamma_nk ./ repmat(sum(gamma_nk,2), 1, K1);

end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function plot_segments(P,Bal,fig,ms,seg1,seg2,seg3,seg4,seg5)

% Plots point cloud segments/subsets defined as subsets of cover sets.
% If the subsets intersect, then assiggnes the common points to the
% segments given first.
%
% Inputs
% P     Point cloud
% Bal   Cover sets. Bal = cover.ball
% fig   figure number
% seg1  Segment/subset 1, color blue
% seg2  (Optional) Segment/subset 2, color red
% seg3  (Optional) Segment/subset 3, color green
% seg4  (Optional) Segment/subset 4, color cyan
% seg5  (Optional) Segment/subset 5, color magenta

if nargin == 5
    S1 = unique(vertcat(Bal{seg1}));
    figure(fig)
    plot3(P(S1,1),P(S1,2),P(S1,3),'b.','Markersize',ms)
    axis equal
elseif nargin == 6
    S1 = unique(vertcat(Bal{seg1}));
    S2 = unique(vertcat(Bal{seg2}));
    S2 = setdiff(S2,S1);
    figure(fig)
    plot3(P(S1,1),P(S1,2),P(S1,3),'b.','Markersize',1.5*ms)
    hold on
    plot3(P(S2,1),P(S2,2),P(S2,3),'r.','Markersize',ms)
    axis equal
    hold off
elseif nargin == 7
    S1 = unique(vertcat(Bal{seg1}));
    S2 = unique(vertcat(Bal{seg2}));
    S3 = unique(vertcat(Bal{seg3}));
    S2 = setdiff(S2,S1);
    S3 = setdiff(S3,S1);
    S3 = setdiff(S3,S2);
    figure(fig)
    plot3(P(S1,1),P(S1,2),P(S1,3),'b.','Markersize',ms)
    hold on
    plot3(P(S2,1),P(S2,2),P(S2,3),'r.','Markersize',ms)
    plot3(P(S3,1),P(S3,2),P(S3,3),'g.','Markersize',ms)
    axis equal
    hold off
elseif nargin == 8
    S1 = unique(vertcat(Bal{seg1}));
    S2 = unique(vertcat(Bal{seg2}));
    S3 = unique(vertcat(Bal{seg3}));
    S4 = unique(vertcat(Bal{seg4}));
    S2 = setdiff(S2,S1);
    S3 = setdiff(S3,S1);
    S3 = setdiff(S3,S2);
    S4 = setdiff(S4,S1);
    S4 = setdiff(S4,S2);
    S4 = setdiff(S4,S3);
    figure(fig)
    plot3(P(S1,1),P(S1,2),P(S1,3),'b.','Markersize',ms)
    hold on
    plot3(P(S2,1),P(S2,2),P(S2,3),'r.','Markersize',ms)
    plot3(P(S3,1),P(S3,2),P(S3,3),'g.','Markersize',ms)
    plot3(P(S4,1),P(S4,2),P(S4,3),'c.','Markersize',ms)
    axis equal
    hold off
elseif nargin == 9
    S1 = unique(vertcat(Bal{seg1}));
    S2 = unique(vertcat(Bal{seg2}));
    S3 = unique(vertcat(Bal{seg3}));
    S4 = unique(vertcat(Bal{seg4}));
    S5 = unique(vertcat(Bal{seg5}));
    S2 = setdiff(S2,S1);
    S3 = setdiff(S3,S1);
    S3 = setdiff(S3,S2);
    S4 = setdiff(S4,S1);
    S4 = setdiff(S4,S2);
    S4 = setdiff(S4,S3);
    S5 = setdiff(S5,S1);
    S5 = setdiff(S5,S2);
    S5 = setdiff(S5,S3);
    S5 = setdiff(S5,S4);
    figure(fig)
    plot3(P(S1,1),P(S1,2),P(S1,3),'b.','Markersize',ms)
    hold on
    plot3(P(S2,1),P(S2,2),P(S2,3),'r.','Markersize',ms)
    plot3(P(S3,1),P(S3,2),P(S3,3),'g.','Markersize',ms)
    plot3(P(S4,1),P(S4,2),P(S4,3),'c.','Markersize',ms)
    plot3(P(S5,1),P(S5,2),P(S5,3),'m.','Markersize',ms)
    axis equal
    hold off
end

author: InverseTampere

repo: TreeQSM


text:
function [rd, boardName] = hdlcoder_ref_design_customization
% Reference design plugin registration file
% 1. The registration file with this name inside of a board plugin folder
%    will be picked up
% 2. Any registration file with this name on MATLAB path will also be picked up
% 3. The registration file returns a cell array pointing to the location of
%    the reference design plugins
% 4. The registration file also returns its associated board name
% 5. Reference design plugin must be a package folder accessible from
%    MATLAB path, and contains a reference design definition file

%   Copyright 2013-2014 The MathWorks, Inc.

rd = {...
    'AnalogDevices.pluto.plugin_rd_rx', ...
    'AnalogDevices.pluto.plugin_rd_tx', ...
    'AnalogDevices.pluto.plugin_rd_rxtx', ...
    };

boardName = 'AnalogDevices ADALM-PLUTO';

end

author: analogdevicesinc


text:
function varargout = trimMesh(varargin)
%TRIMMESH Reduce memory footprint of a polygonal mesh.
%
%   [V2, F2] = trimMesh(V, F)
%   Unreferenced vertices are removed.
%   Following functions are implemented for only numeric faces:
%       Duplicate vertices are removed.
%       Duplicate faces are removed.
%
%   Example
%     [V, F] = createIcosahedron;
%     F(13:20, :) = [];
%     [V2, F2] = trimMesh(V, F);
%     figure; drawMesh(V2, F2)
%     view(3); axis equal;
%     axis([-1 1 -1 1 0 2])
%
%   See also
%     meshes3d, clipMeshVertices

% ------
% Author: David Legland, oqilipo
% e-mail: david.legland@inra.fr
% Created: 2014-08-01,    using Matlab 8.3.0.532 (R2014a)
% Copyright 2014 INRA - Cepia Software Platform.

% parse input data
[vertices, faces] = parseMeshData(varargin{:});

if isnumeric(faces)
    % Delete duplicate vertices
    [tempVertices, ~, tempFaceVertexIdx] = unique(vertices, 'rows');
    tempFaces = tempFaceVertexIdx(faces);
    % Delete unindexed/unreferenced vertices
    usedVertexIdx = ismember(1:length(tempVertices),unique(tempFaces(:)));
    newVertexIdx = cumsum(usedVertexIdx);
    faceVertexIdx = 1:length(tempVertices);
    faceVertexIdx(usedVertexIdx) = newVertexIdx(usedVertexIdx);
    faceVertexIdx(~usedVertexIdx) = nan;
    tempFaces2 = faceVertexIdx(tempFaces);
    tempVertices2 = tempVertices(usedVertexIdx,:);
    % Delete duplicate faces
    [~, uniqueFaceIdx, ~] = unique(tempFaces2, 'rows');
    duplicateFaceIdx=~ismember(1:size(tempFaces2,1),uniqueFaceIdx);
    [vertices2, faces2] = removeMeshFaces(tempVertices2, tempFaces2, duplicateFaceIdx);
elseif iscell(faces)
    % identify vertices referenced by a face
    vertexUsed = false(size(vertices, 1), 1);
    for iFace = 1:length(faces)
        face = faces{iFace};
        vertexUsed(face) = true;
    end
    vertices2 = vertices(vertexUsed, :);
    % compute map from old index to new index
    inds = find(vertexUsed);
    newInds = zeros(size(vertices, 1), 1);
    for iIndex = 1:length(inds)
        newInds(inds(iIndex)) = iIndex;
    end
    % change labels of vertices referenced by faces
    faces2 = cell(1, length(faces));
    for iFace = 1:length(faces)
        faces2{iFace} = newInds(faces{iFace});
    end
else
    error('Unsupported format!')
end

% format output arguments
varargout = formatMeshOutput(nargout, vertices2, faces2);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function coeffs = fitPolynomialTransform2d(pts, ptsRef, degree)
%FITPOLYNOMIALTRANSFORM2D Coefficients of polynomial transform between two point sets.
%
%   COEFFS = fitPolynomialTransform2d(PTS, PTSREF, DEGREE)
%
%   Example
%
%   See also
%     polynomialTransform2d, fitAffineTransform2d

% ------
% Author: David Legland
% E-mail: david.legland@nantes.inra.fr
% Created: 2013-11-05, using Matlab 7.9.0.529 (R2009b)
% Copyright 2013-2022 INRA - Cepia Software Platform

%% Extract data

% ensure degree is valid
if nargin < 3
    degree = 3;
end

% polygon coordinates
xi = pts(:,1);
yi = pts(:,2);
nCoords = size(pts, 1);

% check inputs have same size
if size(ptsRef, 1) ~= nCoords
    error('fitPolynomialTransform2d:sizeError', ...
        'input arrays must have same number of points');
end

%% compute coefficient matrix

% number of coefficients of polynomial transform
nCoeffs = prod(degree + [1 2]) / 2;

% initialize matrix
A1 = zeros(nCoords, nCoeffs);

% iterate over degrees
iCoeff = 0;
for iDegree = 0:degree

    % iterate over binomial coefficients of a given degree
    for k = 0:iDegree
        iCoeff = iCoeff + 1;
        A1(:, iCoeff) = ones(nCoords, 1) .* power(xi, iDegree-k) .* power(yi, k);
    end
end

% concatenate matrix for both coordinates
A = kron(A1, [1 0;0 1]);

%% solve linear system that minimizes least squares

% create the vector of expected values
b = ptsRef';
b = b(:);

% solve the system
coeffs = (A \ b)';

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [position0,position,messages]= delineate3D(w1,w2,w3,intreg,timenew,position0,position,intervalo,position1,position2,position3,indexes,samp,sig,messages)
%
%
% [position0,position,A,V,A2,V2]= delineate3D(w1,w2,w3,intreg,scale,timenew,position0,position,intervalo,position1,samp)
% 3D multilead delineation of T wave
%
% Input Parameters:
%   w1: matrix with WT scales 1 to 5 for lead 1
%   w2: matrix with WT scales 1 to 5 for lead 2
%   w3: matrix with WT scales 1 to 5 for lead 1
%  intreg: interval to be considered in fitting the optimal direction
%  scale to use in delineation
%   timenew:  QRS times in inedexes refering the interval included in the current excerpt (borders excluded)
%   position0: struct vector with the detected points for 3D multilead  step1
%   position: struct vector with the detected points for 3D multilead  step2
%   position1: struct vector with the detected points for lead1 - to be replaced for multilead marks
%   samp: samples included in the current excerpt (borders excluded)
%   intervalo: numeration of the beats processed in this segment
%
%Output Parameters:
%   actualized parameters: position, position0
% A - points of fitted lines in step 1
% A2 - points of fitted lines in step 2
%  V- vectors director to the fitted line
%  V2- vectors director to the fitted line
%
% Rute Almeida: trabalho de doutoramento 2.Dec.2004
% Last update: 05AGO2011
%
% MATLAB Version 6.5.0.180913a (R13)

global recursioncount OPT
A=[];A2=[];
V=[];V2=[];
if isempty(sig)

end
% Aon=[];Von=[];
% Vg=[];
% if nargin==15
% messages = new_empty_mesg;
% end
if ~isfield(messages.setup.wavedet,'rrant_mintol')
    messages.setup.wavedet.rrant_mintol  =  0.5;    % minimum time in sec admited for current RR to be used in T wave delineation or in Exponentially averaged RR
end
if ~isfield(messages.setup.wavedet,'rrant_maxtol')
    messages.setup.wavedet.rrant_maxtol  =  1.5;    % max time in sec admited for current RR to be used in Exponentially averaged RR
end
if ~isfield(messages.setup.wavedet,'rrant_average')
    messages.setup.wavedet.rrant_average  =  [0.8 0.2];    % Exponentially averaged RR weights
end
if ~isfield(messages.setup.wavedet,'scale')
    messages.setup.wavedet.scale=4;
end
if ~isfield(messages.setup.wavedet,'scale2')
    messages.setup.wavedet.scale2=5;
end
if ~isfield(messages.setup.wavedet,'T_CSE_tol')
    messages.setup.wavedet.T_CSE_tol=30.6*2/1000; % based on CSE tolerance
end
if ~isfield(messages.setup.wavedet,'Tconvergence_crit')
    messages.setup.wavedet.Tconvergence_crit=1; % 1 maks are acepted as the same if they difer by Tconvergence_crit
end

rrant_mintol=messages.setup.wavedet.rrant_mintol;
rrant_maxtol=messages.setup.wavedet.rrant_maxtol;
rrant_average=messages.setup.wavedet.rrant_average;
scale=messages.setup.wavedet.scale;
scale2=messages.setup.wavedet.scale2;
Tconvergence_crit=messages.setup.wavedet.Tconvergence_crit;

intervalo=intervalo(1):intervalo(end);

scalei=scale*ones(1,size(intreg,1));
if ~isempty(w3)
    %scalei(sum([position1.Tscale(intervalo);position2.Tscale(intervalo);position3.Tscale(intervalo)]==5)>1)=5;
    %Rute 17Mar06
    aux=[(indexes(1,intervalo));(indexes(2,intervalo));(indexes(3,intervalo))];
    aux(1,~isnan(aux(1,:)))=position1.Tscale(aux(1,~isnan(aux(1,:))));
    aux(2,~isnan(aux(2,:)))=position2.Tscale(aux(2,~isnan(aux(2,:))));
    aux(3,~isnan(aux(3,:)))=position3.Tscale(aux(3,~isnan(aux(3,:))));
    scalei(sum(aux==5)>1)=5;
    %%%%%%
else
    scalei(sum([position1.Tscale(intervalo);position2.Tscale(intervalo)]==scale2)>1)=scale2;
end

for i=1:size(intreg,1)

    if isnan(intreg(i,1))
        position.Tscale(intervalo(i))=NaN;
        position.Toff(intervalo(i))=NaN;
        position.T(intervalo(i))=NaN;
        position.Ttipooff(intervalo(i))=NaN;
        position.Ttipo(intervalo(i))=NaN;
        position.Tprima(intervalo(i))=NaN;
        position.Ton(intervalo(i))=NaN;
        position.Ttipoon(intervalo(i))=NaN;
    else

        scale=scalei(i);
        %marks taken as the last in the 3 leads

        aux2=find(~isnan(indexes(2,intervalo(i))));
        aux3=find(~isnan(indexes(3,intervalo(i))));
        S=max([position1.S(indexes(~isnan(indexes(1,intervalo(i))),intervalo(i))) position2.S(indexes(2*aux2,intervalo(i))) position3.S(indexes(3*aux3,intervalo(i))) ]-samp(1)+1);
        QRSoff=max([position1.QRSoff(indexes(~isnan(indexes(1,intervalo(i))),intervalo(i))) position2.QRSoff(indexes(2*aux2,intervalo(i))) position3.QRSoff(indexes(3*aux3,intervalo(i)))]-samp(1)+1);

        if ~isempty(w3)
            pontos=[w1(intreg(i,1):min(intreg(i,2),length(w1)),scale) w2(intreg(i,1):min(intreg(i,2),length(w2)),scale) w3(intreg(i,1):min(intreg(i,2),length(w3)),scale) ];
        else
            pontos=[w1(intreg(i,1):min(intreg(i,2),length(w1)),scale) w2(intreg(i,1):min(intreg(i,2),length(w2)),scale)];
        end
        weight=ones(size(pontos,1),1);
        [a,v] = optimline(pontos, weight,OPT);

        A=[A;a]; %#ok<AGROW>
        V=[V;v]; %#ok<AGROW>

        %4FEB2010
        if length(intervalo)>1
            rrant=position.qrs(intervalo(2))-position.qrs(intervalo(1));
        elseif intervalo(1)>1
            rrant=position.qrs(intervalo(1))-position.qrs(intervalo(1)-1);
        else
            rrant=messages.setup.wavedet.freq/2; %JUL2011
        end

        if (i>1),
            if (rrant_maxtol*rrant>(timenew(i)-timenew(i-1)))&&(timenew(i)-timenew(i-1)>rrant_mintol*rrant),
                rrmed = rrant_average(1)*rrant + rrant_average(2)*(timenew(i)-timenew(i-1));
            else
                rrmed = rrant;  % Exponentially averaged RR
            end
        else                  % Only for the first in each segment of the ECG
            rrmed = rrant;
        end
        rrant = rrmed;        %#ok<NASGU> % For next segment

        %WT projected: from the previous QRS complex to the following QRS complex
        if ~isempty(w3)
            if i<length(timenew),
                newleadbeat=(w1(timenew(i):timenew(i+1),scale).*v(1)+w2(timenew(i):timenew(i+1),scale)*v(2)+w3(timenew(i):timenew(i+1),scale)*v(3))./norm(v);
                %signew=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2)+sig(timenew(i):timenew(i+1),3)*v(3))./norm(v);
            else                        % if last beat of the segment
                newleadbeat=(w1(timenew(i):end,scale).*v(1)+w2(timenew(i):end,scale)*v(2)+w3(timenew(i):end,scale)*v(3))./norm(v);
                %signew=(sig(timenew(i):end,1).*v(1)+sig(timenew(i):end,2)*v(2)+sig(timenew(i):end,3)*v(3))./norm(v);
            end
        else
            if i<length(timenew),
                newleadbeat=(w1(timenew(i):timenew(i+1),scale).*v(1)+w2(timenew(i):timenew(i+1),scale)*v(2))./norm(v);
                %signew=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2))./norm(v);
            else                        % if last beat of the segment
                newleadbeat=(w1(timenew(i):end,scale).*v(1)+w2(timenew(i):end,scale)*v(2))./norm(v);
                %signew=(sig(timenew(i):end,1).*v(1)+sig(timenew(i):end,2)*v(2))./norm(v);
            end
        end
        [pos,picon,picoff,janelas1,messages]=twave3D( newleadbeat,timenew,i,messages.setup.wavedet.freq,S,QRSoff,samp,rrmed,messages); %#ok<ASGLU>
%         janelas1=janelas1+timenew(i)-1;
        position0.Tscale(intervalo(i))=scale;
        position0.Ton(intervalo(i))=pos.Ton;
        position0.Toff(intervalo(i))=pos.Toff;
        position0.T(intervalo(i))=pos.T;
        position0.Ttipo(intervalo(i))=pos.Ttipo;
        position0.Tprima(intervalo(i))=pos.Tprima;

        %position0.Wavedet(intervalo(i))=newleadbeat; % derivada wavedet Maikel 30 marzo 2009
        %position0.direccion(c)=V; % vector direcci�n Maikel 30 marzo 2009
        %pontos2=[w1(picoff:intreg(i,2),scale) w2(picoff:intreg(i,2),scale) w3(picoff:intreg(i,2),scale) ];

        if ~isnan(picoff)
            if picoff>(intreg(i,2)-1)
                messages.warnings=[messages.warnings {'Tend out of the searching interval, in the fist step of multilead!'}];
                position.Tscale(intervalo(i))=NaN;
                position.Toff(intervalo(i))=NaN;
                position.T(intervalo(i))=NaN;
                position.Ttipooff(intervalo(i))=NaN;
                position.Ttipo(intervalo(i))=NaN;
                position.Tprima(intervalo(i))=NaN;
                recursioncount.Toff(intervalo(i))=NaN;
                position.contadorToff(intervalo(i))=NaN; %Rute
            else
                if (pos.Toff-samp(1)+1)==picoff
                    endaux=min([(position0.Toff(intervalo(i))-samp(1)+1+round(messages.setup.wavedet.T_CSE_tol*messages.setup.wavedet.freq)) length(w1)]); %3JUL

                else
                    endaux=min([intreg(i,2) (position0.Toff(intervalo(i))-samp(1)+1+messages.setup.wavedet.T_CSE_tol*messages.setup.wavedet.freq) ]);
                end
                if ~isempty(w3)
                    pontos2=[w1(picoff:endaux,scale) w2(picoff:endaux,scale) w3(picoff:endaux,scale)];
                    %pontos2=[w1(janelas1(2):janelas1(3),scale) w2(janelas1(2):janelas1(3),scale) w3(janelas1(2):janelas1(3),scale)];
                else
                    pontos2=[w1(picoff:endaux,scale) w2(picoff:endaux,scale)];
                end
                weight2=ones(size(pontos2,1),1);
                [a,v] = optimline(pontos2, weight2,OPT);
                A2=[A2;a]; %#ok<AGROW>
                V2=[V2;v]; %#ok<AGROW>
                if ~isempty(w3)
                    if i<length(timenew),
                        newleadbeat2=(w1(timenew(i):timenew(i+1),scale).*v(1)+w2(timenew(i):timenew(i+1),scale)*v(2)+w3(timenew(i):timenew(i+1),scale)*v(3))./norm(v);
                        %signew2=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2)+sig(timenew(i):timenew(i+1),3)*v(3))./norm(v);
                    else                        % if last beat of the segment
                        newleadbeat2=(w1(timenew(i):end,scale).*v(1)+w2(timenew(i):end,scale)*v(2)+w3(timenew(i):end,scale)*v(3))./norm(v);
                        %signew2=(sig(timenew(i):end,1).*v(1)+sig(timenew(i):end,2)*v(2)+sig(timenew(i):end,3)*v(3))./norm(v);
                    end
                else
                    if i<length(timenew),
                        newleadbeat2=(w1(timenew(i):timenew(i+1),scale).*v(1)+w2(timenew(i):timenew(i+1),scale)*v(2))./norm(v);
                        %signew2=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2))./norm(v);
                    else                        % if last beat of the segment
                        newleadbeat2=(w1(timenew(i):end,scale).*v(1)+w2(timenew(i):end,scale)*v(2))./norm(v);
                        %signew2=(sig(timenew(i):end,1).*v(1)+sig(timenew(i):end,2)*v(2))./norm(v);
                    end
                end
                [pos,picon2,picoff2,janelas2,messages]=twave3D(newleadbeat2,timenew,i,messages.setup.wavedet.freq,S,QRSoff,samp,rrmed,messages); %#ok<ASGLU>
%                 janelasnew=janelas2+timenew(i)-1;
                picoffall=[picoff picoff2];
                if isnan(picoff2)
                    amppicoffall=[abs(newleadbeat(picoff-timenew(i)+1)) NaN];
                else
                    amppicoffall=[abs(newleadbeat(picoff-timenew(i)+1)) abs(newleadbeat2(picoff2-timenew(i)+1))];
                end
                newToff=[position0.Toff(intervalo(i)) pos.Toff];

                if (amppicoffall(end)<amppicoffall(end-1))
                    position.contadorToff(intervalo(i))=-1;
                    %Vg(i,:)=V(end,:);
                    position.Tscale(intervalo(i))=position0.Tscale(intervalo(i));
                    position.Toff(intervalo(i))=position0.Toff(intervalo(i));
                    position.T(intervalo(i))= position0.T(intervalo(i));
                    position.Ttipooff(intervalo(i))=position0.Ttipo(intervalo(i));
                    position.Ttipo(intervalo(i))=NaN;
                    position.Tprima(intervalo(i))=position0.Tprima(intervalo(i));
                    %position.contadorToff=recursioncount.Toff;
                    position.direccionToffopt(:,intervalo(i))=V(end,:);
                    %position.Wavedet(intervalo(i))=newleadbeat2; % derivada wavedet Maikel 30 marzo 2009
                    %position.direccion(:,intervalo(i))=V2; % vector direcci�n Maikel 30 marzo 2009
                else
                    position.contadorToff(intervalo(i))=0;
                    %Vg(i,:)=V2(end,:);
                    position.direccionToffopt(:,intervalo(i))=V2(end,:);
                    %%%%%%%%%%%%%%%%%%%%%%recursion
                    while (~isnan(picoffall(end)) &&  abs(newToff(end)-newToff(end-1))>Tconvergence_crit) && (amppicoffall(end)>amppicoffall(end-1))
                        position.contadorToff(intervalo(i))=position.contadorToff(intervalo(i))+1;
                        %maxTintervalend=maxTintervalend; % teste 14
                        if (pos.Toff-samp(1)+1)==picoff
                            endaux=(pos.Toff-samp(1)+1+round(messages.setup.wavedet.T_CSE_tol*messages.setup.wavedet.freq));
                        else
                            endaux=min([intreg(i,2) (pos.Toff-samp(1)+1+round(messages.setup.wavedet.T_CSE_tol*messages.setup.wavedet.freq))]);
                        end
                        if ~isempty(w3)
                            pontosnew=[w1(picoffall(end):endaux,scale) w2(picoffall(end):endaux,scale) w3(picoffall(end):endaux,scale)];
                            %pontosnew=[w1(picoff(end):endaux,scale) w2(picoff(end):endaux,scale) w3(picoff(end):endaux,scale)];
                            %pontosnew=[w1(janelasnew(2):janelasnew(3),scale) w2(janelasnew(2):janelasnew(3),scale) w3(janelasnew(2):janelasnew(3),scale)];
                        else
                            pontosnew=[w1(picoffall(end):endaux,scale) w2(picoffall(end):endaux,scale)];
                        end

                        if   size(pontosnew,1)>1
                            weightnew=ones(size(pontosnew,1),1);
                            [a,v] = optimline(pontosnew, weightnew,OPT);
%                             if i==6
%                                 title([ecgnr '(Tend-step' num2str(2+position.contadorToff(intervalo(i)))')'])
%                             else
%                                 close
%                             end
                            A2=[A2;a]; %#ok<AGROW>
                            V2=[V2;v]; %#ok<AGROW>
                            if ~isempty(w3)
                                if i<length(timenew),
                                    newleadbeatnew=(w1(timenew(i):timenew(i+1),scale).*v(1)+w2(timenew(i):timenew(i+1),scale)*v(2)+w3(timenew(i):timenew(i+1),scale)*v(3))./norm(v);
                                    %signewnew=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2)+sig(timenew(i):timenew(i+1),3)*v(3))./norm(v);
                                else                        % if last beat of the segment
                                    newleadbeatnew=(w1(timenew(i):end,scale).*v(1)+w2(timenew(i):end,scale)*v(2)+w3(timenew(i):end,scale)*v(3))./norm(v);
                                    %signewnew=(sig(timenew(i):end,1).*v(1)+sig(timenew(i):end,2)*v(2)+sig(timenew(i):end,3)*v(3))./norm(v);
                                end
                            else
                                if i<length(timenew),
                                    newleadbeatnew=(w1(timenew(i):timenew(i+1),scale).*v(1)+w2(timenew(i):timenew(i+1),scale)*v(2))./norm(v);
                                    %signewnew=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2))./norm(v);
                                else                        % if last beat of the segment
                                    newleadbeatnew=(w1(timenew(i):end,scale).*v(1)+w2(timenew(i):end,scale)*v(2))./norm(v);
                                    %signewnew=(sig(timenew(i):end,1).*v(1)+sig(timenew(i):end,2)*v(2))./norm(v);
                                end
                            end

                            [pos,picon2,picoffnew,janelasnew,messages]=twave3D(newleadbeatnew,timenew,i,messages.setup.wavedet.freq,S,QRSoff,samp,rrmed,messages); %#ok<ASGLU>
%                             janelasnew=janelas2+timenew(i)-1;
                        else
                            picoffnew=NaN;
                        end
                        picoffall=[picoffall picoffnew]; %#ok<AGROW>
                        if isnan(picoffnew)
                            amppicoffall=[amppicoffall NaN]; %#ok<AGROW>
                        else
                            amppicoffall=[amppicoffall abs(newleadbeatnew(picoffnew-timenew(i)+1))]; %#ok<AGROW>
                        end
                        newToff=[newToff pos.Toff]; %#ok<AGROW>

                        if (amppicoffall(end)<amppicoffall(end-1))
                            pos.Toff=NaN; % because the lead got worse
                        end
                        if isnan(pos.Toff); %01Jun05
                            pos.Toff=newToff(end-1);
                            position.contadorToff(intervalo(i))=position.contadorToff(intervalo(i))-1;
                        else
                            %Vg(i,:)=V2(end,:);
                            position.direccionToffopt(:,intervalo(i))=V2(end,:);
                            %contador=[contador; recursioncount.Toff];
                        end
                        if  prod(picoffall(1:end-1)-picoffall(end))==0  % ciclo infinito no teste 14
                            picoffall(end)=NaN;
                        end
                    end
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    position.Tscale(intervalo(i))=scale;
                    position.Toff(intervalo(i))=pos.Toff;
                    position.T(intervalo(i))=pos.T;
                    position.Ttipooff(intervalo(i))=pos.Ttipo;
                    position.Ttipo(intervalo(i))=NaN;
                    position.Tprima(intervalo(i))=pos.Tprima;
                    %position.contadorToff=recursioncount.Toff;
                    %position.Wavedet(intervalo(i))=newleadbeatnew; % derivada wavedet Maikel 30 marzo 2009
                    %position.direccion(:,intervalo(i))=V2; % vector direcci�n Maikel 30 marzo 2009

                end
                %
            end
        else
            position.Tscale(intervalo(i))=NaN;
            position.Toff(intervalo(i))=NaN;
            position.T(intervalo(i))=NaN;
            position.Ttipooff(intervalo(i))=NaN;
            position.Ttipo(intervalo(i))=NaN;
            position.Tprima(intervalo(i))=NaN;
            recursioncount.Toff(intervalo(i))=NaN;
            position.contadorToff(intervalo(i))=NaN; %Rute
            %position.Wavedet(intervalo(i))=NaN;  % derivada wavedet Maikel 30 marzo 2009
            %position.direccion(:,intervalo(i))=NaN; % vector direcci�n Maikel 30 marzo 2009
        end

        Tbeginauxiliar
    end
end
% position0.A(intervalo(i))=A;
% position0.A2(intervalo(i))=A2;
% position0.V(intervalo(i))=V;
% position0.V2(intervalo(i))=V2;

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit