function signal = readref (filename, linestoignore)
% signal = readref(filename [, linestoignore])
% Reads the ref file passed in and returns the signal; linestoignore
% defaults to 10, but can be custom set.
if nargin < 2, linestoignore = 10; end
if nargin < 1, error('No filename given to readref().'); end
if ~exist(filename, 'file'), error([filename ' is not a file.']); end
try
fid = fopen(filename);
ignore = textscan(fid, '%s', linestoignore, 'delimiter', '\n');
[ data ] = textscan(fid,'%n');
signal = data{1};
catch e
fclose(fid);
rethrow(e);
end
fclose(fid);
end
function [timeseries_rh timeseries_lh] = extract_timeseries(dataset, subjectID) % % Extract timeseries from .ref file (VoxBo format) % _______________________________________________ % by Marcelo G Mattar (05/24/2013) % mattar@sas.upenn.edu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% EXTRACT TIMESERIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% baseDir = '/jet/aguirre/FaceSpace/'; timeseriesFile_rh = [baseDir dataset '/subjects/' subjectID '/Models/' subjectID '_GLM_nointerests/rh-extracted_timeseries.ref']; timeseries_rh = readref (timeseriesFile_rh,2); timeseriesFile_lh = [baseDir dataset '/subjects/' subjectID '/Models/' subjectID '_GLM_nointerests/lh-extracted_timeseries.ref']; timeseries_lh = readref (timeseriesFile_lh,2);
function FiltVec = hrfconv(signal, SOA)
% FILTVEC = HRFCONV(SIGNAL)
% convolves a signal with a standard 10Hz-sampled HRF
% we assume a 1.5s trial duration and a 10Hz output signal
EndogenousFilter = [5.55112e-17;0.00767224;0.0162726;0.025849;0.0364457;0.0481024;0.0608536;0.0747281;0.0897482;0.10593;0.123281;0.141801;0.161484;0.182312;0.204261;0.227297;0.251376;0.276448;0.30245;0.329315;0.356963;0.38531;0.414261;0.443716;0.473569;0.503705;0.534007;0.564352;0.594614;0.624665;0.654374;0.68361;0.712242;0.740139;0.767176;0.793226;0.818169;0.841889;0.864277;0.885229;0.90465;0.922451;0.938554;0.952889;0.965396;0.976025;0.984737;0.991502;0.996304;0.999136;1;0.998913;0.995897;0.990989;0.984232;0.975681;0.965398;0.953452;0.939923;0.924892;0.908452;0.890696;0.871725;0.851641;0.830549;0.808556;0.78577;0.7623;0.738253;0.713734;0.688848;0.663695;0.638375;0.61298;0.5876;0.562319;0.537218;0.51237;0.487844;0.463702;0.44;0.41679;0.394115;0.372016;0.350523;0.329666;0.309467;0.289942;0.271104;0.25296;0.235516;0.218771;0.202721;0.187361;0.172683;0.158674;0.145324;0.132616;0.120536;0.109067;0.0981926;0.0878945;0.0781552;0.0689569;0.0602817;0.0521123;0.0444312;0.0372216;0.0304668;0.0241507;0.0182571;0.0127706;0.00767572;0.00295715;-0.00140033;-0.00541187;-0.00909271;-0.0124583;-0.0155242;-0.0183063;-0.0208208;-0.0230841;-0.025113;-0.0269245;-0.0285358;-0.0299644;-0.0312277;-0.0323428;-0.0333269;-0.0341965;-0.0349674;-0.0356548;-0.0362724;-0.0368327;-0.0373465;-0.0378227;-0.0382678;-0.038686;-0.0390787;-0.0394442;-0.0397777;-0.0400708;-0.0403115;-0.0404842;-0.0405693;-0.0405433;-0.0403786;-0.0400439;-0.039504;-0.0387199;-0.0376493;-0.0362464;-0.0344627;-0.0322468;-0.0295455;-0.0263035;-0.0224645;-0.0179715;-0.0127673;-0.00679528];
EndoLength=length(signal)*(SOA/100);
EndogenousFilter=[EndogenousFilter; zeros(EndoLength-length(EndogenousFilter),1)];
EndogenousFilter=NormMag( (EndogenousFilter));
EndogenousKernel=fft(EndogenousFilter);
EndogenousKernel(1)=1; % set DC component to pin amplitude
signal=mncnleavezeros(signal);
%signal=signal - mean(signal);
Vec=Geoffcongrid(signal,EndoLength);
FiltVec=ifft(fft(Vec).*EndogenousKernel);
function NewImage = NormMag(Image)
%
% 2007 ddrucker@psych.upenn.edu adapted from gka
FTImage=fft(Image);
MagImage=sqrt(FTImage.*conj(FTImage));
ZeroPoints=find(MagImage==0);
ZeroCount=length(ZeroPoints);
if ZeroCount > 0
MagImage(ZeroPoints)=1;
end
PhaseImage=acos(real(FTImage./MagImage));
NegVals=find(imag(FTImage) < 0);
NegCount=length(NegVals);
if NegCount > 0
PhaseImage(NegVals) = 2*pi - PhaseImage(NegVals);
end
if ZeroCount > 0
PhaseImage(ZeroPoints) = 0;
MagImage(ZeroPoints) = 0;
end
MagImage=MagImage./max(MagImage);
RealPart=MagImage.*cos(PhaseImage);
ImaginaryPart=MagImage.*sin(PhaseImage);
FTImage=complex(RealPart,ImaginaryPart);
NewImage=real(ifft(FTImage));
function meancentered = mncnleavezeros(vec)
% MEANCENTERED = MNCNLEAVEZEROS(VEC)
% mean-centers a vector (or matrix), ignoring 0 values (0s stay 0 and are not included
% as part of the calculation of the mean)
%
% 2007 ddrucker@psych.upenn.edu
nonzeros = find(vec~=0);
nz = numel(nonzeros);
sig = sum(vec(nonzeros));
mean = sig/nz;
meancentered = vec;
meancentered(nonzeros) = meancentered(nonzeros)-mean;
function NewSignal = Geoffcongrid(Signal,NewSize)
% note this only works with vertically oriented ones
%
% 2007 ddrucker@psych.upenn.edu
OrigLength=size(Signal,1);
Ratio=NewSize/OrigLength;
if Ratio ~= floor(Ratio)
error('ratio is not integer');
end
nSimMats=size(Signal,2);
NewSignal=zeros(NewSize,nSimMats);
for i=1:Ratio
NewSignal(i:Ratio:Ratio*(OrigLength-1)+i,:) = Signal;
end
%**************************************************************************
%
% This is a simple extension of the bar plot to include error bars. It
% is called in exactly the same way as bar but with an extra input
% parameter "errors" passed first.
%
% Parameters:
% errors - the errors to be plotted (extra dimension used if assymetric)
% varargin - parameters as passed to conventional bar plot
% See bar and errorbar documentation for more details.
%
% Output:
% H = barwitherr(..) returns a vector of handles to the barseries objects
%
% Symmetric Example:
% y = randn(3,4); % random y values (3 groups of 4 parameters)
% errY = 0.1.*y; % 10% error
% h = barwitherr(errY, y);% Plot with errorbars
%
% set(gca,'XTickLabel',{'Group A','Group B','Group C'})
% legend('Parameter 1','Parameter 2','Parameter 3','Parameter 4')
% ylabel('Y Value')
% set(h(1),'FaceColor','k');
%
%
% Asymmetric Example:
% y = randn(3,4); % random y values (3 groups of 4 parameters)
% errY = zeros(3,4,2);
% errY(:,:,1) = 0.1.*y; % 10% lower error
% errY(:,:,2) = 0.2.*y; % 20% upper error
% barwitherr(errY, y); % Plot with errorbars
%
% set(gca,'XTickLabel',{'Group A','Group B','Group C'})
% legend('Parameter 1','Parameter 2','Parameter 3','Parameter 4')
% ylabel('Y Value')
%
%
% Notes:
% Ideally used for group plots with non-overlapping bars because it
% will always plot in bar centre (so can look odd for over-lapping bars)
% and for stacked plots the errorbars will be at the original y value is
% not the stacked value so again odd appearance as is.
%
% The data may not be in ascending order. Only an issue if x-values are
% passed to the fn in which case their order must be determined to
% correctly position the errorbars.
%
%
% 24/02/2011 Martina F. Callaghan Created
% 12/08/2011 Martina F. Callaghan Updated for random x-values
% 24/10/2011 Martina F. Callaghan Updated for asymmetric errors
% 15/11/2011 Martina F. Callaghan Fixed bug for assymetric errors &
% vector plots
% 14/06/2013 Martina F. Callaghan Returning handle as recommended by
% Eric (see submission comments)
%
%**************************************************************************
function h = barwitherr(errors,varargin)
% Check how the function has been called based on requirements for "bar"
if nargin < 3
% This is the same as calling bar(y)
values = varargin{1};
xOrder = 1:size(values,1);
else
% This means extra parameters have been specified
if isscalar(varargin{2}) || ischar(varargin{2})
% It is a width / property so the y values are still varargin{1}
values = varargin{1};
xOrder = 1:size(values,1);
else
% x-values have been specified so the y values are varargin{2}
% If x-values have been specified, they could be in a random order,
% get their indices in ascending order for use with the bar
% locations which will be in ascending order:
values = varargin{2};
[tmp xOrder] = sort(varargin{1});
end
end
% If an extra dimension is supplied for the errors then they are
% assymetric split out into upper and lower:
if ndims(errors) == ndims(values)+1
lowerErrors = errors(:,:,1);
upperErrors = errors(:,:,2);
elseif isvector(values)~=isvector(errors)
lowerErrors = errors(:,1);
upperErrors = errors(:,2);
else
lowerErrors = errors;
upperErrors = errors;
end
% Check that the size of "errors" corresponsds to the size of the y-values.
% Arbitrarily using lower errors as indicative.
if any(size(values) ~= size(lowerErrors))
error('The values and errors have to be the same length')
end
[nRows nCols] = size(values);
handles.bar = bar(varargin{:}); % standard implementation of bar fn
hold on
h = handles.bar;
if nRows > 1
for col = 1:nCols
% Extract the x location data needed for the errorbar plots:
x = get(get(handles.bar(col),'children'),'xdata');
% Use the mean x values to call the standard errorbar fn; the
% errorbars will now be centred on each bar; these are in ascending
% order so use xOrder to ensure y values and errors are too:
errorbar(mean(x,1),values(xOrder,col),lowerErrors(xOrder,col), upperErrors(xOrder, col), '.k')
end
else
x = get(get(handles.bar,'children'),'xdata');
errorbar(mean(x,1),values,errors,'.k')
end
hold off
function pout=myBinomTest(s,n,p,Sided)
%function pout=myBinomTest(s,n,p,Sided)
%
% Performs a binomial test of the number of successes given a total number
% of outcomes and a probability of success. Can be one or two-sided.
%
% Inputs:
% s- (Scalar or Array) The observed numebr of successful outcomes
% n- (Scalar or Array) The total number of outcomes (successful or not)
% p- (Scalar or Array) The proposed probability of a successful outcome
% Sided- (String) can be 'Two','Greater' or 'Lesser'. A value of
% 'Two' will perform a two-sided test, that the actual number
% of success is different from the expected number of
% successes in any direction. 'Greater' or 'Lesser' will
% perform a 1-sided test to examine if the observed number of
% successes are either significantly greter than or less than
% (respectively) the expected number of successes.
%
% Outputs:
% pout- The probability of observing the resulting value of s or a
% value more extreme (the precise meaning of which depends on
% the value of Sided) given n total outcomes with a
% probability of success of p.
%
% s, n and p can be scalars or arrays of the same size. The
% dimensions and size of pout will match that of these inputs.
%
% For example, the signtest is a special case of this where the value of p
% is equal to 0.5 (and a 'success' is dfeined by whether or not a given
% sample is of a particular sign.), but the binomial test and this code is
% more general allowing the value of p to be any value between 0 and 1.
%
% References:
% http://en.wikipedia.org/wiki/Binomial_test
%
% by Matthew Nelson July 21st, 2009
% matthew.j.nelson.vumail@gmail.com
if nargin<4 || isempty(Sided); Sided='Two'; end
if nargin<3 || isempty(p); p=0.5; end
switch lower(Sided)
case 'two'
[s,n,p]= EqArrayAndScalars(s,n,p);
E=p.*n;
GreaterInds=s>=E;
pout=repmat(0,size(GreaterInds));
dE=pout;
%note that matlab's binocdf(s,n,p) gives the prob. of getting up to AND INCLUDING s # of successes...
%Calc pout for GreaterInds first
pout(GreaterInds)=1-binocdf( s(GreaterInds)-1,n(GreaterInds),p(GreaterInds)); %start with the prob of getting >= s # of successes
%now figure the difference from the expected value, and figure the prob of getting lower than that difference from the expected value # of successes
dE(GreaterInds)=s(GreaterInds)-E(GreaterInds);
pout(GreaterInds)=pout(GreaterInds)+ binocdf(floor(E(GreaterInds)-dE(GreaterInds)),n(GreaterInds),p(GreaterInds)); %the binonmial is a discrete dist. ... so it's value over non-integer args has no menaing... this flooring of E-dE actually doesn't affect the outcome (the result is teh same if the floor was removed) but it's included here as a reminder of the discrete nature of the binomial
%If the expected value is exactly equaled, the above code would have added the probability at that discrete value twice, so we need to adjust (in this case, pout will always = 1 anyways)
if dE==0 pout(GreaterInds)=pout(GreaterInds)- binopdf( E(GreaterInds),n(GreaterInds),p(GreaterInds) ); end
%Calc pout for LesserInds second
pout(~GreaterInds)=binocdf(s(~GreaterInds),n(~GreaterInds),p(~GreaterInds)); %start with the prob of getting <= s # of successes
%now figure the difference from the expected value, and figure the prob of getting greater than that difference from teh expected value # of successes
dE(~GreaterInds)=E(~GreaterInds)-s(~GreaterInds);
pout(~GreaterInds)=pout(~GreaterInds) + 1-binocdf(ceil(E(~GreaterInds)+dE(~GreaterInds))-1,n(~GreaterInds),p(~GreaterInds)); %Here teh ceiling is needed b/c of teh -1 following it, so that integer and non-integer vals of E+dE will bothe give teh correct value with teh same line of code
case 'greater' %one-sided
pout=1-binocdf(s-1,n,p); %just report the prob of getting >= s # of successes
case 'lesser' %one-sided
pout=binocdf(s,n,p); %just report the prob of getting <= s # of successes
otherwise
error(['In myBinomTest, Sided variable is: ' Sided '. Unkown sided value.'])
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout=EqArrayAndScalars(varargin)
%function varargout=EqArrayAndScalars(varargin)
%
% This will compare a collection of inputs that must be either scalars or
% arrays of the same size. If there is at least one array input, all scalar
% inputs will be replicated to be the array of that same size. If there are
% two or more array inputs that have different sizes, this will return an
% error.
%
% created by Matthew Nelson on April 13th, 2010
% matthew.j.nelson.vumail@gmail.com
d=repmat(0,nargin,1);
for ia=1:nargin
d(ia)=ndims(varargin{ia});
end
maxnd=max(d);
s=repmat(1,nargin,maxnd);
for ia=1:nargin
s(ia,1:d(ia))=size(varargin{ia});
end
maxs=max(s);
for ia=1:nargin
if ~all(s(ia,:)==maxs)
if ~all(s(ia,:)==1)
error(['Varargin{' num2str(ia) '} needs to be a scalar or equal to the array size of other array inputs.'])
else
varargout{ia}=repmat(varargin{ia},maxs);
end
else
varargout{ia}=varargin{ia};
end
end