123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134 |
- function [iCSDMat,zAxis]=computeiCSDCubicSplines(LFPMat,cAxis,conductance,diameter,numCSDSpaceSamples,GaussFiltSigma,filtRange,integralTolerance)
- %
- % This function computes iCSD Cubic Splines (Pettersen et al., J Neurosci Methods, 2005)
- %
- % List of INPUTs:
- % LFPMat -> Matrix of LFP recordings (Probe Channels x Time Sequence) [V]
- % cAxis -> Vector of Probe contact points [m]
- % conductance -> Conductance of the Extracellular Fluid [S/m]
- % diameter -> Diameter of disk around the Probe [m]
- % numCSDSpaceSamples -> Number of output samples (Spatial resolution)
- % GaussFiltSigma -> Standard Deviation of the Gaussian Filter
- % filtRange -> Extremity of the Gaussian Filter
- % integralTolerance -> Convergence criteria of the numeric integral approximation
- %
- % List of OUTPUTs:
- % iCSDMat -> Matrix of iCSDs
- % zAxis -> Vector of CSD surface positions
- %
-
- if nargin<8; integralTolerance=1e-6;
- if nargin<7;
- if nargin<6; GaussFiltSigma=1e-4;
- if nargin<5; numCSDSpaceSamples=200; end;
- end;
- filtRange=5*GaussFiltSigma;
- end;
- end;
-
- % Probe Contacts Parameters
- [numContacts,numCSDTimeSamples] = size(LFPMat);
- cAvgDiff=mean(diff(cAxis));
- % Add imaginary contacts on the Border
- cAxisBpad=[cAxis(1)-cAvgDiff cAxis cAxis(end)+cAvgDiff];
- % Output axis with higher spatial resolution
- zAxis=linspace(cAxis(1)-cAvgDiff/2,cAxis(end)+cAvgDiff/2,numCSDSpaceSamples);
- zAvgDiff=mean(diff(zAxis));
- zAxisBpad=[cAxisBpad(1):zAvgDiff:zAxis(1) zAxis(2:end-1) zAxis(end):zAvgDiff:cAxisBpad(end)];
- zAxisIndicesWithinBpad=find(zAxisBpad<=zAxis(1),1,'last'):find(zAxisBpad>=zAxis(end),1,'first');
- effNumCSDSpaceSamples=length(zAxisBpad);
-
- % C Vector and C Matrix
- CVec = 1./diff(cAxisBpad); CMat=diag(CVec);
- % C Matrix zero padded in a Up Frame
- CMat0padUF=zeros(numContacts+2); CMat0padUF(2:end-1,2:end-1)=diag(CVec(1:end-1));
- % C Matrix zero padded in a Down Frame
- CMat0padDF=zeros(numContacts+2); CMat0padDF(2:end-1,2:end-1)=diag(CVec(2:end));
- % C Matrix ones added along Diagonal Margins
- CMat1addDM=CMat0padDF; CMat1addDM(1)=1; CMat1addDM(end)=1;
- % Identity Matrix zero padded Up
- IMat0padU=[zeros(numContacts+1,1) eye(numContacts+1)];
- % Identity Matrix zero padded Down
- IMat0padD=[eye(numContacts+1) zeros(numContacts+1,1)];
- % Identity Matrix zero padded Up Down Left
- IMat0padUDL=[zeros(1,numContacts+2); eye(numContacts+1) zeros(numContacts+1,1)];
- % Identity Matrix zero padded Up Down 2x Right
- IMat0padUD2R=[zeros(1,numContacts+2); zeros(numContacts,2) eye(numContacts); zeros(1,numContacts+2)];
- % Identity Matrix zero padded in a Frame
- IMat0padF=[zeros(1,numContacts+2); zeros(numContacts,1) eye(numContacts) zeros(numContacts,1); zeros(1,numContacts+2)];
- % K Matrix divergency of cubic splines
- KMat = (CMat0padUF*IMat0padUDL+2*CMat0padUF*IMat0padF+2*CMat1addDM+CMat0padDF*IMat0padUD2R)^(-1)*3*...
- (CMat0padUF^2*IMat0padF-CMat0padUF^2*IMat0padUDL+CMat0padDF^2*IMat0padUD2R-CMat0padDF^2*IMat0padF);
- % E Matrix coefficients for cubic spline
- E0th = IMat0padD;
- E1st = IMat0padD*KMat;
- E2nd = 3*CMat^2*(IMat0padU-IMat0padD)-CMat*(IMat0padU+2*IMat0padD)*KMat;
- E3rd = 2*CMat^3*(IMat0padD-IMat0padU)+CMat^2*(IMat0padU+IMat0padD)*KMat;
- clear CVec CMat CMat0padUF CMat0padDF CMat1addDM IMat0padU IMat0padD IMat0padDL IMat0padUD2R IMat0padUDL IMat0padF KMat
- % Define integration matrixes
- F0th = zeros(numContacts,numContacts+1);
- F1st = zeros(numContacts,numContacts+1);
- F2nd = zeros(numContacts,numContacts+1);
- F3rd = zeros(numContacts,numContacts+1);
- for jth = 1:numContacts % Note: Needs to be fixed if conductivity not constant (conductivity~=conductivityTop)!
- for ith = 1:numContacts+1
- F0th(jth,ith) = integral(@(zeta) 1./(2.*conductance).*(sqrt((diameter/2)^2+((cAxisBpad(jth+1)-zeta)).^2)-abs(cAxisBpad(jth+1)-zeta)),...
- cAxisBpad(ith),cAxisBpad(ith+1),'RelTol',0,'AbsTol',integralTolerance);
- F1st(jth,ith) = integral(@(zeta) 1./(2.*conductance).*(zeta-cAxisBpad(ith)).*(sqrt((diameter/2)^2+((cAxisBpad(jth+1)-zeta)).^2)-abs(cAxisBpad(jth+1)-zeta)),...
- cAxisBpad(ith),cAxisBpad(ith+1),'RelTol',0,'AbsTol',integralTolerance);
- F2nd(jth,ith) = integral(@(zeta) 1./(2.*conductance).*(zeta-cAxisBpad(ith)).^2.*(sqrt((diameter/2)^2+((cAxisBpad(jth+1)-zeta)).^2)-abs(cAxisBpad(jth+1)-zeta)),...
- cAxisBpad(ith),cAxisBpad(ith+1),'RelTol',0,'AbsTol',integralTolerance);
- F3rd(jth,ith) = integral(@(zeta) 1./(2.*conductance).*(zeta-cAxisBpad(ith)).^3.*(sqrt((diameter/2)^2+((cAxisBpad(jth+1)-zeta)).^2)-abs(cAxisBpad(jth+1)-zeta)),...
- cAxisBpad(ith),cAxisBpad(ith+1),'RelTol',0,'AbsTol',integralTolerance);
- end;
- end;
- FMat=zeros(numContacts+2); FMat(2:numContacts+1,:)=F0th*E0th+F1st*E1st+F2nd*E2nd+F3rd*E3rd; FMat(1)=1; FMat(end)=1;
- clear F0th F1st F2nd F3rd ith jth
- VMat0padUD=[zeros(1,numCSDTimeSamples); LFPMat; zeros(1,numCSDTimeSamples)];
- % CSDLowRes=F^(-1)*VMat0padUD;
- iCSDLowRes=FMat\VMat0padUD;
- % The cubic spline polynomial coeffecients
- A0th = E0th*iCSDLowRes;
- A1st = E1st*iCSDLowRes;
- A2nd = E2nd*iCSDLowRes;
- A3rd = E3rd*iCSDLowRes;
- % Preallocate Unfiltered CSD Matrix
- uiCSDMat = zeros(effNumCSDSpaceSamples,numCSDTimeSamples);
- for zth=1:effNumCSDSpaceSamples;
- iCAxisBpad=find([cAxisBpad(1:end-1) cAxisBpad(end)+eps]>zAxisBpad(zth),1,'first')-1;
- uiCSDMat(zth,:) = A0th(iCAxisBpad,:) + A1st(iCAxisBpad,:).*(zAxisBpad(zth)-cAxisBpad(iCAxisBpad)) + ...
- A2nd(iCAxisBpad,:)*(zAxisBpad(zth)-cAxisBpad(iCAxisBpad)).^2 + A3rd(iCAxisBpad,:)*(zAxisBpad(zth)-cAxisBpad(iCAxisBpad)).^3;
- end;
-
- clear E0th E1st E2nd E3rd A0th A1st A2nd A3rd VMat0padUD CSDLowRes F iCAxisBpad zth
- % Create Gaussian Filter
- filtAxes = -filtRange/2:zAvgDiff:filtRange/2;
- numFiltSamples = length(filtAxes);
- GaussFiltSamples = 1/(GaussFiltSigma*sqrt(2*pi))*exp(-filtAxes.^2/(2*GaussFiltSigma^2));
- % Zero pad CSD Matrix of filter size Up and Down
- uiCSDMat0padfUD = [zeros(numFiltSamples,numCSDTimeSamples); uiCSDMat; zeros(numFiltSamples,numCSDTimeSamples)];
- % Filter padded CSD Matrix
- iCSDMat0padfUD = filter(GaussFiltSamples/sum(GaussFiltSamples),1,uiCSDMat0padfUD);
- % Remove extra samples due to the Zero padding and Filtering function (which adds 1/2 filter size samples)
- iCSDMat=iCSDMat0padfUD(round(1.5*numFiltSamples)+1:round(1.5*numFiltSamples)+effNumCSDSpaceSamples,:);
- clear uiCSDMat iCSDMat0padfUD uiCSDMat0padfUD filtAxes GaussFiltSamples numFiltSamples
- % Recenter the iCSDMat to zAxis (remove imaginary contacts)
- iCSDMat=iCSDMat(zAxisIndicesWithinBpad,:);
- end
|