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