text:
function [signal,H] = sload(FILENAME,varargin)
% SLOAD loads signal data of various data formats
% [signal,header] = sload(FILENAME)
% [signal,header] = sload(FILENAME,CHAN)
%       read selected channels in list CHAN
%       CHAN=0 [default], reads all channels
%
% [signal,header] = sload(dir('f*.emg') ... )
% [signal,header] = sload('f*.emg' ...)
%  	loads all files 'f*.emg'
%
% [signal,header] = sload(..., PropertyName1,PropertyValue1,...)
%  	PropertyName(s)		PropertyValue
%	'UCAL'			-		data uncalibrated (not scaled)
%	'OVERFLOWDETECTION'	'On'		[default]
%				'Off'		no overflow detection
%	'OUTPUT'		'single'	single precision data [default: 'double']
%	'NUMBER_OF_NAN_IN_BREAK'   N		inserts N NaN's between two concatanated segments
%						default: N=100
%	'SampleRate'		Fs		target sampling rate (supports resampling)
%	'EOG_CORRECTION'	'On'		uses two-channel regression analysis - if possible
%				'Off'		no correction of EOG artifacts [default]
%       'CNT32', '32bit'			force CNT 32bit format
%
% The list of supported formats is available here:
% http://pub.ist.ac.at/~schloegl/biosig/TESTED
%
%    SLOAD loads all the data (of the selected channels)
%    at once. In case of large data files, This can be
%    a problem. Instead, You can use
%       HDR = sopen(FILENAME,'r');
%       [s,HDR]=sread(HDR,duration_segment1);
%       [s,HDR]=sread(HDR,duration_segment2);
%       ....
%       [s,HDR]=sread(HDR,duration_segmentM);
%       HDR = sclose(HDR);
%
% see also: SVIEW, SOPEN, SREAD, SCLOSE, SAVE2BKR, TLOAD
%
% In order to increase the speed, install mexSLOAD.mex from biosig4c++
%
% Reference(s):

%	$Id: sload.m 2656 2011-03-18 22:30:26Z schloegl $
%	Copyright (C) 1997-2007,2008,2009,2010,2011 by Alois Schloegl
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

%
%    BioSig 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.
%
%    BioSig 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 BioSig.  If not, see <http://www.gnu.org/licenses/>.

if length(varargin)<2;
	MODE = '';
else
	MODE = varargin{2};
end;

CHAN = 0;
STATE.CNT32 = 0;
STATE.UCAL = 0;
STATE.OVERFLOWDETECTION = 1;
STATE.NUMBER_OF_NAN_IN_BREAK = 100;
STATE.EOG_CORRECTION = 0 ;
STATE.OUTPUT = 'double';

Fs = NaN;
k = 1;
while (k<=length(varargin))
	if isnumeric(varargin{k})
		if (k==1), CHAN=varargin{k}; end;
	elseif ~isempty(strfind(varargin{k},'UCAL'))
		MODE = varargin{k};
		if strcmpi(varargin{k+1},'on');
			STATE.UCAL = 1;
			k=k+1;
		elseif strcmpi(varargin{k+1},'off');
			STATE.UCAL = 0;
			k=k+1;
		else
			STATE.UCAL = 1;
		end;
	elseif ~isempty(strfind(varargin{k},'OVERFLOWDETECTION:OFF'))
		MODE = varargin{k};
		STATE.OVERFLOWDETECTION = 0;
	elseif strcmpi(varargin{k},'OVERFLOWDETECTION')
		if strcmpi(varargin{k+1},'on');
			STATE.OVERFLOWDETECTION = 1;
		elseif strcmpi(varargin{k+1},'off');
			STATE.OVERFLOWDETECTION = 0;
		elseif isnumeric(varargin{k+1})
			STATE.OVERFLOWDETECTION = varargin{k+1};
		end;
		k=k+1;
	elseif ~isempty(strfind(varargin{k},'OUTPUT:SINGLE'))
		MODE = varargin{k};
		STATE.OUTPUT = 'single';
	elseif strcmpi(varargin{k},'OUTPUT')
		STATE.OUTPUT = varargin{k+1};
	elseif strcmpi(varargin{k},'EOG_CORRECTION')
		if strcmpi(varargin{k+1},'on');
			STATE.EOG_CORRECTION = 1;
		elseif strcmpi(varargin{k+1},'off');
			STATE.EOG_CORRECTION = 0;
		elseif isnumeric(varargin{k+1})
			STATE.EOG_CORRECTION = varargin{k+1};
		end;
		k=k+1;
	elseif strcmpi(varargin{k},'Number_of_nan_in_break')
		STATE.NUMBER_OF_NAN_IN_BREAK = varargin{k+1};
		k=k+1;
	elseif strcmpi(varargin{k},'SampleRate')
		Fs = varargin{k+1};
		k=k+1;
	elseif strcmpi(varargin{k},'CNT32') || strcmpi(varargin{k},'32bit')
		STATE.CNT32 = 1;
	elseif exist(varargin{k},'file') && (k==1)
	        CHAN = varargin{k};
	end;
	k=k+1;
end;

%%% resolve wildcards %%%
if ischar(FILENAME)
if any(FILENAME=='*')
        p = fileparts(FILENAME);
        f = dir(FILENAME);
        EOGix = zeros(1,length(f));
        for k = 1:length(f);
                f(k).name = fullfile(p,f(k).name);
                [p,g,e]   = fileparts(f(k).name);
                lg = length(g);
                if (lg>2) && strcmp(upper(g(lg+(-2:0))),'EOG')
                        EOGix(k) = 1;
                end
        end;
        FILENAME={f([find(EOGix),find(~EOGix)]).name};
end;
end;

if nargout>0,
	signal = [];
end;

if (iscell(FILENAME) && (length(FILENAME)==1)),
	FILENAME = FILENAME{k};
end;

if ((iscell(FILENAME) || isstruct(FILENAME))),
	signal = [];
	for k = 1:length(FILENAME),
		if iscell(FILENAME(k))
			f = FILENAME{k};
		else
			f = FILENAME(k);
		end

                [s,h] = sload(f,varargin{:});
		if k==1,
			H = h;
			signal = s;
			H.SegLen = [0,size(s,1)];
                        %H.EVENT.POS = [H.EVENT.POS; 0];
                        %H.EVENT.TYP = [H.EVENT.TYP; hex2dec('7ffe')];
                        %if isfield(H.EVENT,'CHN');
                        %        H.EVENT.CHN = [H.EVENT.CHN; 0];
                        %end;
                        %if isfield(H.EVENT,'DUR');
                        %        H.EVENT.DUR = [H.EVENT.DUR; size(s,1)];
                        %end;
                        %if isfield(H.EVENT,'Desc');	% TFM-Excel-Beat-to-Beat
                        %        H.EVENT.Desc = [H.EVENT.Desc; {'New Segment'}; h.EVENT.Desc];
                        %end;
		else
			H.FILE(k) = h.FILE;
                        H.T0(k,1:6) = h.T0;
			if ~isnan(h.SampleRate) && (H.SampleRate ~= h.SampleRate),
				fprintf(2,'Warning SLOAD: sampling rates of multiple files differ %i!=%i.\n',H.SampleRate, h.SampleRate);
			end;
                        if size(s,2)==size(signal,2), %(H.NS == h.NS)
				signal = [signal; repmat(NaN,STATE.NUMBER_OF_NAN_IN_BREAK,size(s,2)); s];
				H.SegLen = [H.SegLen,size(signal,1)];
			else
				error('ERROR SLOAD: incompatible channel numbers %i!=%i of multiple files\n',H.NS,h.NS);
			end;

                        if ~isempty(h.EVENT.POS),
                                H.EVENT.POS = [H.EVENT.POS; size(signal,1); h.EVENT.POS+size(signal,1)-size(s,1)];
                                H.EVENT.TYP = [H.EVENT.TYP; hex2dec('7ffe'); h.EVENT.TYP];
                                if isfield(H.EVENT,'CHN');
        	                	H.EVENT.CHN = [H.EVENT.CHN; 0; h.EVENT.CHN];
                                end;
                                if isfield(H.EVENT,'DUR');
      					H.EVENT.DUR = [H.EVENT.DUR; size(s,1); h.EVENT.DUR];
                                end;
                                if isfield(H.EVENT,'Desc');	% TFM-Excel-Beat-to-Beat
                                        H.EVENT.Desc = [H.EVENT.Desc; {'New Segment'}; h.EVENT.Desc];
                                end;
                        end;

                        if ~isequal(H.Label,h.Label) || ~isequal(H.PhysDimCode,h.PhysDimCode)
				warning('Labels and PhysDim of multiple files differ!\n');
                                for k2 = 1:length(H.InChanSelect),
                                	k1 = H.InChanSelect(k2);
                                        if ~strcmp(H.Label{k1},h.Label{k1}) || (H.PhysDimCode(k1)~=h.PhysDimCode(k1)),
                                                warning('#%02i:  %s | %s | (%i)%s | (%i)%s\n',k1, H.Label{k1},h.Label{k1}, H.PhysDimCode(k1),H.PhysDim{k1},h.PhysDimCode(k1),h.PhysDim{k1});
                                        end;
                                end;
                        end;
                        if ~isequal(H.Calib,h.Calib) && H.FLAG.UCAL,
				error('SLOAD: concatanating uncalibrated data with different scaling factors does not make sense!');
                        end;

                        if isfield(H,'TriggerOffset'),
                                if H.TriggerOffset ~= h.TriggerOffset,
                                        fprintf(2,'Warning SLOAD: Triggeroffset %f does not fit.%f \n',H.TriggerOffset,h.TriggerOffset);
                                end;
                        end;

                        if isfield(h,'TRIG'),
                                if ~isfield(H,'TRIG'),
                                        H.TRIG = [];
                                end;
                                H.TRIG = [H.TRIG(:); h.TRIG(:)+size(signal,1)-size(s,1)];
                        else h.TRIG=[];
                        end;

                        if isfield(h,'Classlabel'),
				if isfield(h,'ArtifactSelection'),
                                	if (any(h.ArtifactSelection>1) || (length(h.ArtifactSelection) < length(h.Classlabel)))
                                        	sel = zeros(size(h.Classlabel));
                                                sel(h.ArtifactSelection) = 1;
                                        else
                                                sel = h.ArtifactSelection(:);
                                        end;
                                else
                                	h.ArtifactSelection = repmat(logical(0),length(h.Classlabel),1);
                                end;

                                if isfield(H,'ArtifactSelection')
                                	H.ArtifactSelection = [H.ArtifactSelection; h.ArtifactSelection(:)];
                                else
                                	H.ArtifactSelection = [repmat(logical(0),length(H.Classlabel),1); h.ArtifactSelection(:)];
                                end;
                                H.Classlabel = [H.Classlabel(:); h.Classlabel(:)];

                        elseif isfield(H,'Classlabel') && (length(h.TRIG)>0),
			        warning(sprintf('field classlabel missing in %s. This might corrupt the resulting classlabels!',h.FileName));

                        end;
			clear s;
                end;
	end;

	% fprintf(1,'SLOAD: data segments are concatenated with NaNs in between.\n');
	return;
end;
%%% end of multi-file section

%%%% start of single file section
%%%%%%%%%% --------- EOG CORRECTION -------------- %%%%%%
if STATE.EOG_CORRECTION,
try
        h = get_regress_eog(FILENAME,'REG fft16');
catch,
	fprintf(H.FILE.stderr,'Error: SLOAD (EOG_CORRECTION): %s\n',lasterr);
	H.FLAG.EOG_CORRECTION = 0;
end;
end;

FlagLoaded = 0;
if exist('mexSLOAD','file')==3,
	try
		valid_rerefmx = 1;
		if ischar(CHAN)
		        HDR = sopen(CHAN,'r'); HDR=sclose(HDR);
		        if isfield(HDR.Calib),
		                ReRefMx = HDR.Calib;
		        else
        			valid_rerefmx = 0;
		        end;
		elseif all(size(CHAN)>1) || any(floor(CHAN)~=CHAN) || any(CHAN<0) || (any(CHAN==0) && (numel(CHAN)>1));
		        ReRefMx = CHAN;
	        	CHAN = find(any(CHAN,2));
		elseif all(CHAN>0) && all(floor(CHAN)==CHAN),
			[tmp,ix]= sort(CHAN);
		        ReRefMx = sparse(CHAN,1:length(CHAN),1);
		else
		        ReRefMx = [];
			valid_rerefmx=0;
		end

		if STATE.EOG_CORRECTION,
			ReRefMx = h.r0*ReRefMx;
			valid_rerefmx=1;
		end
		if STATE.OVERFLOWDETECTION,
			arg1 = 'OVERFLOWDETECTION:ON';
		else
			arg1 = 'OVERFLOWDETECTION:OFF';
		end
		if STATE.UCAL,
			arg2 = {'UCAL:ON'};
		else
			arg2 = {'UCAL:OFF'};
		end
		if STATE.CNT32
			arg2 = {arg2,'CNT32'};
		end
		if ~valid_rerefmx,
			[signal,HDR] = mexSLOAD(FILENAME,0,arg1,arg2{:});
			if isfield(HDR,'ErrNum') && (HDR.ErrNum>0),
				fprintf(1,'%s\n',HDR.ErrMsg);
				H = HDR;
				x = HDR.nonsense; 	% hack: triggers CATCH branch in TRY-CATCH-END;
			end

			if isfield(HDR.FLAG,'ROW_BASED_CHANNELS') && HDR.FLAG.ROW_BASED_CHANNELS, signal = signal.'; end;
			FlagLoaded   = isfield(HDR,'NS');
			HDR.InChanSelect = 1:HDR.NS;
		else
			InChanSelect = find(any(ReRefMx,2));
			[signal,HDR] = mexSLOAD(FILENAME,InChanSelect,arg1,arg2{2});
			if isfield(HDR,'ErrNum') && HDR.ErrNum==3,
				%% file not found - fopen failed
				H = HDR;
				return;
			end
			if isfield(HDR.FLAG,'ROW_BASED_CHANNELS') && HDR.FLAG.ROW_BASED_CHANNELS, signal = signal.'; end;
			FlagLoaded   = isfield(HDR,'NS');
			HDR.InChanSelect = InChanSelect(InChanSelect <= HDR.NS);
			signal = signal*ReRefMx(HDR.InChanSelect,:); %% can be sparse if just a single channel is loaded
			signal = full(signal);  %% make sure signal is not sparse
		end;

		HDR.T0 = datevec(HDR.T0);
		HDR.Patient.Birthday = datevec(HDR.Patient.Birthday);
		HDR.Calib = [HDR.Off(:)';diag(HDR.Cal)];
		[HDR.FILE.Path,HDR.FILE.Name,HDR.FILE.Ext] = fileparts(FILENAME);
		HDR.FileName = FILENAME;
		HDR = leadidcodexyz(HDR);
		HDR.EVENT.POS = HDR.EVENT.POS+1; % convert from 0-based to 1-based index

		H=HDR;
		H.FLAG.EOG_CORRECTION = STATE.EOG_CORRECTION;

		if isfield(HDR,'Patient') && isfield(HDR.Patient,'Weight') && isfield(HDR.Patient,'Height')
			%% Body Mass Index
			HDR.Patient.BMI = HDR.Patient.Weight * HDR.Patient.Height^-2 * 1e4;

		       	%% Body Surface Area
			% DuBois D, DuBois EF. A formula to estimate the approximate surface area if height and weight be known. Arch Intern Medicine. 1916; 17:863-71.
			% Wang Y, Moss J, Thisted R. Predictors of body surface area. J Clin Anesth. 1992; 4(1):4-10.
		       	HDR.Patient.BSA = 0.007184 * HDR.Patient.Weight^0.425 * HDR.Patient.Height^0.725;
		end;

		if strcmp(H.TYPE,'GDF');
                        % Classlabels according to
                        % http://biosig.cvs.sourceforge.net/*checkout*/biosig/biosig/doc/eventcodes.txt

                        % sort event table before extracting HDR.Classlabel and HDR.TRIG
			[H.EVENT.POS,ix] = sort(H.EVENT.POS);
			H.EVENT.TYP = H.EVENT.TYP(ix);
			if isfield(H.EVENT,'CHN')
				H.EVENT.CHN = H.EVENT.CHN(ix);
			end;
			if isfield(H.EVENT,'DUR')
				H.EVENT.DUR = H.EVENT.DUR(ix);
			end;

			if (length(H.EVENT.TYP)>0)
                                ix = (H.EVENT.TYP>hex2dec('0300')) & (H.EVENT.TYP<hex2dec('030d'));
                                ix = ix | ((H.EVENT.TYP>=hex2dec('0320')) & (H.EVENT.TYP<=hex2dec('037f')));
                                ix = ix | (H.EVENT.TYP==hex2dec('030f')); % unknown/undefined cue
                                H.Classlabel = mod(H.EVENT.TYP(ix),256);
                                H.Classlabel(H.Classlabel==15) = NaN; % unknown/undefined cue
                        end;

                        % Trigger information and Artifact Selection
                        ix = find(H.EVENT.TYP==hex2dec('0300'));
                        H.TRIG = H.EVENT.POS(ix);
                        ArtifactSelection = repmat(logical(0),length(ix),1);
                        for k = 1:length(ix),
                                ix2 = find(H.EVENT.POS(ix(k))==H.EVENT.POS);
                                if any(H.EVENT.TYP(ix2)==hex2dec('03ff'))
                                        ArtifactSelection(k) = logical(1);
                                end;
                        end;
                        if any(ArtifactSelection), % define only if necessary
                                H.ArtifactSelection = ArtifactSelection;
                        end;

			% apply channel selections to EVENT table
			if valid_rerefmx && ~isempty(H.EVENT.POS) && (isfield(H.EVENT,'CHN')),	% only if channels are selected.
				sel = (H.EVENT.CHN(:)==0);	% memory allocation, select all general events
				for k = find(~sel'),		% select channel specific elements
					sel(k) = any(H.EVENT.CHN(k)==InChanSelect);
				end;
				H.EVENT.POS = H.EVENT.POS(sel);
				H.EVENT.TYP = H.EVENT.TYP(sel);
				H.EVENT.DUR = H.EVENT.DUR(sel);	% if EVENT.CHN available, also EVENT.DUR is defined.
				H.EVENT.CHN = H.EVENT.CHN(sel);
				% assigning new channel number
				a = zeros(1,HDR.NS);
				for k = 1:length(InChanSelect),		% select channel specific elements
					a(InChanSelect(k)) = k;		% assigning to new channel number.
				end;
				ix = H.EVENT.CHN>0;
				H.EVENT.CHN(ix) = a(H.EVENT.CHN(ix));	% assigning new channel number
			end;

		elseif strcmp(H.TYPE,'BDF');
			H.BDF.ANNONS = zeros(HDR.NRec*HDR.SPR,1);
			ix  = HDR.EVENT.TYP==hex2dec('7ffe');
			ix1 = find(~ix);
			POS = HDR.EVENT.POS(ix1);
			TYP = HDR.EVENT.TYP(ix1);
			[s,ix2] = sort([POS; HDR.NRec*HDR.SPR+1]);
			for k = 1:length(s)-1,
				H.BDF.ANNONS(s(k):s(k+1)-1) = TYP(ix2(k));
			end;
			H.BDF.ANNONS(HDR.EVENT.POS(ix)) = H.BDF.ANNONS(HDR.EVENT.POS(ix)) | hex2dec('10000');

		elseif strcmp(H.TYPE,'BKR');
	                H.Classlabel = [];
	                H.TRIG = [];
		        tmp=fullfile(H.FILE.Path,[H.FILE.Name,'.mat']);
        		if ~exist(tmp,'file'),
        		        tmp=fullfile(H.FILE.Path,[H.FILE.Name,'.MAT']);
       			end
		        x = [];
        		if exist(tmp,'file'),
                		x = load('-mat',tmp);
        		end;

                        if isfield(x,'header'),
                                H.MAT  = x.header;
                                if isfield(x.header,'Setup'),
                                        if isfield(x.header.Setup,'Bits'),
                                                H.Bits = x.header.Setup.Bits;
                                                [datatyp, limits, datatypes] = gdfdatatype(H.Bits+255);
                                                % THRESHOLD for Overflow detection
                                                if ~isfield(H,'THRESHOLD')
                                                    H.THRESHOLD = repmat(limits, H.NS, 1);
                                                end;
                                        end;
                                end;
                                if isfield(x.header,'Result') && isfield(x.header.Result,'Classlabel'),
                                        H.Classlabel = x.header.Result.Classlabel;
                                end;
                                if isfield(x.header,'Paradigm')
                                        if isempty(H.Classlabel) && isfield(x.header.Paradigm,'Classlabel')
                                                H.Classlabel = x.header.Paradigm.Classlabel;
                                        end;
                                        H.BCI.Paradigm = x.header.Paradigm;
                                        if isfield(H.BCI.Paradigm,'TriggerOnset');
                                                H.TriggerOffset = H.BCI.Paradigm.TriggerOnset;
                                        elseif isfield(H.BCI.Paradigm,'TriggerTiming');
                                            %    H.BCI.Paradigm.TriggerTiming,
                                                H.TriggerOffset = H.BCI.Paradigm.TriggerTiming;
                                                fprintf(2,'Warning BKROPEN: Paradigm.TriggerOnset is unknown. Paradigm.TriggerTiming= %f ms is used instead\n',H.TriggerOffset);
                                        end;
                                end;

                                if isfield(x.header,'PhysioRec'), % R. Leeb's data
                                       H.Label = cellstr(x.header.PhysioRec);
                                end;
                                if isfield(x.header,'BKRHeader'), % R. Scherer Data
                                        if isfield(x.header.BKRHeader,'TO'),
                                                H.T0 = x.header.BKRHeader.TO;
                                        end;
                                        if isfield(x.header.BKRHeader,'Label'),
                                                H.Label = cellstr(x.header.BKRHeader.Label);
                                                ns = H.NS-length(H.Label);
                                                if ns == 1;
                                                        H.Label = strvcat(H.Label,'TRIGGER');
                                                elseif ns > 1;
                                                        H.Label = strvcat(H.Label,char(repmat('n.a.',ns,1)));
                                                end;
                                        end;
                                end;
                                if isfield(x.header,'Model'), % More
                                        if isfield(x.header.Model,'AnalogInput'),
                                                for k = 1:length(x.header.Model.AnalogInput),
                                                        H.Filter.HighPass(k) = x.header.Model.AnalogInput{k}{5};
                                                        H.Filter.LowPass(k)  = x.header.Model.AnalogInput{k}{6};
                                                        H.Filter.Notch(k)    = strcmpi(x.header.Model.AnalogInput{k}{7},'on');

                                                        H.MAT.Cal(k) = x.header.Model.AnalogInput{k}{3};
                                                end
                                        end;
                                end;
                                if ~isempty(strmatch('TRIGGER',H.Label))
                                        H.AS.TRIGCHAN = H.NS; %strmatch('TRIGGER',H.Label);
                                end;
                        end;
                        if 1; ~isfield(H,'Classlabel');
                                tmp=fullfile(H.FILE.Path,[H.FILE.Name,'.par']);
                                if ~exist(tmp,'file'),
                                    tmp=fullfile(H.FILE.Path,[H.FILE.Name,'.PAR']);
                                end
                                if exist(tmp,'file'),
                                        H.Classlabel = load(tmp);
                                end;
                        end;

                        %%% Artifact Selection files
                        tmp1=fullfile(H.FILE.Path,[H.FILE.Name,'.sel']);
                        if ~exist(tmp1,'file'),
                                tmp1=fullfile(H.FILE.Path,[H.FILE.Name,'.SEL']);
                        end
                        tmp2 = fullfile(H.FILE.Path,[H.FILE.Name,'_artifact.mat']);
                        SW   = (exist(tmp1,'file')>0) + 2*(exist(tmp2,'file')>0);
                        if SW == 0,
                        elseif SW == 1,
                                if exist('OCTAVE_VERSION','builtin')
                                        H.ArtifactSelection = load('-ascii',tmp1);
                                else
                                        H.ArtifactSelection = load(tmp1);
                                end;
                        elseif SW == 2,
                                if exist('OCTAVE_VERSION','builtin')
                                        tmp = load('-mat',tmp2);
                                else
                                        tmp = load(tmp2);
                                end;
                                H.ArtifactSelection = tmp.artifact(:);
                        elseif SW == 3,
                                fprintf(H.FILE.stderr,'Warning BKROPEN: more than one ArtifactSelection files. File %s is used.\n',tmp1);
                                if exist('OCTAVE_VERSION')>5
                                        H.ArtifactSelection = load('-ascii',tmp1);
                                else
                                        H.ArtifactSelection = load(tmp1);
                                end;
                        end;
                        if isfield(H,'ArtifactSelection'),
                                if any(H.ArtifactSelection>1) || (length(H.ArtifactSelection)<length(H.Classlabel))
                                        sel = zeros(size(H.Classlabel));
                                        sel(H.ArtifactSelection) = 1;
                                        H.ArtifactSelection = sel(:);
                                end;
                                H.ArtifactSelection = H.ArtifactSelection(:);
                        end;

                        if isfield(H.AS,'TRIGCHAN') % & isempty(H.EVENT.POS)
                                if H.AS.TRIGCHAN <= H.NS, %size(H.data,2),
                                        H.THRESHOLD(H.AS.TRIGCHAN,1:2) = [-1-2^15,2^15]; % do not apply overflow detection for Trigger channel
                                        data = mexSLOAD(H.FileName,H.AS.TRIGCHAN,'UCAL:ON','OVERFLOWDETECTION:OFF');
                                        TRIGon = gettrigger(data);
                                        %TRIGoff = gettrigger(-double(H.data(:,H.AS.TRIGCHAN)));
                                        if isfield(H,'TriggerOffset')
                                                TRIGon  = TRIGon - round(H.TriggerOffset/1000*H.SampleRate);
                                        %        TRIGoff = TRIGoff - round(H.TriggerOffset/1000*H.SampleRate);
                                        end;
                                end;
                                H.TRIG = TRIGon(:);
                                H.EVENT.POS = TRIGon(:); %[TRIGon(:); TRIGoff(:)];
                                H.EVENT.TYP = repmat(hex2dec('0300'),numel(TRIGon),1); %repmat(hex2dec('8300'),numel(TRIGoff),1)];
                        end;
                        if length(H.TRIG)~=length(H.Classlabel),
                                % hack to deal with BCI22 data
                                fprintf(2,'Warning BKROPEN: Number of triggers (%i) and number of Classlabels (%i) do not fit\n',length(H.TRIG),length(H.Classlabel));
                                H.TRIG = [];
                                H.Classlabel = [];
                                H.ArtifactSelection = [];
                        end;
			%% end of BKR

		elseif strcmp(H.TYPE,'BrainVision');
        	        try
        	        	H = bv2biosig_events(H);
       	        	catch
       	        		%warning('bv2biosig_events not executed');
       	        	end;

		elseif strcmp(H.TYPE,'EDF') && isfield(H.EVENT,'Desc');
			if length(H.EVENT.TYP)==length(H.EVENT.Desc)
	                        [H.EVENT.CodeDesc, CodeIndex, H.EVENT.TYP] = unique(H.EVENT.Desc);
	                end;

			%% end of BKR

		elseif strcmp(H.TYPE,'BrainVision');
			HDR = sopen(fullfile(H.FILE.Path,[H.FILE.Name,'.vmrk']));
			HDR = sclose(HDR);
			H.EVENT = HDR.EVENT;

		elseif strcmp(H.TYPE,'PDP');
			signal = repmat(NaN,max(HDR.EVENT.POS),HDR.NS);
                end;

		if isempty(signal);
			signal = repmat(NaN,round(max(HDR.EVENT.POS)*HDR.SampleRate/HDR.EVENT.SampleRate),HDR.NS);
			for k = 1:HDR.NS,
				ix = find(HDR.EVENT.CHN==k);
				signal(round(HDR.EVENT.POS(ix)*HDR.SampleRate/HDR.EVENT.SampleRate),k)=HDR.EVENT.DUR(ix);
			end;
                end;

	        H.CHANTYP = repmat(' ',1,H.NS);
		for k=1:H.NS,
			if     ~isempty(strfind(lower(H.Label{k}),'eeg')) 	H.CHANTYP(k) = 'E';
			elseif ~isempty(strfind(lower(H.Label{k}),'meg')) 	H.CHANTYP(k) = 'E';
			elseif ~isempty(strfind(lower(H.Label{k}),'emg')) 	H.CHANTYP(k) = 'M';
			elseif ~isempty(strfind(lower(H.Label{k}),'eog')) 	H.CHANTYP(k) = 'O';
			elseif ~isempty(strfind(lower(H.Label{k}),'ecg')) 	H.CHANTYP(k) = 'C';
			elseif ~isempty(strfind(lower(H.Label{k}),'air')) 	H.CHANTYP(k) = 'R';
			elseif ~isempty(strfind(lower(H.Label{k}),'trig')) 	H.CHANTYP(k) = 'T';
			end;
		end;

	catch
		%fprintf(1,lasterr);
		fprintf(1, 'SLOAD: mexSLOAD failed - the slower M-function is used.\n');
	end;
else
	global FLAG_HINT_mexSLOAD;
	if isempty(FLAG_HINT_mexSLOAD)
		fprintf(1, 'Hint: the performance of SLOAD can be improved with mexSLOAD.mex which is part of biosig4c++.\n');
		FLAG_HINT_mexSLOAD = 1;	 	% turn off hint
	end;
end;

if ~FlagLoaded,
	H = getfiletype(FILENAME);

	if isempty(H)
		fprintf(2,'Warning SLOAD: no file found\n');
		return;
	else
		% FILENAME can be fn.name struct, or HDR struct.
		FILENAME = H.FileName;
	end;
	H.FLAG.UCAL = STATE.UCAL;
	H.FLAG.OVERFLOWDETECTION = STATE.OVERFLOWDETECTION;
	H.FLAG.EOG_CORRECTION = STATE.EOG_CORRECTION;
	H.FLAG.OUTPUT = STATE.OUTPUT;

if strncmp(H.TYPE,'IMAGE:',5)
	[H] = iopen(H);
	if H.FILE.OPEN,
		signal = iread(H);
		H.FILE.OPEN = 0;
		fclose(H.FILE.FID);
	end;
	return;
end;

%%%%%%%%%%%%%%% --------- Load single file ------------%%%%%%%%%%%
H = sopen(H,'r',CHAN,MODE);
if 0, ~isnan(H.NS),
%------ ignore 'NaC'-channels
	NS = size(H.Calib,2);
	SelMx = speye(NS);
	ch = 1:NS; %ch(strmatch('NaC',H.Label))=[];
	if length(ch)<NS,
		fprintf(2,'Warning SLOAD: Some NaC channels have been removed %s\n',H.FileName);
	end;
	SelMx = SelMx(:,ch);
	H.Calib = H.Calib*SelMx;

	% generate correction matrix
	if H.FLAG.EOG_CORRECTION,
		H.Calib = H.Calib*h.REGRESS.r0;
		if ~isequal(H.Label(ch),h.Label)
			fprintf(2, 'Warning SLOAD (EOG correction): Channel labels in %s do not fit with arti* recording!\n',H.FileName);
		end;
	elseif STATE.EOG_CORRECTION,
		fprintf(2, 'Warning: EOG correction not supported for this file (%s)!\n',H.FileName);
	end;

	%----- generate HDR.Calib -----
	if all(size(CHAN)>1) || any(floor(CHAN)~=CHAN) || any(CHAN<0) || (any(CHAN==0) && (numel(CHAN)>1));
        	ReRefMx = CHAN;
        	CHAN = find(any(CHAN,2));
	elseif all(CHAN>0) && all(floor(CHAN)==CHAN),
		if any(diff(CHAN)<=0),
		%	fprintf(HDR.FILE.FID,'Warning SOPEN: CHAN-argument not sorted - header information like Labels might not correspond to data.\n');
		end;
        	ReRefMx = sparse(CHAN,1:length(CHAN),1,H.NS,length(CHAN));
	else %if (CHAN==0)
		ReRefMx = [];
	end
	if ~isempty(ReRefMx),
		%[size(H.Calib),size(SelMx),size(h.REGRESS.r0),size(ReRefMx)]
		ReRefMx = [ReRefMx(1:min(size(ReRefMx,1),size(H.Calib,2)),:); zeros(max(0,size(H.Calib,2)-size(ReRefMx,1)),size(ReRefMx,2))];
		H.Calib = H.Calib*ReRefMx;
	end;

	H.InChanSelect = find(any(H.Calib(2:end,:),2));
	H.Calib = H.Calib([1;1+H.InChanSelect(:)],:);
end

if 0,

elseif isfield(H,'data') && all(diag(H.Calib(2:end,1:end))==1) && ~H.FLAG.OVERFLOWDETECTION,
	% only a single copy of the data
	% important for large data sets, close to the available memory
	signal = H.data;
	H = rmfield(H,'data');
	if (CHAN>0)
		signal = signal(:,CHAN);
	end;

elseif any(strmatch(H.TYPE,{'native','TFM_EXCEL_Beat_to_Beat','EEProbe-CNT','EEProbe-AVR'}));
	[signal,H] = sread(H);
        H = sclose(H);

elseif (H.FILE.OPEN > 0)
	[signal,H] = sread(H);
        H = sclose(H);

elseif (H.FILE.OPEN > 0)
        signal = repmat(NaN,H.SPR*H.NRec,size(H.Calib,2));
	k1 = 0;
        while ~seof(H),
	        [s,H] = sread(H,100);
        	if 1,
			k2 = size(s,1);
			signal(k1+1:k1+k2,:)=s;
			k1 = k1+k2;
		else
	        	 %fprintf('%i/%i\n',stell(H),H.SPR*H.NRec);
		        signal=[signal;s];
		end;
	end;
        H = sclose(H);

elseif strncmp(H.TYPE,'EVENT',5)
        signal = H.EVENT;
        return;

elseif strncmp(H.TYPE,'ELPOS',5)
        signal = H.ELEC.XYZ;
	return;

elseif strcmp(H.TYPE,'DAQ')
	fprintf(1,'Loading a matlab DAQ data file - this can take a while.\n');
	tic;
        [signal, tmp, H.DAQ.T0, H.DAQ.events, DAQ.info] = daqread(H.FileName);
        fprintf(1,'Loading DAQ file finished after %.0f s.\n',toc);
        [H.SPR,H.NS] = size(signal);

        H.SampleRate = DAQ.info.ObjInfo.SampleRate;
        sz     = size(signal);
        if length(sz)==2, sz=[1,sz]; end;
        H.NRec = sz(1);
        H.Dur  = sz(2)/H.SampleRate;
        H.NS   = sz(3);
        H.FLAG.TRIGGERED = H.NRec>1;
        H.FLAG.UCAL = 1;

        H.PhysDim = {DAQ.info.ObjInfo.Channel.Units};
        H.DAQ   = DAQ.info.ObjInfo.Channel;

        H.Cal   = diff(cat(1,DAQ.info.ObjInfo.Channel.InputRange),[],2).*(2.^(-DAQ.info.HwInfo.Bits));
        H.Off   = cat(1,DAQ.info.ObjInfo.Channel.NativeOffset);
        H.Calib = sparse([H.Off';eye(H.NS)]*diag(H.Cal));

        if CHAN<1,
                CHAN = 1:H.NS;
        end;
        if ~H.FLAG.UCAL,
                Calib = H.Calib;	% Octave can not index sparse matrices within a struct
                signal = [ones(size(signal,1),1),signal]*Calib(:,CHAN);
        end;

elseif 0, strcmp(H.TYPE,'BIFF'),
	try,
                [H.TFM.S,H.TFM.E] = xlsread(H.FileName,'Beat-To-Beat');
	        if size(H.TFM.S,1)+1==size(H.TFM.E,1),
                        H.TFM.S = [repmat(NaN,1,size(H.TFM.S,2));H.TFM.S];
                end;
                H.TYPE = 'TFM_EXCEL_Beat_to_Beat';
	catch
		try
	                [H.TFM.S,H.TFM.E] = xlsread(H.FileName,'Beat-to-Beat');
        	        H.TYPE = 'TFM_EXCEL_Beat_to_Beat';
        	catch
			0;
        	end;
	end;

	if strcmp(H.TYPE, 'TFM_EXCEL_Beat_to_Beat');
                if ~isempty(strfind(H.TFM.E{3,1},'---'))
                        H.TFM.S(3,:) = [];
                        H.TFM.E(3,:) = [];
                end;

                H.Label   = strvcat(H.TFM.E{4,:});
                H.PhysDim = strvcat(H.TFM.E{5,:});

                H.TFM.S = H.TFM.S(6:end,:);
                H.TFM.E = H.TFM.E(6:end,:);

                ix = find(isnan(H.TFM.S(:,2)) & ~isnan(H.TFM.S(:,1)));
                H.EVENT.Desc = H.TFM.E(ix,2);
                H.EVENT.POS  = ix;

                if any(CHAN),
			H.TFM.S = H.TFM.S(:,CHAN);
			H.TFM.E = H.TFM.E(:,CHAN);
		end;
		[H.SPR,H.NS] = size(H.TFM.S);
		H.NRec = 1;
		H.THRESHOLD  = repmat([0,NaN],H.NS,1);

		signal  = H.TFM.S;
		signal(signal==0) = NaN;
        end;

elseif strcmp(H.TYPE,'MatrixMarket'),
        H.FILE.FID = fopen(H.FileName,'rt','ieee-le');

    	line = fgetl(H.FILE.FID);

	H.FLAG.Coordinate = ~isempty(strfind(line,'coordinate'));
	H.FLAG.Array 	  = ~isempty(strfind(line,'array'));

	H.FLAG.Complex = ~isempty(strfind(line,'complex'));
	H.FLAG.Real = ~isempty(strfind(line,'real'));
	H.FLAG.Integer = ~isempty(strfind(line,'integer'));
	H.FLAG.Pattern = ~isempty(strfind(line,'pattern'));

	H.FLAG.General = ~isempty(strfind(line,'general'));
	H.FLAG.Symmetric = ~isempty(strfind(line,' symmetric'));
	H.FLAG.SkewSymmetric = ~isempty(strfind(line,'skew-symmetric'));
	H.FLAG.Hermitian = ~isempty(strfind(lower(line),'hermitian'));

	while strncmp(line,'%',1)
        	line = fgetl(H.FILE.FID);
	end;

	[tmp,status] = str2double(line);
	if any(status)
		fprintf(H.FILE.stderr,'SLOAD (MM): invalid size %s\n',line);
	else
		H.MATRIX.Size = tmp;
	end;

	if length(H.MATRIX.Size)==3,
		H.Length = tmp(3);
		signal = sparse([],[],[],tmp(1),tmp(2),tmp(3));
		for k = 1:H.Length,
	        	line = fgetl(H.FILE.FID);
			[tmp,status] = str2double(line);
			if any(status)
				fprintf(H.FILE.stderr,'SLOAD (MM): invalid size %s\n',line);
			elseif length(tmp)==4,
		    		val = tmp(3) + i*tmp(4);
			elseif length(tmp)==3,
				val = tmp(3);
			elseif length(tmp)==2,
				val = 1;
			else
				fprintf(H.FILE.stderr,'SLOAD (MM): invalid size %s\n',line);
			end;

			if H.FLAG.General,
				signal(tmp(1),tmp(2)) = val;
			elseif H.FLAG.Symmetric,
				signal(tmp(1),tmp(2)) = val;
				signal(tmp(2),tmp(1)) = val;
			elseif H.FLAG.SkewSymmetric,
				signal(tmp(1),tmp(2)) = val;
				signal(tmp(2),tmp(1)) =-val;
			elseif H.FLAG.Hermitian,
				signal(tmp(1),tmp(2)) = val;
				signal(tmp(2),tmp(1)) = conj(val);
			else
				fprintf(H.FILE.stderr,'SLOAD (MM): invalid size %s\n',line);
			end;
		end;

	elseif length(H.MATRIX.Size)==2
		H.Length = prod(tmp);
		signal = zeros(H.MATRIX.Size);
		if H.FLAG.General==1,
			[IX,IY]=find(ones(H.MATRIX.Size));
		else
			[IX,IY]=find(cumsum(eye(H.MATRIX.Size)));
		end;

		for k = 1:H.Length,
	        	line = fgetl(H.FILE.FID);
			[tmp,status] = str2double(line);
			if any(status)
				error('SLOAD (MM)');
			elseif length(tmp)==2,
				val=tmp(1) + i*tmp(2);
			elseif length(tmp)==1,
				val=tmp(1);
			else
				fprintf(H.FILE.stderr,'SLOAD (MM): invalid size %s\n',line);
			end;

			signal(IX(k),IY(k)) = val;
			if H.FLAG.Symmetric,
				signal(IY(k),IX(k)) = val;
			elseif H.FLAG.SkewSymmetric,
				signal(IY(k),IX(k)) =-val;
			elseif H.FLAG.Hermitian,
				signal(IY(k),IX(k)) = conj(val);
			else
				fprintf(H.FILE.stderr,'SLOAD (MM): invalid size %s\n',line);
			end;
		end;
        end;
        fclose(H.FILE.FID);

elseif strcmp(H.TYPE,'OFF'),
	        H.FILE.FID = fopen(H.FileName,'rt','ieee-le');

                line1 = fgetl(H.FILE.FID);
                line2 = fgetl(H.FILE.FID);
		while ~feof(H.FILE.FID) && (line2(1)=='#'),
	                line2 = fgetl(H.FILE.FID);
		end;
		[tmp,status] = str2double(line2);
		if status || (size(tmp,2)~=3),
			fclose(H.FILE.FID);
			error('SOPEN (OFF)');
		else
			H.VertexCount = tmp(1);
			H.FaceCount = tmp(2);
			H.EdgeCount = tmp(3);
		end

		H.Vertex = repmat(NaN,H.VertexCount,3);
		for k = 1:H.VertexCount,
			line = '';
			while isempty(line) || strncmp(line,'#',1)
				line = fgetl(H.FILE.FID);
			end;
			len = min(length(line),min(find(line=='#')));
			tmp = str2double(line(1:len));
			H.Vertex(k,:) = tmp(1:H.ND);
		end;

%		H.Face = repmat(NaN,H.FaceCount,3);
		for k = 1:H.FaceCount,
			line = '';
			while isempty(line) || strncmp(line,'#',1)
				line = fgetl(H.FILE.FID);
			end;
			len = min(length(line),min(find(line=='#')));
			tmp = str2double(line(1:len));
			H.Ngon(k) = tmp(1);
			H.Face{k} = tmp(2:tmp(1)+1) + 1;
		end;
		if all(H.Ngon(1)==H.Ngon),
			H.Face = cat(1,H.Face{:});
		end;
                fclose(H.FILE.FID);

elseif strcmp(H.TYPE,'POLY'),
        H.FILE.FID = fopen(H.FileName,'rt','ieee-le');

	K = 0;
	while ~feof(H.FILE.FID)
        	line = fgetl(H.FILE.FID);
		if isempty(line),
		elseif line(1)=='#',
		else
			K = K + 1;

		end;
        end;
        fclose(H.FILE.FID);

elseif strcmp(H.TYPE,'SMF'),
	        H.FILE.FID = fopen(H.FileName,'rt','ieee-le');

		VertexCount = 0;
		FaceCount = 0;
		PalLen = 0;
		K = 1;
		while ~feof(H.FILE.FID)
	                line = fgetl(H.FILE.FID);
			if isempty(line)
			elseif line(1)=='#';

			elseif line(1)=='v';
				[tmp,status] = str2double(line(3:end));
				if ~any(status)
					VertexCount = VertexCount + 1 ;
					H.Vertex(VertexCount,:) = tmp;
				else
					fprintf(H.FILE.stderr,'Warning SLOAD: could not read line %i in file %s\n',K,H.FileName);
				end;

			elseif line(1)=='f';
				[tmp,status] = str2double(line(3:end));
				if ~any(status)
					FaceCount  = FaceCount + 1;
					H.Ngon(FaceCount) = length(tmp);
					H.Face{FaceCount} = tmp;
				else
					fprintf(H.FILE.stderr,'Warning SLOAD: could not read line %i in file %s\n',K,H.FileName);
				end;

			elseif line(1)=='n';
				[tmp,status] = str2double(line(3:end));
				if ~any(status)
					H.NormalVector = tmp;
				else
					fprintf(H.FILE.stderr,'Warning SLOAD: could not read line %i in file %s\n',K,H.FileName);
				end;

			elseif line(1)=='c';
				[tmp,status] = str2double(line(3:end));
				if ~any(status)
					PalLen = PalLen +1;
					H.Palette(PalLen,:)= tmp;
				else
					fprintf(H.FILE.stderr,'Warning SLOAD: could not read line %i in file %s\n',K,H.FileName);
				end;
			else

			end;
			K = K+1;
		end;
		fclose(H.FILE.FID);
		if all(H.Ngon(1)==H.Ngon),
			H.Face = cat(1,H.Face{:});
		end;

elseif strcmp(H.TYPE,'TVF 1.1A'),
        H.FILE.FID = fopen(H.FileName,'rt');

	tmp = fgetl(H.FILE.FID);
	H.TVF.Name = fgetl(H.FILE.FID);
	H.TVF.Desc = fgetl(H.FILE.FID);
	[tmp,status] = str2double(fgetl(H.FILE.FID));
	H.TVF.NTR = tmp(1);
	H.TVF.NV  = tmp(2);
	H.FLAG.CULL = ~~tmp(3);
	[tmp,status] = str2double(fgetl(H.FILE.FID));
	H.TVF.NTRC  = tmp(1);
	H.TVF.NTRM  = tmp(2);
	[tmp,status] = str2double(fgetl(H.FILE.FID));
	H.TVF.NVC   = tmp(1);
	H.TVF.NVN   = tmp(2);
	[tmp,status] = str2double(fgetl(H.FILE.FID));
	H.TVF.GlobalColor = tmp;
	[tmp,status] = str2double(fgetl(H.FILE.FID));
	H.TVF.GlobalMtrlProps = tmp;

	H.TVF.Triangles = repmat(NaN,[H.TVF.NTR,3]);
	for k = 1:H.TVF.NTR,
		[tmp, status] = str2double(fgetl(H.FILE.FID));
		H.TVF.Triangles(k,:) = tmp;
	end;
	H.TVF.TrColorSets = reshape(NaN,[H.TVF.NTRC,H.TVF.NTR]);
	for k = 1:H.TVF.NTRC,
		[tmp, status] = str2double(fgetl(H.FILE.FID));
		H.TVF.TrColorSets(k,:) = tmp;
	end;
	H.TVF.TrMtrlSets = repmat(NaN, [H.TVF.NTRM,H.TVF.NTR]);
	for k = 1:H.TVF.NTRM,
		[tmp, status] = str2double(fgetl(H.FILE.FID));
		H.TVF.TrMtrlSets(k,:) = tmp;
	end;

	H.TVF.Vertices   = repmat(NaN, [H.TVF.NV,3]);
	for k = 1:H.TVF.NV,
		[tmp, status] = str2double(fgetl(H.FILE.FID));
		H.TVF.Vertices(k,:) = tmp;
	end;
	H.TVF.VColorSets = repmat(NaN, [H.TVF.NVC,H.TVF.NV]);
	for k = 1:H.TVF.NVC,
		[tmp, status] = str2double(fgetl(H.FILE.FID));
		H.TVF.VColorSets(k,:) = tmp;
	end;
	H.TVF.VNrmlSets  = repmat(NaN, [H.TVF.NVN,3]);
	for k = 1:H.TVF.NVN,
		[tmp, status] = str2double(fgetl(H.FILE.FID));
		H.TVF.VNrmlSets(k,:) = tmp;
	end;

	fclose(H.FILE.FID);

elseif strcmp(H.TYPE,'TVF 1.1B'),
        H.FILE.FID = fopen(H.FileName,'rb',H.Endianity);

	tmp = fread(H.FILE.FID,12,'uchar');
	H.TVF.Name = char(fread(H.FILE.FID,32,'uchar'));
	H.TVF.Desc = char(fread(H.FILE.FID,80,'uchar'));
	tmp = fread(H.FILE.FID,7,'uint32');
	H.TVF.NTR = tmp(1);
	H.TVF.NV  = tmp(2);
	H.FLAG.CULL = ~~tmp(3);
	H.TVF.NTRC  = tmp(4);
	H.TVF.NTRM  = tmp(5);
	H.TVF.NVC   = tmp(6);
	H.TVF.NVN   = tmp(7);
	tmp = fread(H.FILE.FID,[4,2],'float32');
	H.TVF.GlobalColor = tmp(:,1)';
	H.TVF.GlobalMtrlProps = tmp(:,2)';

	H.TVF.Triangles   = fread(H.FILE.FID, [3,H.TVF.NTR],'uint32')';
	H.TVF.TrColorSets = fread(H.FILE.FID, [H.TVF.NTR,H.TVF.NTRC],'uint32')';
	H.TVF.TrMtrlSets  = fread(H.FILE.FID, [H.TVF.NTR,H.TVF.NTRM],'uint32')';

	H.TVF.Vertices   = fread(H.FILE.FID, [3,H.TVF.NV],'float32')';
	H.TVF.VColorSets = fread(H.FILE.FID, [H.TVF.NV,H.TVF.NVC],'uint32')';
	H.TVF.VNrmlSets  = fread(H.FILE.FID, [3,H.TVF.NVN],'float32')';

	fclose(H.FILE.FID);

elseif strcmp(H.TYPE,'VTK'),
                H.FILE.FID = fopen(H.FileName,'rt','ieee-le');

                H.VTK.version = fgetl(H.FILE.FID);
                H.VTK.Title   = fgetl(H.FILE.FID);
                H.VTK.type    = fgetl(H.FILE.FID);

		while ~feof(H.FILE.FID),
			tline = fgetl(H.FILE.FID);

			if 0,
			elseif strncmp(tline,'CELLS',5);
			elseif strncmp(tline,'CELL_DATA',9);
			elseif strncmp(tline,'DATASET',7);
		                H.VTK.DATASET = tline(8:end);
			elseif strncmp(tline,'LINES',6);
			elseif strncmp(tline,'POINTS',6);
			elseif strncmp(tline,'POINT_DATA',10);
			elseif strncmp(tline,'POLYGON',7);
			elseif strncmp(tline,'SCALARS',7);
				[t1,r]=strtok(tline);
				[dataName,r]=strtok(r);
				[dataType,r]=strtok(r);
				[numComp ,r]=strtok(r);
				if isempty(numComp), numComp=1;
				else numComp = str2double(numComp); end;
				tline = fgetl(fid);
				if strcmp(tline,'LOOKUP_TABLE');
	    				[t1,r]=strtok(tline);
	    				[tableName,r]=strtok(tline);
				else
					tline = fgetl(fid);
				end;
			end;

		end;
                fclose(H.FILE.FID);

                fprintf(H.FILE.stderr,'Warning SOPEN: VTK-format not supported, yet.\n');

elseif strcmp(H.TYPE,'unknown')
	[p,f,e]=fileparts(H.FileName);
        TYPE = upper(e);
        if 0,
        elseif strcmp(TYPE,'.RAW')
                loadraw;
        elseif strcmp(TYPE,'.RDT')
                [signal] = loadrdt(FILENAME,CHAN);
                fs = 128;
        elseif strcmp(TYPE,'.XLS')
                loadxls;
        elseif strcmp(TYPE,'.DA_')
                fprintf('Warning SLOAD: Format DA# in testing state and is not supported\n');
                loadda_;
        elseif strcmp(TYPE,'.RG64')
                [signal,H.SampleRate,H.Label,H.PhysDim,H.NS]=loadrg64(FILENAME,CHAN);
                %loadrg64;
        else
                fprintf('Error SLOAD: Unknown Data Format\n');
                signal = [];
                return;
        end;
end;

end;

%%%%%%%%%% Post-Processing %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(H.TYPE,'CNT');
        f = fullfile(H.FILE.Path, [H.FILE.Name,'.txt']);
        if exist(f,'file'),
                fid = fopen(f,'r');
		tmp = fread(fid,inf,'uint8');
		fclose(fid);
		[tmp,v] = str2double(char(tmp'));
		if ~any(v),
            		H.Classlabel=tmp(:);
	        end;
        end
        f = fullfile(H.FILE.Path, [H.FILE.Name,'.par']);
        if exist(f,'file'),
                fid = fopen(f,'r');
		tmp = fread(fid,inf,'uint8');
		fclose(fid);
		[tmp,v] = str2double(char(tmp'));
		if ~any(v),
            		H.Classlabel=tmp(:);
	        end;
        end
        f = fullfile(H.FILE.Path, [H.FILE.Name,'.mat']);
        if exist(f,'file'),
                try,
                        tmp = load(f);
                catch
                        tmp = [];
                end
                if isfield(tmp,'classlabel') && ~isfield(H,'Classlabel')
                        H.Classlabel=tmp.classlabel(:);
                elseif isfield(tmp,'classlabel') && isfield(tmp,'header') && isfield(tmp.header,'iniFile') && strcmp(tmp.header.iniFile,'oom.ini'), %%% for OOM study only.
                        H.Classlabel=tmp.classlabel(:);
                end;
        end;
        f = fullfile(H.FILE.Path, [H.FILE.Name,'c.mat']);
        if exist(f,'file'),
                try,
                        tmp = load(f);
                catch
                        tmp = [];
                end
                if isfield(tmp,'classlabel') && ~isfield(H,'Classlabel')
                        H.Classlabel=tmp.classlabel(:);
                end;
        end;
        f = fullfile(H.FILE.Path, [H.FILE.Name,'_classlabel.mat']);
        if exist(f,'file'),
                try,
                        tmp = load(f);
                catch
                        tmp = [];
                end
                if isfield(tmp,'Classlabel') && (size(tmp.Classlabel,2)==4)
                        [x,H.Classlabel] = max(tmp.Classlabel,[],2);
                end;
                if isfield(tmp,'classlabel') && (size(tmp.classlabel,2)==4)
                        [x,H.Classlabel] = max(tmp.classlabel,[],2);
                end;
        end;

        f=fullfile(H.FILE.Path,[H.FILE.Name,'.sel']);
        if ~exist(f,'file'),
                f=fullfile(H.FILE.Path,[H.FILE.Name,'.SEL']);
        end
        if exist(f,'file'),
                fid = fopen(f,'r');
		tmp = fread(fid,inf,'uint8');
		fclose(fid);
		[tmp,v] = str2double(char(tmp'));
		if any(isnan(tmp)) || any(tmp~=ceil(tmp)) || any(tmp<0) || (any(tmp==0) && any(tmp>1))
                        fprintf(2,'Warning SLOAD(CNT): corrupted SEL-file %s\n',f);
                else
                        if ~isfield(H,'Classlabel'),
                                H.Classlabel = H.EVENT.TYP;
                        end;
                        n = length(H.Classlabel);

                        H.ArtifactSelection = zeros(n,1);
                        if all((tmp==0) | (tmp==1)) && (length(tmp)>1) && (sum(diff(sort(tmp))~=0) ~= length(tmp)-1)
                                H.ArtifactSelection = logical(tmp);
                        else
                                H.ArtifactSelection(tmp) = 1;
                        end;
                end;
        end;
end;

if ~isempty(strfind(upper(MODE),'TSD'));
        f = fullfile(H.FILE.Path, [H.FILE.Name,'.tsd']);
        if ~exist(f,'file'),
                        fprintf(2,'Warning SLOAD-TSD: file %s.tsd found\n',H.FILE(1).Name,H.FILE(1).Name);
        else
                fid = fopen(f,'rb');
                tsd = fread(fid,inf,'float');
                fclose(fid);
                nc = size(signal,1)\size(tsd,1);
                if (nc == round(nc)),
                        signal = [signal, reshape(tsd,nc,size(tsd,1)/nc)'];
                else
                        fprintf(2,'Warning SLOAD: size of %s.tsd does not fit to size of %s.bkr\n',H.FILE(1).Name,H.FILE(1).Name);
                end;
        end;
end;

if strcmp(H.TYPE,'GDF')
        if isfield(H.EVENT,'TYP')
                if isempty(H.EVENT.TYP)
                        %%%%% if possible, load Reinhold's configuration files
                        f = fullfile(H.FILE.Path, [H.FILE.Name,'.mat']);
                        if exist(f,'file'),
                                try
                                        x = load(f,'header');
                                catch;
                                        x=[];
                                end;
                                if isfield(x,'header'),
                                        if isfield(x.header,'Paradigm')
                                                fprintf(2,'Warning SLOAD (GDF): Obsolete feature is used (header information loaded from MAT-file %s).\n',H.FileName);
                                                fprintf(2,'This feature will be removed in future. Contact <a.schloegl@ieee.org> if you need this.\n ');
                                                H.BCI.Paradigm = x.header.Paradigm;
                                                if isfield(H.BCI.Paradigm,'TriggerTiming');
                                                        H.TriggerOffset = H.BCI.Paradigm.TriggerTiming;
                                                elseif isfield(H.BCI.Paradigm,'TriggerOnset');
                                                        H.TriggerOffset = H.BCI.Paradigm.TriggerOnset;
                                                end;

                                                if isfield(H,'Classlabel') && isempty(H.Classlabel),
                                                        H.Classlabel = x.header.Paradigm.Classlabel;
                                                end;
                                        end;
                                end;
                        end;
                end;
        end;

        fid=fopen(fullfile(H.FILE.Path,[H.FILE.Name,'.sel']),'r');
        if fid<0,
	        fid=fopen(fullfile(H.FILE.Path,[H.FILE.Name,'.SEL']),'r');
        end
        if fid>0,
                [tmp,c] = fread(fid,[1,inf],'uint8');
                fclose(fid);
                [tmp,v,sa] = str2double(tmp);
                if isempty(sa) || isempty(sa{1})
                        H.ArtifactSelection = repmat(0,length(H.TRIG),1);
                elseif any(isnan(tmp)) || any(tmp~=ceil(tmp)) || any(tmp<0) || (any(tmp==0) && any(tmp>1))
                        fprintf(2,'Warning SLOAD(GDF): corrupted SEL-file %s\n',fullfile(H.FILE.Path,[H.FILE.Name,'.sel']));
		else
                        H.ArtifactSelection = zeros(length(H.TRIG),1);
                        if all((tmp==0) | (tmp==1)) && (length(tmp)>1) && (sum(diff(sort(tmp))~=0) ~= length(tmp)-1)
                                H.ArtifactSelection = logical(tmp);
                        else
                                H.ArtifactSelection(tmp) = 1;
                        end;
                end;
        end;
end;

% resampling
if ~isnan(Fs) && (H.SampleRate~=Fs);
        tmp = ~mod(H.SampleRate,Fs) || ~mod(Fs,H.SampleRate);
        tmp2= ~mod(H.SampleRate,Fs*2.56);
        if tmp,
                signal = rs(signal,H.SampleRate,Fs);
                H.EVENT.POS = H.EVENT.POS/H.SampleRate*Fs;
                if isfield(H.EVENT,'DUR');
                        H.EVENT.DUR = H.EVENT.DUR/H.SampleRate*Fs;
                end;
                if isfield(H,'TRIG');
                        H.TRIG = H.TRIG/H.SampleRate*Fs;
                end;
                H.SampleRate = Fs;
        elseif tmp2,
                x = load('resample_matrix.mat');
                signal = rs(signal,x.T256100);
                if H.SampleRate*100~=Fs*256,
                        signal = rs(signal,H.SampleRate/(Fs*2.56),1);
                end;
                H.EVENT.POS = H.EVENT.POS/H.SampleRate*Fs;
                if isfield(H.EVENT,'DUR');
                        H.EVENT.DUR = H.EVENT.DUR/H.SampleRate*Fs;
                end;
                if isfield(H,'TRIG');
                        H.TRIG = H.TRIG/H.SampleRate*Fs;
                end;
                H.SampleRate = Fs;
        else
                fprintf(2,'Warning SLOAD: resampling %f Hz to %f Hz not implemented.\n',H.SampleRate,Fs);
        end;
end;

return;
ratiomissing = mean(isnan(signal));
if any(ratiomissing>.1)
	fprintf(2,'Warning SLOAD: ratio of missing samples exceeds 10%% in file %s.\n',H.FileName);
	ix = find(ratiomissing);
	fprintf(1,'#%3i:  %4.1f%%\n',[ix;ratiomissing(ix)*100])
end;


text:
% MP_SOLVE_MAXWELL_SRC_2D: Solve the 2d Maxwell source problem with a B-spline discretization, in a multipatch domain.
%
% Example to solve the problem
%
%    curl ( epsilon(x) curl (u)) + mu u = f      in Omega
%              (epsilon(x) curl(u)) x n = g      on Gamma_N
%                                 u x n = h x n  on Gamma_D
%
% where the domain \Omega is formed by several patches of the form F((0,1)^n).
%
% USAGE:
%
%  [geometry, msh, space, u] = mp_solve_maxwell_src (problem_data, method_data)
%
% INPUT:
%
%  problem_data: a structure with data of the problem. It contains the fields:
%    - geo_name:     name of the file containing the geometry
%    - nmnn_sides:   sides with Neumann boundary condition (may be empty)
%    - drchlt_sides: sides with Dirichlet boundary condition
%    - c_elec_perm:  electric permittivity (epsilon in the equation)
%    - c_magn_perm:  magnetic permeability (mu in the equation)
%
%  method_data : a structure with discretization data. Its fields are:
%    - degree:     degree of the spline functions.
%    - regularity: continuity of the spline functions.
%    - nsub:       number of subelements with respect to the geometry mesh
%                   (nsub=1 leaves the mesh unchanged)
%    - nquad:      number of points for Gaussian quadrature rule
%
% OUTPUT:
%
%  geometry:  array of geometry structures (see mp_geo_load)
%  msh:       multipatch mesh, consisting of several Cartesian meshes (see msh_multipatch)
%  space:     multipatch space, formed by several tensor product spaces plus connectivity and orientation (see sp_multipatch)
%  u:         the computed degrees of freedom
%
% Copyright (C) 2010-2015 Rafael Vazquez
%
%    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 [geometry, msh, space, u] = mp_solve_maxwell_src (problem_data, method_data)

% Extract the fields from the data structures into local variables
data_names = fieldnames (problem_data);
for iopt  = 1:numel (data_names)
  eval ([data_names{iopt} '= problem_data.(data_names{iopt});']);
end
data_names = fieldnames (method_data);
for iopt  = 1:numel (data_names)
  eval ([data_names{iopt} '= method_data.(data_names{iopt});']);
end

% Construct geometry structure
[geometry, boundaries, interfaces, ~, boundary_interfaces] = mp_geo_load (geo_name);
npatch = numel (geometry);

msh = cell (1, npatch);
sp = cell (1, npatch);
for iptc = 1:npatch
  [knots, zeta] = ...
         kntrefine (geometry(iptc).nurbs.knots, nsub-1, degree, regularity);
  [knots_hcurl, degree_hcurl] = knt_derham (knots, degree);

% Construct msh structure
  rule      = msh_gauss_nodes (nquad);
  [qn, qw]  = msh_set_quad_nodes (zeta, rule);
  msh{iptc} = msh_cartesian (zeta, qn, qw, geometry(iptc));

% Construct space structure
  scalar_spaces = cell (msh{iptc}.ndim, 1);
  for idim = 1:msh{iptc}.ndim
    scalar_spaces{idim} = sp_bspline (knots_hcurl{idim}, degree_hcurl{idim}, msh{iptc});
  end
  sp{iptc} = sp_vector (scalar_spaces, msh{iptc}, 'curl-preserving');
  clear scalar_spaces
end

msh = msh_multipatch (msh, boundaries);
space = sp_multipatch (sp, msh, interfaces, boundary_interfaces);
clear sp

stiff_mat = op_curlu_curlv_mp (space, space, msh, c_stiff);
mass_mat = op_u_v_mp (space, space, msh, c_mass);
rhs = op_f_v_mp (space, msh, f);

% Apply Neumann boundary conditions
Nbnd = cumsum ([0, boundaries.nsides]);
for iref = nmnn_sides
  iref_patch_list = Nbnd(iref)+1:Nbnd(iref+1);
  gref = @(varargin) g(varargin{:},iref);
  rhs_nmnn = op_f_v_mp (space.boundary, msh.boundary, gref, iref_patch_list);
  rhs(space.boundary.dofs) = rhs(space.boundary.dofs) + rhs_nmnn .* space.boundary.boundary_orientation.';
end

% Apply Dirichlet boundary conditions
u = zeros (space.ndof, 1);
[u_drchlt, drchlt_dofs] = sp_drchlt_l2_proj (space, msh, h, drchlt_sides);
u(drchlt_dofs) = u_drchlt;
int_dofs = setdiff (1:space.ndof, drchlt_dofs);

% Solve the linear system
rhs(int_dofs) = rhs(int_dofs) - stiff_mat(int_dofs, drchlt_dofs)*u_drchlt ...
                              - mass_mat(int_dofs, drchlt_dofs)*u_drchlt;

u(int_dofs) = (stiff_mat(int_dofs, int_dofs) + ...
                mass_mat(int_dofs, int_dofs)) \ rhs(int_dofs);

end

%!demo
%! ex_maxwell_src_Lshaped_mp

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
clear all, close all, clc
A = imread('../../CH01_SVD/DATA/dog.jpg');
B = rgb2gray(A);

%% Wavelet decomposition (2 level)
n = 2; w = 'db1'; [C,S] = wavedec2(B,n,w);

% LEVEL 1
A1 = appcoef2(C,S,w,1); % Approximation
[H1 V1 D1] = detcoef2('a',C,S,k); % Details
A1 = wcodemat(A1,128);
H1 = wcodemat(H1,128);
V1 = wcodemat(V1,128);
D1 = wcodemat(D1,128);

% LEVEL 2
A2 = appcoef2(C,S,w,1); % Approximation
[H2 V2 D2] = detcoef2('a',C,S,k); % Details
A2 = wcodemat(A2,128);
H2 = wcodemat(H2,128);
V2 = wcodemat(V2,128);
D2 = wcodemat(D2,128);

dec2 = [A2 H2; V2 D2];
dec1 = [imresize(dec2,size(H1)) H1 ; V1 D1];
image(dec1);
colormap gray
axis off, axis tight
set(gcf,'Position',[100 100 600 800])

author: dynamicslab


text:
function value = month_length_julian ( y, m )

%*****************************************************************************80
%
%% MONTH_LENGTH_JULIAN returns the number of days in a Julian month.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year in which the month occurred.
%
%    Input, integer M, the number of the month.
%
%    Output, integer VALUE, the number of days
%    in the month.
%
  mdays = [ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 ];
%
%  Copy the input.
%
  m2 = m;
  y2 = y;
%
%  Check the input.
%
  [ y2, m2, ierror ] = ym_check_julian ( y2, m2 );

  if ( ierror ~= 0 )
    value = 0;
    return
  end
%
%  Get the number of days in the month.
%
  value = mdays(m2);
%
%  If necessary, add 1 day for February 29.
%
  if ( m2 == 2 && year_is_leap_julian ( y2 ) )
    value = value + 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_bug905

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_freqanalysis ft_specest_mtmfft

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug905.mat'));

timeunitinsec=10^(-4); % in seconds
fs=1/(4.8*timeunitinsec);
timestep=1/fs;

windowlength=2;
windowlengthinsamples=ceil(windowlength/timestep);
windowlengthinsec=timestep*windowlengthinsamples;

% determine the fourier frequencies
ralfreq=1/windowlengthinsec;
minfreq=ralfreq;
maxfreq=150;
freqs=minfreq:ralfreq:maxfreq;
nfreq=length(freqs);

smoothfraction=0.50;
smoothfreqs=freqs*smoothfraction;
ntaps=floor(smoothfreqs*windowlengthinsec - 1);
% identify the first frequency bin
selvec=(ntaps<=1);
freqsthisbin=freqs(selvec);
freqbins{1}=[freqsthisbin(1),freqsthisbin(end)];
% identify the remaining frequency bins
% we make frequency bins of increasing width
nfreqsinfreqbin=length(freqsthisbin);
lastfreqindxthisbin=nfreqsinfreqbin;
finished=false;
freqbinindx=1;
while ~finished
    freqbinindx=freqbinindx+1;
    nfreqsinfreqbin=nfreqsinfreqbin+2;
    firstfreqindxthisbin=lastfreqindxthisbin+1;
    lastfreqindxthisbin=lastfreqindxthisbin+nfreqsinfreqbin;
    if lastfreqindxthisbin>=nfreq
        lastfreqindxthisbin=nfreq;
        finished=true;
    end
    freqsthisbin=freqs(firstfreqindxthisbin:lastfreqindxthisbin);
    freqbins{freqbinindx}=[freqsthisbin(1),freqsthisbin(end)];
end
nfreqbins=length(freqbins);
tapsmofreqs=zeros(1,nfreqbins);
freqbinindcs=zeros(1,nfreq);
timeresol=1/fs;
mintapsmofreq=1/(windowlengthinsec - timeresol);
for freqbinindx=1:nfreqbins
    tmpval=mean(freqbins{freqbinindx})*smoothfraction/2;
    tapsmofreqs(freqbinindx)=max([tmpval,mintapsmofreq]);
    freqindx1=find(freqs==freqbins{freqbinindx}(1));
    freqindx2=find(freqs==freqbins{freqbinindx}(2));
    freqbinindcs(freqindx1:freqindx2)=freqbinindx;
end

% perform a frequency analysis
allfreqs=[];
freqout_mtmfft=cell(1,nfreqbins);
for freqbinindx=1:nfreqbins
  freqcfg=[];
  freqcfg.method = 'mtmfft';
  freqcfg.output = 'fourier';
  freqcfg.keeptrials = 'yes';
  freqcfg.foilim = freqbins{freqbinindx};
  freqcfg.tapsmofrq = tapsmofreqs(freqbinindx);
  freqcfg.taper = 'dpss';
  freqcfg.pad='maxperlen';
  freqout_mtmfft{freqbinindx}=ft_freqanalysis(freqcfg,datapart);
  allfreqs=[allfreqs,freqout_mtmfft{freqbinindx}.freq];
end

ntrials=length(datapart.trial);
trialindcs=1:ntrials;

freqdescrout_mtmfft=cell(1,nfreqbins);
for freqbinindx=1:nfreqbins
  freqdescrcfg=[];
  freqdescrcfg.trials=trialindcs;
  freqdescrout_mtmfft{freqbinindx}=ft_freqdescriptives(freqdescrcfg,freqout_mtmfft{freqbinindx});
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [statslice,cesslice] = er_stxgslice(hAvgFile,contrast,varargin);
% [statslice,cesslice] = er_stxgslice(hAvgFile,contrast,[options]);
%
% er_stxgslice: fs-fast stx slice grinder script, ported for use with
% mrLoadRet
%
% Computes one (inplane) slice of a contrast map, across one or more input
% scans. Returns the specified statistic as statslice.
%
% hAvgFile: path to the tSeries for this slice (e.g.
% Inplane/Original/Tseries/Scan1/tSeries1.mat)
%
% contrast: a structure with the same format as that produced by the
% FS-FAST function mkcontrast-sess. It has at least the following fields:
%       ContrastMtx_0: specifies the condition contrast. For block-design
%       (non-deconvolved) experiments, it should have the following format:
%       1 row, nConds columns, where nConds does not include the fixation/
%       baseline condition (usually specified as condition 0). Active
%       conditions should have positive sign, control conditions should
%       have negative sign, conditions which are not to be considered
%       should be 0. The sum of the active conditions should be 1, control
%       conditions should be -1. E.g. for 5 conditions, comparing 1 and 2
%       vs. 4 and 5 and ignoring condition 3:
%
%               contrast.ContrastMtx_0 = [0.5 0.5 0 -0.5 -0.5];
%
% Right now this doesn't save anything, just hands off two matrices,
% statslice and cesslice. statslice contains the results of a statistic, by
% default the -log10 p value from a t test between the active Versus
% control conditions. cesslice is a measure (still not sure the
% mathematical basis) of the size of the contrast effect. In mrLoadRet
% terms, these might have analogous functions to the coherence (statslice)
% and amplitude (cesslice) of the corAnal.
%
% Options: hand in options in pairs, e.g. 'ihDelta',1.25. Can be strings or
% numeric values (if the parameter is numeric). You can set any of the
% parameters set at the top of the code (sorry this is so hands on; in exchange
% I hope to make the code much more legible as I go over several iterations).

% $Id: er_stxgslice.m,v 1.12 2005/08/09 22:58:27 sayres Exp $
% er_stxgslice: ported from fmri_stxgslice by ras, 10/1/03
% 03/05/04 ras: changed from a script to a function. Doesn't save any
% output files, but instead returns a matrix containing the requested stat
% matrix for the given slice.
%
% See also: computeContrastMap,er_stxgrinder,er_mkcontrast.
% ras 04/03/04: fixed the scaling issue (wasn't taking estimated
% variance.^2)
statslice = []; cesslice = [];

% fprintf(1,' ---------- StXGSlice.m : Starting ------------\n');

%%%%% PARAMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
QuitOnError = 1;
OutputFormat = 1; % 0 = log(p) (natural log)
                  %  1 = log10(p)
                  %  2 = p
                  %  3 = test value
ActiveCond = contrast.condNums(find(contrast.WCond > 0));
ControlCond = contrast.condNums(find(contrast.WCond < 0));
CorrIdeal = 0;
CmpIdeal = 0;
TestType = 'tm'; % for list of test types, help fmri_mrestriction
ihDelta = 1.25;
ihTau = 2.5;
dataSize = [];  % [X Y] size of inplanes; enter through varargin if functional is not square (e.g., cropped)
datfile = fullfile(fileparts(hAvgFile),'h.dat'); % this should be produced by selxavg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% parse the option flags
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:length(varargin)
    if (i < length(varargin)) & (ischar(varargin{i}))
		if (isnumeric(varargin{i+1}))
	        cmd = sprintf('%s = [%s];',varargin{i},num2str(varargin{i+1}));
		elseif ischar(varargin{i+1})
	        cmd = sprintf('%s = ''%s'';',varargin{i},varargin{i+1})
		elseif iscell(varargin{i+1});
			tmp = unNestCell(varargin{i+1});
	        cmd = sprintf('%s = {',varargin{i});
			for j = 1:length(tmp)-1
				cmd = [cmd '''' tmp{j} ''','];
			end
			cmd = [cmd '''' tmp{end} '''};']
	    end
        eval(cmd);
	end
end

if (CmpIdeal) | (CorrIdeal) % still trying to figure out what CorrIdeal is -- it's for deconvolved data, at least
    if isempty(ihDelta)
        error('er_stxgslice: need to input ihDelta if correlating w/ Ideal HRF.')
    end

    if isempty(ihTau)
        error('er_stxgslice: need to input ihTau if correlating w/ Ideal HRF.')
    end
end

%%%% ----- Check that variables are the proper values ---- %%%%%
% Check the Test Type %
if ( isempty( strmatch(upper(TestType),{'T ','TM','FM','F0','FD','FC','FCD','FDC'},'exact')))
    fprintf(2,'Error: Unkown TestType %s',TestType);
    if (QuitOnError) quit;
    else            return;
    end
end

if ( OutputFormat < 0 | OutputFormat > 3)
    fprintf(2,'Error: OutputFormat=%d, must be within 0 to 3\n',OutputFormat);
    if (QuitOnError) quit;
    else            return;
    end
end

%% -- check that input files exist -- %%
if ~exist(hAvgFile,'file')
    error(['Can''t find file ' hAvgFile '.']);
end
if ~exist(datfile,'file')
    error(['Can''t find file ' datfile '.']);
end

%% --- Read the h.dat Header File ---- %%%
% fprintf(1,'Reading h.dat Header \n');
hd = fmri_lddat3(datfile);
if isempty(dataSize) % if not entered through varargin
    dataSize = [hd.Nrows hd.Ncols];
end

%% --- Read the hAvg File ---- %%%
% fprintf(1,'Reading hAvg File \n');
load(hAvgFile)
if ~exist('tSeries','var')
    error([hAvgFile ' does not contain a tSeries variable!']);
end

%% --- This should always be 1 until I add deconvolution support ---- %%%
if (hd.GammaFit > 0)
    HDelTest = 1;
else
   newHDelMin = hd.TER*round(HDelMin/hd.TER)+hd.TER*round(hd.TPreStim/hd.TER);
   newHDelMax = hd.TER*round(HDelMax/hd.TER)+hd.TER*round(hd.TPreStim/hd.TER);
   fprintf(1,'Info: newHDelMin = %g, newHDelMax = %g\n',newHDelMin,newHDelMax);
   HDelTest = round([newHDelMin/hd.TER:newHDelMax/hd.TER]') + 1; %'
   fprintf(1,'Info: HDelTest ');
   fprintf(1,' %d',HDelTest);
   fprintf(1,'\n');
end

%% --- From the sxa format of the tSeries, grab hAvg and eVar ---- %%%
% These are the average planes for each condition and the estimated
% residual variance. I may be incorrect--I need to verify this--but I
% believe the format selxavg saves as contains the residual variance as the
% estimated variances for the null, 0, condition. This is the first plane.
% So for block design expts (Nh = 1), the format of the rows of the tSeries
% is: zeros for baseline, eVar, hAvg for Cond 1, variance for Cond 1, hAvg
% for Cond 2, variance for cond 2, etc. er_stxgrinder requires hAvg and
% eVar:
eVar = reshape(tSeries(2,:),dataSize);
eVar = eVar .* eVar;
ind = 3:2:size(tSeries,1);
for i = 1:length(ind)
	hAvg(:,:,i) = reshape(tSeries(ind(i),:),dataSize);
end

if ( exist('HDelMin') & hd.GammaFit > 0)
    msg = 'Cannot specify a delay range with gamma-fit average files';
    qoe(msg); error(msg);
end

if (CorrIdeal & hd.GammaFit > 0)
    msg = 'Cannot correlate ideal HDR with gamma-fit average files';
    qoe(msg); error(msg);
end

nVoxels = hd.Nrows*hd.Ncols;

nPreStim = floor(hd.TPreStim/hd.TER);

% set up params for er_stxgrinder, which does the actual GLM contrast
% computation:
%Ch = fmri_ldbfile(hcovFile);
Ch = hd.hCovMtx;              % covariance matrix
RM = contrast.ContrastMtx_0;  % restriction matrix
nRM = size(RM,2);
nh  = size(hAvg,3);
if (nRM ~= nh)
    msg = sprintf('er_stxgslice: hAvg size (%d) is inconsistent with CMtx (%d)',nh,nRM);
    error(msg);
else
    RM = contrast.ContrastMtx_0;
end

q = zeros(size(hAvg,3),1);

% The meat of the calculation -- create the maps
[vSig pSig ces] = er_stxgrinder(TestType,hAvg,eVar,Ch,hd.DOF,RM,q);

if nargout > 1
    % ---- Compute ces as percent of baseline --- %
    inputdir = fileparts(hAvgFile);
    a = findstr('hAvg',hAvgFile)+4;
    b = findstr('.mat',hAvgFile)-1;
    sliceno = str2num(hAvgFile(a:b));
    hoffsetname = fullfile(inputdir,sprintf('mean_%03d.mat',sliceno));
    hoffset = er_ldtfile(hoffsetname);
    if (isempty(hoffset))
        fprintf('ERROR: could not load %s\n',hoffsetname);
        return;
    elseif size(hoffset,2)==1
        % somehow the mean image (hoffset) gets saved
        % as a column rather than row vector
        hoffset = hoffset';
    end
    indz = find(hoffset==0);
    hoffset(indz) = 10^10;
    nces = size(ces,3);
    cesslice = ces./repmat(reshape(hoffset,size(ces)),[1 1 nces]);
end

if (CmpIdeal)
    pSig = 1 - pSig ;
end

if (strncmp('T',upper(TestType),1))
    SignMask = ((vSig>=0) - (vSig<0));
    pSig = pSig .* SignMask;
else
    SignMask = [];
end

% Dont let pSig = zero (messes up log10)
iz = find(abs(pSig) < 10^-300);
pSig(iz) = sign(pSig(iz)) * 10^-300;
iz = find(pSig == 0);
pSig(iz) = 10^-300; % Have to do this because sign(0) = 0

if (OutputFormat == 0)     pSig = -log(abs(pSig))  .*sign(pSig);
elseif (OutputFormat == 1) pSig = -log10(abs(pSig)).*sign(pSig);
end

statslice = pSig;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function transition_net_multi(varargin)

opts.expDir = fullfile('../data','DCFNet-net-7') ;
opts = vl_argparse(opts, varargin);

modelPath = dir(fullfile(opts.expDir, 'net-epoch-*.mat'));
modelPath = sort({modelPath.name});

[~,out_name,~] = fileparts(opts.expDir);

for i = 1:numel(modelPath)
    load(fullfile(opts.expDir, modelPath{i}));
    net = dagnn.DagNN.loadobj(net) ;
    net = deployDCFNet(net);
    if ~exist('../model_multi', 'dir')
        mkdir('../model_multi');
    end
    save(fullfile('../model_multi', [out_name num2str(i, 'epoch_%d.mat')]), 'net') ;
end
end

author: foolwood

repo: DCFNet

save_path: github-repos/MATLAB/foolwood-DCFNet


text:
function p = prior_corrunif(varargin)
% * PRIOR_R  Correlation prior structure
%
% * Description:
%
%  - REFERENCE: prior for correlation matrix
%    Barnard, J.; McCulloch, R. & li Meng, X.
%    Modelling covariance matrices in terms of standart deviations and correlations
%    with applications to shrinkage. Statistical Sinica, 2000
%
%  - P = PRIOR_CORRUNIF('PARAM1', VALUE1, ...)
%    creates the prior structure in which the
%    named parameters have the specified values. Any unspecified
%    parameters are set to default values.
%
%  - P = PRIOR_CORRUNIF(P, 'PARAM1', VALUE1, ...)
%    modify a prior structure with the named parameters altered
%    with the specified values.
%
%  - Parameters for correlation prior [default]
%     nu       - degree of freedom [15]
%     prior_nu - prior for nu [prior_fixed]
%
%  - some inverse-Wishart properties
%    if W ~ InvWish_d (v, A) then E[W] = A/(v-d-1)
%    nu > = d
%    d = number os species (dimension of the square matrix)
%    the construction for this prior assumes A = I.
%
%  - dimension of the correlation vector: (numberClass^2 - numberClass)/2
%
% * See also
%    PRIOR_*
%
% Copyright (c) 2000-2001,2010 Aki Vehtari
% Copyright (c) 2009,2015 Jarno Vanhatalo
% Copyright (c) 2010 Jaakko Riihimäki
% ------------ 2015 Marcelo Hartmann

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip = inputParser;
  ip.FunctionName = 'PRIOR_CORRUNIF';
  ip.addOptional('p', [], @isstruct);
  ip.addParamValue('nu', 15, @(x) isscalar(x) && x > 0);
  ip.addParamValue('prior_nu', [], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('numberClass', 2, @(x) mod(x, 1) == 0 && x > 1);
  ip.addParamValue('aValue', 1, @(x) isreal(x) &&  x > 0);
  ip.parse(varargin{:});
  p = ip.Results.p;

  if isempty(p)
    init = true;
    p.type = 'CORRUNIF';

  else
    if ~isfield(p, 'type') && ~isequal(p.type, 'CORRUNIF')
      error('First argument does not seem to be a valid prior structure')
    end

    init = false;
  end

  % transformation to the real line (aValue stretches or squeeze the real line).
  p.aValue = ip.Results.aValue;

  % Initialize parameters
  % check the condition nu > numberSpecies
  if init || ~ismember('nu', ip.UsingDefaults)
      p.nu = ip.Results.nu;
  end

  if init || ~ismember('numberClass', ip.UsingDefaults)
      p.numberClass = ip.Results.numberClass;
      p.vectSize = (p.numberClass^2 - p.numberClass)/2;
  end

  if p.nu < (p.numberClass - 1)
      error('Degrees of freedom (nu) should be greater than number of classes')
  end

  % Initialize prior structure
  if init
      p.p = [];
  end
  if init || ~ismember('nu_prior', ip.UsingDefaults)
      p.p.nu = ip.Results.prior_nu;
  end

  if init
    % set functions
    p.fh.pak = @prior_corrunif_pak;
    p.fh.unpak = @prior_corrunif_unpak;
  % p.fh.RealToRho = @prior_corrunif_RealToRho;
    p.fh.lp = @prior_corrunif_lp;
    p.fh.lpg = @prior_corrunif_lpg;
    p.fh.recappend = @prior_corrunif_recappend;
  end

end

function [w, s, h] = prior_corrunif_pak(p)
% This is a mandatory subfunction used for example
% in energy and gradient computations.

  w = [];  s = {};  h = [];

  if ~isempty(p.p.nu)
      w = [w p.nu];
      s = [s; 'R.nu'];
      h = [h 1];
  end

end

function [p, w] = prior_corrunif_unpak(p, w)
% This is a mandatory subfunction used for example
% in energy and gradient computations.

  if ~isempty(p.p.nu)
      i1 = 1;
      p.nu = w(i1);
      w = w(i1 + 1:end);
  end

end

function lp = prior_corrunif_lp(x, p)
% This is a mandatory subfunction used for example
% in energy computations.

  % Evaluating log-prior(R)
  % correlation vector
  rho = x';

  % create entries
  seq = 1:p.vectSize;
  i = ceil(0.5 + 0.5 * sqrt(1 + 8*(seq)));
  j = seq - (i - 2).*(i - 1)/2;
  ind1 = (j - 1) * p.numberClass + i;
  ind2 = (i - 1) * p.numberClass + j;

  % build correlation matrix
  R = eye(p.numberClass);

  % filling elements in lower and upper part
  R([ind1, ind2]) = [rho; rho];
  detR = det(R);

  if eps < detR && ~any(abs(rho) > 1)
      % cholesk decompostion
      L = chol(R, 'lower');

      % parameters of the distribution
      a = 0.5*(p.nu - 1)*(p.numberClass - 1) - 1;
      b = - p.nu/2;

      % building principal submatrices and evaluating log determinant
      sDetlogSub = 0;
      for k = 1:p.numberClass
          A = R;
          A(:, k) = []; A(k, :) = [];
          Lsub = chol(A, 'lower');
          sDetlogSub = sDetlogSub + sum(log(diag(Lsub)));
      end

      % evaluating unnormalized log-prior
      lp = 2 * (a*sum(log(diag(L))) + b*sDetlogSub);

  else
      lp = -Inf;

  end

% adding log-hyperprior(nu)
if ~isempty(p.p.nu)
    lp = lp + p.p.nu.fh.lp(p.nu, p.p.nu) + log(p.nu);
end

end

function lpg = prior_corrunif_lpg(x, p)
% This is a mandatory subfunction used for example
% in gradient computations.

 % taking the correlation vector
 rho = x';

 % creating entries
 seq = 1:p.vectSize;
 i = ceil(0.5 + 0.5 * sqrt(1 + 8*(seq)));
 j = seq - (i - 2).*(i - 1)/2;
 ind1 = (j - 1) * p.numberClass + i;
 ind2 = (i - 1) * p.numberClass + j;

 % building corr matrix
 R = eye(p.numberClass);

 % filling elements in lower and upper part
 R([ind1, ind2]) = [rho; rho];

 % all(eig(R) > 0)
  if ~any(abs(rho) > 1)

     % grad vector
     lpg = ones(1, p.vectSize);

     % cholesk decompostion
     % L = chol(R, 'lower');

     % parameters of the distribution
     a = 0.5*(p.nu - 1)*(p.numberClass - 1) - 1;
     b = - p.nu/2;

     % building principal submatrices and evaluating log determinant
     invR = inv(R);

     % COULD WE USE COFACTOR MATRIX ?
     % Id1 = eye(p.numberClass-1);
     aux = 0;

     for j = 2:p.numberClass
         for i = 1:(j-1)
             sumtrDer = 0;
             sumtrDer = sumtrDer + 2*a*invR(j, i);
             for k = 1:p.numberClass
                 if k == i || k == j
                     % sumtrDer = sumtrDer + 0;
                     continue
                 else
                     A = R;
                     A(:, k) = []; A(k, :) = [];
                     % L = chol(A, 'lower');
                     invAk = inv(A);
                     if k < j && (i == 1 || k > i)
                         m = j-1;
                         sumtrDer = sumtrDer + 2*b*invAk(m, i);
                     elseif k < j && k < i
                         l = i-1;
                         m = j-1;
                         sumtrDer = sumtrDer + 2*b*invAk(m, l);
                     else
                         sumtrDer = sumtrDer + 2*b*invAk(j, i);
                     end
                 end
             end
             aux = aux + 1;
             lpg(aux) = sumtrDer;
         end
     end

 else
     lpg = repmat(-Inf, 1, p.vectSize);
 end

 if ~isempty(p.p.nu)
     lpgnu = p.p.nu.fh.lpg(p.nu, p.p.nu).*p.nu + 1;
     lpg = [lpg lpgnu];
  end
end

function rec = prior_corrunif_recappend(rec, ri, p)
% This subfunction is needed when using MCMC sampling (gp_mc).
% The parameters are not sampled in any case.

rec = rec;
if ~isempty(p.p.nu)
    rec.nu(ri,:) = p.nu;
end

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
% Batch script for sparse prediction
clear all

% --- training method for variable sample size data
method_id = 1;
%case	1 linear_sparse_space_sw
%case	2 linear_sparse_space_vec_sw
%case	3 linear_sparse_stepwise_sw
%case	4 linear_sparse_stepwise_em_sw

% --- Basic Learning Parameters
parm.Ntrain = 2000; % # of total training iteration
parm.Nskip  = 100;	% skip steps for display info

% --- Time delay embedding parameter
parm.Tau   = 1 ;   % Lag time steps
parm.Dtau  = 3 ;   % Number of embedding dimension
parm.Tpred = 0 ;   % Prediction time step : y(t+Tpred) = W * x(t)

% --- Normalization parameter
parm.data_norm  = 1;	% Normalize input and output

% File name of old training result for retraining
old_file = [];

% File name for test data
datafile  = ['./test/test.mat'];
modelfile = ['./test/model'];
newdata   = 1; % = 1: make new data

%parm.Trial  = [1000 200 100];% # of samples in each trial
parm.Trial  = [1000 600];% # of samples in each trial

% --- Make or Load training & test data
if newdata == 0 && exist(datafile,'file')
	load(datafile, ...
		'xdata','ydata','xtest','ytest')
else
	% Training & test data setting
	Tdata  = 1000 ; % # of training data
	Ttest  = 1000 ; % # of test data

	parm.Ydim = 3;
	parm.Xdim = 200;   % Input dim
	parm.Meff = 30;    % Effective Input
	parm.Sdim = 10; % number of sinusoidal input
	parm.Tmin = 50; % minimum period
	parm.Tmax = 200;% maximum period
	parm.SY   = [0.3; 0.1; 0.5];% output moise variance
	parm.Ntrial = length(parm.Trial);

	parm.sy = parm.SY;
	[xdata,ydata,Wout,Weff,Wrdn,Win,Wirr] = ...
		make_train_data(parm, Tdata);

	parm.sy = 0.0;
	[xtest,ytest] = ...
		make_train_data(parm, Ttest, Wout,Wrdn,Win,Wirr);

	if ~isempty(datafile)
		save(datafile, ...
			'xdata','ydata','xtest','ytest','Wout','Weff','parm')
	end
end

T = size(ydata,2);
if T~=size(xdata,2), error('Time of X and Y is not match'); end

if isfield(parm,'data_norm') && parm.data_norm==2
	% Add Bias term
	[M,T,K] = size(xdata);
	xdata = [xdata; ones(1,T,K)];
end;

% Time alignment for prediction using embedding input
[tx,ty] = pred_time_index(xdata,parm);
xdata = xdata(:,tx,:);
ydata = ydata(:,ty,:);
% Adjust sample length for each trial by time embedding
if isfield(parm,'Trial'), parm.Trial = parm.Trial - (T - length(ty)); end;

% Normalize input data
[X,nparm] = normalize_data(xdata, parm.data_norm);
parm.xmean = nparm.xmean;
parm.xnorm = nparm.xnorm;

% Normalize output data
[Y,nparm] = normalize_data(ydata, parm.data_norm);
parm.ymean = nparm.xmean;
parm.ynorm = nparm.xnorm;

% --- Initialization of Model Parameters
if ~isempty(old_file)
	% Start from old result
	load([old_file], 'Model')
else
	Model = [];
end

%profile_on = 1;
%profile_start(profile_on);

%
% --- Sparse estimation
%
switch	method_id
case	1
	[Model, Info] = linear_sparse_space_sw(X, Y, Model, parm);
case	2
	[Model, Info] = linear_sparse_space_vec_sw(X, Y, Model, parm);
case	3
	[Model, Info] = linear_sparse_stepwise_sw(X, Y, Model, parm);
case	4
	[Model, Info] = linear_sparse_stepwise_em_sw(X, Y, Model, parm);
end

%
% --- Estimate prediction error for test data
%

% Time alignment for prediction using embedding input
[tx,ty] = pred_time_index(xtest,parm);
xtest = xtest(:,tx,:);
ytest = ytest(:,ty,:);

if isfield(parm,'data_norm') && parm.data_norm==2
	% Add Bias term
	[M,T,K] = size(xtest);
	xtest = [xtest; ones(1,T,K)];
end;

% Use normalization constant calculated by training data
xtest = normalize_data(xtest, parm.data_norm, parm);

% --- Prediction for test data
ypred = predict_output(xtest, Model, parm);
err   = sum((ytest(:)-ypred(:)).^2)/sum(ytest(:).^2)

if ~isempty(modelfile)
	fsave = [modelfile sprintf('_id%d.mat',method_id)];
	save(fsave, 'Model', 'Info', 'parm');
end

plot_predict

%profile_end(profile_on)
return

author: KamitaniLab

repo: GenericObjectDecoding


text:
function svd_snowfall_test04 ( x )

%*****************************************************************************80
%
%% SVD_SNOWFALL_TEST04 looks at the first 6 modes in the U matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(8,123), the snowfall data.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SVD_SNOWFALL_TEST04\n' );
  fprintf ( 1, '  Look at the first 6 modes in the U matrix.\n' );
  fprintf ( 1, '  Each of these represents a pattern for snowfall over a year.\n' );
  fprintf ( 1, '  The first mode is the pattern that is strongest in the data.\n' );
%
%  Compute the SVD.
%
  [ u, s, v ] = svd ( x );
%
%  Normalize the patterns so that each column has maximum entry 1.
%
  us = r8col_normalize_li ( 8, 8, u );
%
%  Plot the modes.
%
  t = 1 : 8;

  figure ( 4 )

  subplot ( 2, 3, 1 )
  hold on
  plot ( [1,8], [0,0], 'r-' )
  plot ( t, us(:,1), 'b-', 'Linewidth', 3 )
  axis ( [1,8,-1.0,+1.0] )
  grid on
  title ( 'U Mode #1' )
  hold off

  subplot ( 2, 3, 2 )
  hold on
  plot ( [1,8], [0,0], 'r-' )
  plot ( t, us(:,2), 'b-', 'Linewidth', 3 )
  axis ( [1,8,-1.0,+1.0] )
  grid on
  title ( 'U Mode #2' )
  hold off

  subplot ( 2, 3, 3 )
  hold on
  plot ( [1,8], [0,0], 'r-' )
  plot ( t, us(:,3), 'b-', 'Linewidth', 3 )
  axis ( [1,8,-1.0,+1.0] )
  grid on
  title ( 'U Mode #3' )
  hold off

  subplot ( 2, 3, 4 )
  hold on
  plot ( [1,8], [0,0], 'r-' )
  plot ( t, us(:,4), 'b-', 'Linewidth', 3 )
  axis ( [1,8,-1.0,+1.0] )
  grid on
  title ( 'U Mode #4' )
  hold off

  subplot ( 2, 3, 5 )
  hold on
  plot ( [1,8], [0,0], 'r-' )
  plot ( t, us(:,5), 'b-', 'Linewidth', 3 )
  axis ( [1,8,-1.0,+1.0] )
  grid on
  title ( 'U Mode #5' )
  hold off

  subplot ( 2, 3, 6 )
  hold on
  plot ( [1,8], [0,0], 'r-' )
  plot ( t, us(:,6), 'b-', 'Linewidth', 3 )
  axis ( [1,8,-1.0,+1.0] )
  grid on
  title ( 'U Mode #6' )
  hold off

  print ( '-dpng', 'u_modes.png' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plotted U modes 1 to 6 in "u_modes.png".\n' );

  return
end


text:
## Copyright (C) 2013 CarnĂŤ Draug <carandraug@octave.org>
##
## 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/>.

## -*- texinfo -*-
## @deftypefn  {Function File} {@var{cmap} =} ycbcr2rgb (@var{YCbCrmap})
## @deftypefnx {Function File} {@var{RGB} =} ycbcr2rgb (@var{YCbCr})
## @deftypefnx {Function File} {@dots{} =} ycbcr2rgb (@dots{}, [@var{Kb} @var{Kr}])
## @deftypefnx {Function File} {@dots{} =} ycbcr2rgb (@dots{}, @var{standard})
## Convert YCbCr color space to RGB.
##
## The convertion changes the image @var{YCbCr} or colormap @var{YCbCrmap},
## from the YCbCr (luminance, chrominance blue, and chrominance red)
## color space to RGB values.  @var{YCbCr} must be of class double, single,
## uint8, or uint16.
##
## The formula used for the conversion is dependent on two constants, @var{Kb}
## and @var{Kr} which can be specified individually, or according to existing
## standards:
##
## @table @asis
## @item "601" (default)
## According to the ITU-R BT.601 (formerly CCIR 601) standard.  Its values
## of @var{Kb} and @var{Kr} are 0.114 and 0.299 respectively.
## @item "709" (default)
## According to the ITU-R BT.709 standard.  Its values of @var{Kb} and
## @var{Kr} are 0.0722 and 0.2116 respectively.
## @end table
##
## @seealso{hsv2rgb, ntsc2rgb, rgb2hsv, rgb2ntsc, rgb2ycbcr}
## @end deftypefn

function rgb = ycbcr2rgb (ycbcr, standard = "601")
  if (nargin < 1 || nargin > 2)
    print_usage ();
  endif
  rgb = ycbcrfunc ("ycbcr2rgb", ycbcr, standard);
endfunction

%!assert (ycbcr2rgb (rgb2ycbcr (jet (10))), jet (10), 0.00001);

## Copyright (C) 2013 CarnĂŤ Draug <carandraug+dev@gmail.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 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/>.

## Private function for ycbcr2rgb and rgb2ycbcr functions which are
## very similar

function out = ycbcrfunc (func, in, standard)

  img = false; # was input an image?
  if (iscolormap (in))
    ## do nothing, it's a colormap
  elseif (isrgb (in))
    img = true;
    ## we shape it as a colormap (2D matrix) so we can use matrix multiplcation
    nRows = rows (in);
    nCols = columns (in);
    in    = reshape (in, [nRows*nCols 3]);
  else
    error ("%s: input must be a colormap (Nx3) or RGB image (NxMx3)", func);
  endif

  if (ischar (standard))
    if (strcmpi (standard, "601")) # for ITU-R BT.601
      Kb = 0.114;
      Kr = 0.299;
    elseif (strcmpi (standard, "709")) # for ITU-R BT.709
      Kb = 0.0722;
      Kr = 0.2126;
    else
      error ("%s: unknown standard `%s'", func, standard);
    endif
  elseif (isnumeric (standard) && numel (standard) == 2)
    Kb = standard(1);
    Kr = standard(2);
  else
    error ("%s: must specify a standard (string), or Kb and Kr values", func);
  endif

  ## the color matrix for the conversion. Derived from:
  ##    Y  = Kr*R + (1-Kr-Kb)*G + kb*B
  ##    Cb = (1/2) * ((B-Y)/(1-Kb))
  ##    Cr = (1/2) * ((R-Y)/(1-Kr))
  ## It expects RGB values in the range [0 1], and returns Y in the
  ## range [0 1], and Cb and Cr in the range [-0.5 0.5]
  cmat = [  Kr            (1-Kr-Kb)            Kb
          -(Kr/(2-2*Kb)) -(1-Kr-Kb)/(2-2*Kb)   0.5
            0.5          -(1-Kr-Kb)/(2-2*Kr) -(Kb/(2-2*Kr)) ];

  cls = class (in);
  in  = im2double (in);

  ## note that these blocks are the inverse of one another. Changes
  ## in one will most likely require a change on the other
  if (strcmp (func, "rgb2ycbcr"))
    ## convert to YCbCr colorspace
    out = (cmat * in')'; # transpose at the end to get back colormap shape
    ## rescale Cb and Cr to range [0 1]
    out(:, [2 3]) += 0.5;
    ## footroom and headroom will take from the range 16/255 each for Cb and Cr,
    ## and 16/255 and 20/255 for Y. So we have to compress the values of the
    ## space, and then shift forward
    out(:,1)     = (out(:,1) * 219/255) + 16/255;
    out(:,[2 3]) = (out(:,[2 3]) * 223/255) + 16/255;

  elseif (strcmp (func, "ycbcr2rgb"))
    ## just the inverse of the rgb2ycbcr conversion
    in(:,[2 3])  = (in(:,[2 3]) - 16/255) / (223/255);
    in(:,1)      = (in(:,1) - 16/255) / (219/255);
    in(:,[2 3]) -= 0.5;
    out          = (inv (cmat) * in')';
  else
    error ("internal error for YCbCr conversion. Unknown function %s", func);
  endif

  switch (cls)
    case {"single", "double"}
      ## do nothing. All is good
    case "uint8"
      out = im2uint8 (out);
    case "uint16"
      out = im2uint16 (out);
    otherwise
      error ("%s: unsupported image class %s", func, cls);
  endswitch

  if (img)
    ## put the image back together
    out = reshape (out, [nRows nCols 3]);
  endif

endfunction

author: fumin

repo: pencil

save_path: github-repos/MATLAB/fumin-pencil