classdef TRData
% TRData class
% Container for time-resolved spectroscopy data
% Petar Lambrev, 2011
%
% A TRData dataset contains a 2D array of data (Z) dependent on two variables
% (X and Y) plus additional identifying information, such as ID, units, etc.
% The class is intended for use with time-resolved spectroscopy data.
%
% X is a 1*m array containing the spectral axis (wavelength) values,
% Y is an n*1 array containing the time axis values,
% Z is an n*m array storing the measured quantity (abs, fluorescence, etc.).
%
% Objects of class TRData contain the data themselves as well as
% a number of methods (functions) for data manipulation and presentation.
properties
XType = 'Wavelength'; % type of X-axis data, default - 'Wavelength'
XUnit = 'cm-1'; % unit for X-axis data, default - 'cm-1'
YType = 'Time'; % type of Y-axis data, default - 'Time'
YUnit = 'ps'; % unit for Y-axis data, default - 'ps'
ZType = '\DeltaA'; % type of Z data (measured quantity - Abs, fluorescence, etc.)
ZUnit = 'mOD'; % unit for Z data, default - 'mOD'
X = []; % X-axis data (usually wavelength), 1 x m vector
Y = []; % Y-axis data (usually time), n x 1 vector
Z = []; % Z data, m x n matrix
ID = ''; % dataset ID
ExpID = ''; % experiment ID
Comment = ''; % comments
DateTime = now; % date and time
end
methods
function SP = TRData(X,Y,Z,ID)
% Create a dataset (TRData object)
%
% Usage:
% dat = TRData
% Creates a blank dataset
%
% dat = TRData(X, Y, Z)
% Creates a dataset dat with supplied X, Y and Z values
%
% dat = TRData(X, Y, Z, ID)
% Creates a dataset dat with X, Y, Z values and ID
%
% see also: load
if (nargin == 2),
if (length(X)~=size(Z,2)),
error('Number of elements in X does not match the data');
end
if (length(Y)~=size(Z,1)),
error('Number of elements in Y does not match the data');
end
SP.X = X;
SP.Y = Y;
SP.Z = Z;
end %if
if (nargin == 4),
if (length(X)~=size(Z,2)),
error('Number of elements in X does not match the data');
end
if (length(Y)~=size(Z,1)),
error('Number of elements in Y does not match the data');
end
SP.X = X;
SP.Y = Y;
SP.Z = Z;
SP.ID = ID;
end %if
end %Constructor
%% Queries %%
% --------- %
function n = Xn(SP)
% number of X channels (X data points) in the dataset(s)
n = zeros(1,length(SP));
for i = 1:length(SP),
n(i) = length(SP(i).X);
end %for
end %size
function n = Yn(SP)
% number of Y channels (Y data points) in the dataset(s)
n = zeros(1,length(SP));
for i = 1:length(SP),
n(i) = length(SP(i).Y);
end %for
end %size
function res = Xind(SP,Xval)
% index of the X value(s) specified by Xval.
% Xind = dat.Xind(Xval)
% If dat is an array, Xind will also be an array containing indices
% for each dataset in dat.
res = zeros(1,length(SP));
for i = 1:length(SP),
res(i) = find(SP.X == Xval, 1);
end %for
end %Xind
function res = Yind(SP,Yval)
% index of the Y value(s) specified by Yval.
% Yind = dat.Yind(Yval)
% If dat is an array, Yind will also be an array containing indices
% for each dataset in dat.
res = zeros(1,length(SP));
for i = 1:length(SP),
res(i) = find(SP.Y == Yval, 1);
end %for
end %Xind
function res = Transients(SP,Xval)
% Transients (kinetic traces)
% Transients = Dat.Transients(Xval)
% Returns the kinetic traces Z = f(X) corresponding to wavelengths
% specified by Xval.
n = length(SP);
if (n > 1)
res = cell(1,length(SP));
else
res = [];
end
for i = 1:n,
Xind = zeros(1,length(Xval));
for j = 1:length(Xval)
Xind(j) = find(SP(i).X >= Xval(j), 1);
if (isempty(Xind)), error('X value not found in data'); end
end
if n>1
tres = SP(i).Z(:,Xind);
res{i} = [tres];
else
res = SP(i).Z(:,Xind);
end
end %for
end %Transient
function res = Spectra(SP,Yval)
% Transient spectra
% Spectra = Dat.Spectra(Yval)
% Returns the transient spectra Z = f(X) corresponding to time points
% specified by Yval.
n = length(SP);
if (n > 1)
res = cell(1,length(SP));
else
res = [];
end
for i = 1:n,
Yind = zeros(1,length(Yval));
for j = 1:length(Yval)
Yind(j) = find(SP(i).Y >= Yval(j), 1);
if (isempty(Yind)), error('Y value not found in data'); end
end
if n>1
tres = SP(i).Z(Yind,:)';
res{i} = [tres];
else
res = SP(i).Z(Yind,:)';
end
end %for
end %Spectra
%% Arithmetic operators
% ---------------------
function res = plus(SP1,SP2)
% addition of datasets
%
% Usage:
% new = dat1 + dat2
%
% new contains all values of parameters of dat1 except that
% the Z values are added
%
% see also: minus, times, rdivide
if (~isa(SP1,'Spectrum')),
error('First operand must be a Spectrum object');
end %if
res = SP1;
singlemode = 0;
if (length(SP1)~=length(SP2)),
if(length(SP2)==1),
singlemode = 1;
else
error ('Objects must have the same number of elements');
end %if
end %if
if(isa(SP2,'TRData')),
for i = 1:length(SP1)
if (singlemode), j = 1; else j = i; end
if (SP1(i).dim~=SP2(j).dim), error ('Objects must have the same size'); end
if (SP1(i).X~=SP2(j).X), error ('The X-axes do not match'); end
if (SP1(i).Y~=SP2(j).Y), error ('The Y-axes do not match'); end
res(i).Z = SP1(i).Z + SP2(j).Z;
res(i).Comment = sprintf('%s + %s',SP1(i).Comment,SP2(j).ID);
end %for
elseif(isreal(SP2)),
for i = 1:length(SP1)
if (singlemode), j = 1; else j = i; end
res(i).Z = SP1(i).Z + SP2(j);
res(i).Comment = sprintf('%s + %d',SP1(i).Comment,SP2(j));
end %for
else
error('Unrecognized operand data type');
end %if
end %plus
function res = minus(SP1,SP2)
% subtraction of datasets
%
% Usage:
% new = dat1 - dat2
%
% new contains all values of parameters of dat1 except that
% the Z values are subtracted
%
% see also: plus, times, rdivide
if (~isa(SP1,'TRData')),
error('First operand must be a TRData object');
end %if
res = SP1;
singlemode = 0;
if (length(SP1)~=length(SP2)),
if(length(SP2)==1),
singlemode = 1;
else
error ('Objects must have the same number of elements');
end %if
end %if
if(isa(SP2,'TRData')),
for i = 1:length(SP1)
if (singlemode), j = 1; else j = i; end
if (SP1(i).dim~=SP2(j).dim), error ('Objects must have the same size'); end
if (SP1(i).X~=SP2(j).X), error ('The X-axes do not match'); end
if (SP1(i).Y~=SP2(j).Y), error ('The Y-axes do not match'); end
res(i).Z = SP1(i).Z - SP2(j).Z;
res(i).Comment = sprintf('%s - %s',SP1(i).Comment,SP2(j).ID);
end %for
elseif(isreal(SP2)),
for i = 1:length(SP1)
if (singlemode), j = 1; else j = i; end
res(i).Z = SP1(i).Z - SP2(j);
res(i).Comment = sprintf('%s - %d',SP1(i).Comment,SP2(j));
end %for
else
error('Unrecognized operand data type');
end %if
end %minus
function res = uminus(SP)
% unary minus
% Usage:
% new = -dat1
res = SP;
for i =1:length(SP)
res(i).Z = -SP(i).Z;
end
end %uminus
function res = times(SP1,SP2)
% multiplication of datasets
%
% Usage:
% new = dat1 * dat2
%
% new contains all values of parameters of dat1 except that
% the Z values are multiplied
%
% see also: plus, minus, uminus, rdivide
if (~isa(SP1,'TRData')),
error('First operand must be a TRData object');
end %if
res = SP1;
singlemode = 0;
if (length(SP1)~=length(SP2)),
if(length(SP2)==1),
singlemode = 1;
else
error ('Objects must have the same number of elements');
end %if
end %if
if(isa(SP2,'TRData')),
for i = 1:length(SP1)
if (singlemode), j = 1; else j = i; end
if (SP1(i).dim~=SP2(j).dim), error ('Objects must have the same size'); end
if (SP1(i).X~=SP2(j).X), error ('The X-axes do not match'); end
if (SP1(i).Y~=SP2(j).Y), error ('The Y-axes do not match'); end
res(i).Z = SP1(i).Z .* SP2(j).Z;
res(i).Comment = sprintf('%s * %s',SP1(i).Comment,SP2(j).ID);
end %for
elseif(isreal(SP2)),
for i = 1:length(SP1)
if (singlemode), j = 1; else j = i; end
res(i).Z = SP1(i).Z .* SP2(j);
res(i).Comment = sprintf('%s * %d',SP1(i).Comment,SP2(j));
end %for
else
error('Unrecognized operand data type');
end %if
end %times
function res = mtimes(SP1,SP2)
% see also: times
res = times(SP1,SP2);
end %mtimes
function res = rdivide(SP1,SP2)
% division of datasets
%
% Usage:
% new = dat1 / dat2
%
% new contains all values of parameters of dat1 except that
% the Z values are divided
%
% see also: plus, minus, uminus, times
if (~isa(SP1,'TRData')),
error('First operand must be a TRData object');
end %if
res = SP1;
singlemode = 0;
if (length(SP1)~=length(SP2)),
if(length(SP2)==1),
singlemode = 1;
else
error ('Objects must have the same number of elements');
end %if
end %if
if(isa(SP2,'TRData')),
for i = 1:length(SP1)
if (singlemode), j = 1; else j = i; end
if (SP1(i).dim~=SP2(j).dim), error ('Objects must have the same size'); end
if (SP1(i).X~=SP2(j).X), error ('The X-axes do not match'); end
if (SP1(i).Y~=SP2(j).Y), error ('The Y-axes do not match'); end
res(i).Z = SP1(i).Z ./ SP2(j).Z;
res(i).Comment = sprintf('%s / %s',SP1(i).Comment,SP2(j).ID);
end %for
elseif(isreal(SP2)),
for i = 1:length(SP1)
if (singlemode), j = 1; else j = i; end
res(i).Z = SP1(i).Z ./ SP2(j);
res(i).Comment = sprintf('%s / %d',SP1(i).Comment,SP2(j));
end %for
else
error('Unrecognized operand data type');
end %if
end %rdivide
function res = mrdivide(SP1,SP2)
% see also: rdivide
res = rdivide(SP1,SP2);
end %mrdivide
function res = sum(SP)
% sum of datasets
res = SP(1);
if (length(SP)>1),
for i = 2:length(SP)
res.Z = res.Z + SP(i).Z;
end %for
end %if
end %sum
function res = mean(SP)
% average datasets
%
% see also: average
res = average(SP);
end %mean
function res = max(SP)
% maximal Z value
res = max(SP.Z);
end %max
%% Other calculations
% -------------------
function res = smooth(SP,varargin)
% smooth spectra in the datasets
res = SP;
for i = 1:length(SP)
for x = 1:SP.Xn
res(i).Z(:,x) = smooth(SP(i).Z(:,x),varargin{:});
end
end
end %smooth
function res = brushX(SP,Xi)
% discard data in deffective channels
% and replace with interpolated value
%
% Usage:
% NewData = Dat.brushX(Xi)
%
% Erases data in the selected X channel(s) but does not remove
% the channels. Instead, the data are interpolated (linearly)
% from the neighbouring X channels.
% In contrast to dumpX, brushX does not change the X axis.
%
% see also: dumpX, dumpY
res = SP;
for i = 1:length(SP)
% for each channel
for j = 1:length(Xi)
% find 2 nearest neighbours
if (Xi(j) == 1),
j1 = 2; j2 = 3;
elseif (Xi(j) == SP(i).Xn),
j1 = Xi(j)-1; j2 = Xi(j)-2;
else
j1 = Xi(j)-1; j2 = Xi(j)+1;
end %if
% draw line between 2 points
xi = SP(i).X(Xi(j));
x1 = SP(i).X(j1);
x2 = SP(i).X(j2);
y1 = SP(i).Z(:,j1);
y2 = SP(i).Z(:,j2);
lina = (y2 - y1) ./ (x2 - x1);
linb = y1 - lina .* x1;
% find interpolated values
res(i).Z(:,Xi(j)) = lina .* xi + linb;
end %for
end
end %brushX
function res = dumpX(SP,dumpchn)
% delete X (wavelength) channels
%
% Usage:
% NewData = Dat.dumpX(Xi)
%
% Deletes data in the selected X channel(s).
% Useful when specific detector channels are unreliable or unused.
%
% see also: dumpY, brushX
res = SP;
for i = 1:length(SP)
dumpchn = sort(dumpchn,'descend');
res(i).Z(:,dumpchn) = [];
res(i).X(dumpchn) = [];
end
end %dumpX
function res = dumpY(SP,dumpchn)
% delete Y (time) channels
%
% Usage:
% NewData = Dat.dumpY(Yi)
%
% Deletes data in the selected Y channel(s).
% Useful when specific detector channels are unreliable or unused.
%
% see also: dumpX, brushX
res = SP;
for i = 1:length(SP)
dumpchn = sort(dumpchn,'descend');
res(i).Z(dumpchn,:) = [];
res(i).Y(dumpchn) = [];
end
end %dumpY
function res = setXlim(SP,Xlim)
% set X-axis limits
%
% Usage:
% NewSpectrum = Dat.setXlim([Xmin Xmax])
%
% Sets minimum and maximum limits for the X values.
% Any datapoints lying outside the range specified by Xmin and Xmax
% are deleted.
%
% see also: binX
res = SP;
for i = 1:length(SP)
j = 1;
while (j <= res(i).Xn)
if (res(i).X(j) < Xlim(1)) || (res(i).X(j) > Xlim(2)),
res(i).X(j) = [];
res(i).Z(:,j) = [];
else
j = j + 1;
end %if
end %while
end %for
end %setXlim
function res = setYlim(SP,Ylim)
% set Y-axis limits
%
% Usage:
% NewDat = Dat.setYlim([Ymin Ymax])
%
% Sets minimum and maximum limits for the Y values.
% Any datapoints lying outside the range specified by Ymin and Ymax
% are deleted.
%
% see also: binY
res = SP;
for i = 1:length(SP)
j = 1;
while (j <= res(i).Yn)
if (res(i).Y(j) < Ylim(1)) || (res(i).Y(j) > Ylim(2)),
res(i).Y(j) = [];
res(i).Z(j,:) = [];
else
j = j + 1;
end %if
end %while
end %for
end %setYlim
function res = binX(SP,step)
% bin data across X channels
%
% Usage:
% NewData = Dat.binX(step)
%
% Bins (aggregates) datapoints using the specified X channel step.
% For example if Dat has X datapoints at every 0.1 nm,
% Dat.binX(10) will return a dataset where spectra have 1 nm resolution
% and each new datapoint is an average of 10 original datapoints.
% The binX function is useful to reduce the file size and noise %
% in the data but is also decreases spectral resolution.
%
% see also: setXlim, binY
res = SP;
for i = 1:length(SP)
bins = fix(SP(i).Xn/step);
res(i).X = zeros(1,bins);
res(i).Z = zeros(SP(i).Yn,bins);
for j = 1:bins
ia = (j-1)*step + 1;
ib = j*step - 1;
res(i).X(j) = mean(SP(i).X(ia:ib));
res(i).Z(:,j) = mean(SP(i).Z(:,ia:ib),2);
end %for
end %for
end %binZ
function res = binY(SP,step)
% bin data across Y channels
%
% Usage:
% NewData = Dat.binY(step)
%
% Bins (aggregates) datapoints using the specified Y channel step.
% The binY function is useful to reduce the file size and noise
% in the data but it also decreases time resolution.
%
% see also: binX
res = SP;
for i = 1:length(SP)
bins = fix(SP(i).Yn/step);
res(i).Y = zeros(bins,1);
res(i).Z = zeros(bins,SP(i).Xn);
for j = 1:bins
ia = (j-1)*step + 1;
ib = j*step - 1;
res(i).Y(j) = mean(SP(i).Y(ia:ib));
res(i).Z(j,:) = mean(SP(i).Z(ia:ib,:),1);
end %for
end %for
end %binY
function res = average(SP)
% create average dataset
%
% Usage:
% NewData = Dat.average
% NewData = average(Dat1, Dat2, )
%
% Creates a new dataset which is an average of all datasets
% in the object Dat (or the objects Dat1, Dat2, ).
% Datasets can only be averaged if their X and Y axes are identical.
if (length(SP) < 2)
fprintf('Warning! Number of datasets is less than two. Nothing to average.\n');
res = SP; return;
else
fprintf('Averaging %g datasets... ',length(SP));
Zsum = SP(1).Z;
for i = 2:length(SP)
if (~isequal(SP(1).X,SP(i).X)), error('X values do not match'); end
if (~isequal(SP(1).Y,SP(i).Y)), error('Y values do not match'); end
if (~isequal(size(SP(1)),size(SP(2)))), error('data dimensions do not match'); end
Zsum = Zsum + SP(i).Z;
end %for
res = SP(1);
res.Z = Zsum ./ length(SP);
end %if
fprintf('Done.\n');
end %average
function res = baselineY(SP,Xchn)
% subtract transient baseline
%
% Usage:
% NewData = Dat.baselineY(Xi)
%
% Subtracts baseline from the time traces.
% The baseline time trace is an average the time traces at wavelength channels Xi.
% The transient baseline can be taken at wavelengths where the signal %
% (absorbance/fluorescence/etc) from the sample is negligible.
%
% see also: baselineX
res = SP;
for i = 1:length(SP)
baseoff = mean(SP(i).Z(:,Xchn),2);
for j = 1:SP(i).Xn
res(i).Z(:,j) = res(i).Z(:,j) - baseoff;
end
end %for
end %baselineY
function res = baselineX(SP,Ychn)
% subtract spectral baseline
%
% Usage:
% NewData = Dat.baselineX(Yi)
%
% Subtracts baseline from the spectra.
% The baseline spectrum is the average of the spectra at Y channels Yi.
% Usually the data at times pre zero contain the spectral baseline.
%
% see also: baselineY
res = SP;
for i = 1:length(SP)
baseoff = mean(SP(i).Z(Ychn,:),1);
for j = 1:SP(i).Yn
res(i).Z(j,:) = res(i).Z(j,:) - baseoff;
end
end %for
end %baselineY
%% Save
function saveTIMP(SP,filename)
% Save spectra as ASCII file formatted for TIMP (wavelength-explicit).
%
% Usage:
% dat.save(filename)
%
% If dat contains more than one dataset, they are saved as separate files.
% filename can be a single string or a cell array of strings:
% {file1, file2, } for each dataset.
% If a single file name is given, numbers will be appended to the filename
% to identify each dataset.
%
% To save the data as a MATLAB file (.mat) instead,
% use the MATLAB built-in save command:
% save filename dat
for i=1:length(SP)
s = size(SP(i).Z,2);
if(nargin>=2)
if(iscell(filename))
fout = filename{i};
else
if (i > 1)
fout = strcat(filename,sprintf('%g.ascii',i));
else
fout = filename;
end
end
else
fout = strcat(SP(i).ID,'.ascii');
end
fprintf('Saving %s... ',fout);
fid = fopen(fout,'wt');
fprintf(fid,'%s\n',SP(i).ID);
fprintf(fid,'%s\n',SP(i).Comment);
fprintf(fid,'%s\n','Wavelength explicit');
fprintf(fid,'Intervalnr %d\n\t',s);
fprintf(fid,'%0g\t',SP(i).X(1:s-1));
fprintf(fid,'%0g\n',SP(i).X(s));
for j = 1:length(SP(i).Y)
fprintf(fid,'%0g\t',SP(i).Y(j));
fprintf(fid,'%0g\t',SP(i).Z(j,1:s-1));
fprintf(fid,'%0g\n',SP(i).Z(j,s));
end
fclose(fid);
fprintf('Done.\n');
end
end
%% Plotting
function plotSurface(SP,viewmode)
% plot data as a 2D colour map
if(length(SP)>1)
fprintf('Warning! Only first dataset will be plotted.');
SP = SP(1);
end
Cmax = max(abs(SP.Z(:)))/1.5;
cl = length(colormap)/2;
C = fix((SP.Z ./ Cmax + 1) .* cl);
surf(SP.X,SP.Y,SP.Z,C,'CDataMapping','direct');
if (SP.X(end)<SP.X(1)), set(gca,'Xdir','rev'); end
%set(gca,'YScale','log');
shading interp;
if (exist('viewmode')), view(viewmode), else view(2), end;
xlabel(sprintf('%s (%s)',SP.XType, SP.XUnit));
ylabel(sprintf('%s (%s)',SP.YType, SP.YUnit));
end %surface
function plotSpectra(SP,SelY,varargin)
% plot spectra at selected time points
%
% Usage:
% dat.plotSpectra(Y,Options)
%
% Time-resolved spectra are plotted for each time point in the array Y.
% If exact match is not found, the next larger time point is plotted.
%
% Optional parameters:
% 'Smooth' or 'Smooth','yes' - spectra are plotted as smooth lines
% 'Smooth', 'no' - spectra are plotted as straight lines (default)
%
% see also: plotTransients
if (length(SP)>1), SP = SP(1); end;
% Optional parameters
ShowExpID = false;
NoTeX = false;
SmoothL = false;
if nargin > 2,
for i = 1:nargin-2
switch(lower(varargin{i}))
case 'smooth'
% if 'smooth' is the last argument then
% assume smooth = yes
if (nargin>(i+2)),
smootharg = lower(varargin{i+1});
else
smootharg = 'yes';
end
switch smootharg
case 'no'
SmoothL = false;
case 'false'
SmoothL = false;
otherwise
SmoothL = true;
end
end %switch
end %for
end %if
if (~isempty(SelY))
for i = 1:length(SelY)
fy = find(SP.Y >= SelY(i),1);
if(isempty(fy)), error 'The selected Y values are not found in the dataset.'; end
if (SmoothL == true),
[plX,plY] = smoothLine(SP.X,SP.Z(fy,:));
else
plX = SP.X; plY = SP.Z(fy,:);
end
plot(plX,plY);
hold all
l{i} = sprintf('%g %s', SelY(i),SP.YUnit);
end
if (SP.X(end)<SP.X(1)), set(gca,'Xdir','rev'); end
legend(l,'location','NorthEast','FontSize',11);
legend boxoff
xlabel(sprintf('%s (%s)',SP.XType, SP.XUnit));
ylabel(sprintf('%s (%s)',SP.ZType, SP.ZUnit));
end
end %spectra
function plotTransients(SP,SelX)
% plot time traces at selected wavelengths
%
% Usage:
% dat.plotSpectra(X)
%
% Time traces are plotted for each wavelength in the array X.
% If exact wavelength match is not found, the next larger
% wavelength is plotted.
%
% see also: plotSpectra
if (length(SP)>1), SP = SP(1); end;
for i = 1:length(SelX)
plot(SP.Y,SP.Z(:,find(SP.X >= SelX(i),1),:));
hold all
l{i} = sprintf('%g %s', SelX(i),SP.XUnit);
end
legend(l,'location','NorthEast','FontSize',11);
xlabel(sprintf('%s (%s)',SP.YType, SP.YUnit));
ylabel(sprintf('%s (%s)',SP.ZType, SP.ZUnit));
end %Transients
end %methods
methods (Static)
function spect = load(filepattern,options, firstrecord, lastrecord)
% load data from ASCII text file(s)
%
% Usage:
% dat = TRData.load(filepattern, options)
%
% The files must contain only numerical data arranged in columns.
% The first column of the data must contain time points
% The following columns contain the measured data,
% where each column represents a wavelength channel.
%
% If time points occur several times (e.q. in repeated scan
% measurements), the data are averaged.
% filepattern can specify one or more filenames (e.g. *.DAT).
% If several files match the criteria, then dat is a cell
% array of datasets dat(1), dat(2), - for each file.
% options a structure specifying optional parameters (as fields):
% XType, XUnit, ExpID, Comment.
% The X axis also can be specified here,
% otherwise it will contain consecutive number
% (X values are not read from the file).
% IDs are not specified here but will carry the file names.
if (iscell(filepattern)),
files = filepattern;
else
files = getdir(filepattern);
end %if
if (~iscell(files)) , error('No files selected.'); end
if (length(files)<1), error('No files selected.'); end
fn = length(files);
%data = cell(1,fn);
for i = 1:fn
fprintf('Loading file %s... ',files{i});
store = load(files{i});
% check number of time points
timepoints = unique(store(:,1));
tn = length(timepoints);
cn = size(store,2); % number of columns
rn1 = size(store,1)/tn; % number of records (scans)
rn = fix(rn1);
if (rn1 > rn)
fprintf('Warning: Partial record will be ignored.\n');
end
if (rn < 1),
error('No valid records found in file.');
end
if (rn > 1),
% average if many records
jstart = 1;
jend = rn;
if (nargin>=3),
jstart = firstrecord;
end
if (nargin>=4),
jend = lastrecord;
end;
idata = zeros(tn,cn);
idata(:,1) = store(1:tn,1);
for j = jstart:jend
idata(:,2:cn) = idata(:,2:cn) + store(j*tn-tn+1:j*tn,2:cn);
end;
idata(:,2:cn) = idata(:,2:cn)/rn;
data = idata;
else
% only one record
data = store;
end;
fprintf('Done, %g record(s) read, %g time points.\n',rn, tn);
cols = size(data,2);
Y = data(:,1);
Z = data(:,2:cols);
X = 1:cols-1;
spect(i) = TRData(X,Y,Z,files{i});
if (nargin < 2),
options = struct;
end
if (isfield(options,'XType')),
spect(i).XType = options.XType;
end
if (isfield(options,'YType')),
spect(i).YType = options.YType;
end
if (isfield(options,'ZType')),
spect(i).ZType = options.ZType;
end
if (isfield(options,'XUnit')),
spect(i).XUnit = options.XUnit;
end
if (isfield(options,'YUnit')),
spect(i).YUnit = options.YUnit;
end
if (isfield(options,'ZUnit')),
spect(i).ZUnit = options.ZUnit;
end
if (isfield(options,'Comment')),
spect(i).Comment = options.Comment;
end
if (isfield(options,'ExpID')),
spect(i).ExpID = options.ExpID;
end
if (isfield(options,'X')),
if (~isequal(size(options.X),size(spect(i).X))),
error('Number of X values does not match the dataset.');
end
spect(i).X = options.X;
else
spect(i).XUnit = '';
spect(i).XType = 'Channel No.';
end
end
end %loadData
end %Static
end
