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