Upload files to ''

SyntheticLFPtoCSD.m

@@ -0,0 +1,179 @@
+clear all;
+close all;
+clc;
+
+%% LFP EXTRACTION SETTINGS
+% Time periods to extract from the data
+LFP_START_TIME = -1100;
+LFP_END_TIME   = 1100;
+% Account for delays between NLX events and actual visual presentations.
+%STIM_DELAY = 30.5 * 10^3;
+SelectedChannels=[1:16];
+NumOfChs=length(SelectedChannels);
+
+%% LFP GENERATION SETTINGS
+lfpNumSamples = ceil((LFP_END_TIME - LFP_START_TIME));
+LfpTimestamps = linspace(LFP_START_TIME, LFP_END_TIME, lfpNumSamples);
+
+% Main Peaks 
+PeaksLocations=[7 8 3 2 4 11]; % channel number where you expect to find peak
+PeaksAmplitudes=[-13 -11 40 28 32 -8]; % signal peak amplitude (not visible in z-scored plots)
+PeaksDelays=[50 45 65 55 70 59]; % [ms] where you expect to find peaks
+PeaksDurations=[10 4 12 8 4 8];
+% Side Peaks Randomly set
+SidePeaksAmplitudeStDev=2; % Std Dev of integer (positive/negative) non-zero values
+SidePeaksMaxDelays=900;    % Support of integer random shift (ranging from 0 to MaxDelays)
+SidePeaksDurations=4*ones(1,NumOfChs); % Duration of side peaks
+% Additive White Gaussian Noise Power set accordingly to
+AverageSNR=10^(25/10); % Average SNR you expect to see
+NumOfTrials=1e2;
+
+%% CSD EXTRACTION SETTINGS
+CSD_START_BASELINE_TIME = -200;
+CSD_END_BASELINE_TIME = 30;
+CSD_START_TIME = 0;
+CSD_END_TIME = 150;
+
+%% CSD GENERATION SETTINGS
+contactSpacing=1.5e-4;
+contactPositions = (1:NumOfChs).*contactSpacing;
+conductance = 0.4; % Taken from Logothetis NK, Kayser C, Oeltermann A. In vivo measurement of cortical impedance spectrum in monkeys: implications for signal propagation. Neuron 55: 809?823, 2007.
+diameter = 10e-4; % Taken from Mountcastle 1957
+filterRange = 5e-4; % Taken from Petterson et al Journal of Neurosci Methods 2006
+gaussSigma = 1e-4; % Taken from Petterson et al Journal of Neurosci Methods 2006
+
+csdStartIndex = find(LfpTimestamps >= CSD_START_TIME, 1, 'first');
+csdEndIndex = find(LfpTimestamps >= CSD_END_TIME, 1, 'first');
+csdZeroIndex = find(LfpTimestamps >= 0, 1, 'first');
+CSDNumTimeSamples = csdEndIndex - csdStartIndex + 1;
+CSDTimestamps = linspace(CSD_START_TIME, CSD_END_TIME, CSDNumTimeSamples);
+
+csdStartBaselineIndex = find(LfpTimestamps >= CSD_START_BASELINE_TIME, 1, 'first');
+csdEndBaselineIndex = find(LfpTimestamps >= CSD_END_BASELINE_TIME, 1, 'first');
+
+
+%% GENERATING SYNTHETIC LFPs
+
+
+LfpMat = zeros(NumOfChs, NumOfTrials, lfpNumSamples);
+
+% Generate Clean Signal
+%KernelSignal=1/sqrt(2*pi*PeakDuration^2)*exp(-(LfpTimestamps-PeakDelay).^2/(2*PeakDuration^2));
+PeaksDurationsVec=SidePeaksDurations; 
+PeaksDurationsVec(PeaksLocations)=PeaksDurations;
+CleanSignal=exp(-bsxfun(@rdivide, (LfpTimestamps).^2,(2*PeaksDurationsVec.^2)'));
+
+for tr=1:NumOfTrials
+    
+    % Randomly Modulate Side Peaks
+    % PeakAmps=randi(SidePeaksAmplitude,1,NumOfChs);
+    PeakAmps=zeros(1,NumOfChs);
+    while ~isempty(find(PeakAmps==0,1))
+        Ppos=find(PeakAmps==0,1);
+        PeakAmps(Ppos)=ceil(SidePeaksAmplitudeStDev*randn);
+    end
+    PeakAmps(PeaksLocations)=PeaksAmplitudes;
+    ModulatedSignal=diag(PeakAmps)*CleanSignal;
+
+    TimeShiftedSignal=ModulatedSignal;
+    % Delay the peaks
+    for ch=PeaksLocations
+        TimeShiftedSignal(ch,:)=circshift(ModulatedSignal(ch,:),[1 PeaksDelays(PeaksLocations==ch)]);
+    end
+    
+    % Randomly Scramble the timing of Side Peaks
+    SidePeaksToRandomlyShift=1:NumOfChs; 
+    SidePeaksToRandomlyShift(PeaksLocations)=[];
+    for ch=SidePeaksToRandomlyShift;
+        TimeShiftedSignal(ch,:)=circshift(ModulatedSignal(ch,:),[1 randi(SidePeaksMaxDelays)]);                
+    end
+    
+    % Add NOISE with same SNR for all channels, wrt their signal power
+    SigPow=sum(abs(TimeShiftedSignal).^2,2);
+    LfpNoise=diag(sqrt(SigPow/AverageSNR*1/lfpNumSamples))*randn(NumOfChs,lfpNumSamples);
+    NoisePow=sum(abs(LfpNoise).^2,2);
+    % Uncomment to check that average SNR is properly set
+    % EmpiricalSNR(:,tr,:)=SigPow./NoisePow;
+    LfpMat(:,tr,:)=TimeShiftedSignal+LfpNoise;
+end
+
+
+%% COMPUTE LFP TRIAL AVERAGE
+
+        LfpTimestampsCsdTimeWindow=LfpTimestamps(csdStartIndex:csdEndIndex);
+        
+        TotalBaseline=mean(LfpMat(:,:,csdStartBaselineIndex:csdEndBaselineIndex),3);
+        LfpMinusBaseline=LfpMat(:,:,csdStartBaselineIndex:csdEndBaselineIndex)-repmat(TotalBaseline,1,1,length(csdStartBaselineIndex:csdEndBaselineIndex));
+        TotalStdDev=std(LfpMinusBaseline,0,3);
+        DivStdDevMinusBaseline=(LfpMat(:,:,csdStartIndex:csdEndIndex)-repmat(TotalBaseline,1,1,CSDNumTimeSamples))./repmat(TotalStdDev+eps,1,1,CSDNumTimeSamples);
+        LfpTrialAvgCsdTimeWindow=squeeze(mean(DivStdDevMinusBaseline,2));
+        % Uncomment to see the true trial average signal non z-scored
+        % LfpTrialAvgCsdTimeWindow=squeeze(mean(LfpMat(:,:,csdStartIndex:csdEndIndex),2));
+        
+        EnergiesVec=sum(abs(LfpTrialAvgCsdTimeWindow),2);
+        ScaleFactorsVec=max(EnergiesVec)./EnergiesVec;
+        
+        LfpTrialAvgCsdTimeWindowNormalized=diag(ScaleFactorsVec)*LfpTrialAvgCsdTimeWindow;
+        
+        hFig4=figure('visible','on');
+        set(hFig4, 'Position', [0 0 1000 750]); cmp=jet(16);
+        for chf=1:16; subplot(2,2,1); plot(LfpTimestampsCsdTimeWindow,LfpTrialAvgCsdTimeWindow(chf,:),'Color',cmp(chf,:)); hold on; end;
+        axFig4=gca; axFig4.FontSize=8; axFig4.TickLabelInterpreter='latex'; %axFig4.Color=.3*ones(1,3); axFig4.GridColor=.7*ones(1,3);
+        xlabel('Time aligned to stimulus onset [ms]','Interpreter','latex'); grid on; grid minor;
+        xlim([LfpTimestampsCsdTimeWindow(1) LfpTimestampsCsdTimeWindow(end)]);
+        legend(arrayfun(@(x) num2str(x,'ch %.0f'),1:length(SelectedChannels),'unif',0),'Position',axFig4.Position*[1 0 1 0; 0 1 0 0; 0 0 0 0; 0 0 0 1]'+[0.032 0 0.01 0],'interpreter','latex');
+        title('Baseline z-scored LFPs','Interpreter','latex'); 
+        
+        for chf=1:16; subplot(2,2,2); plot(LfpTimestampsCsdTimeWindow,LfpTrialAvgCsdTimeWindowNormalized(chf,:),'Color',cmp(chf,:)); hold on; end;
+        axFig4=gca; axFig4.FontSize=8; axFig4.TickLabelInterpreter='latex'; %axFig4.Color=.3*ones(1,3); axFig4.GridColor=.7*ones(1,3);
+        xlabel('Time aligned to stimulus onset [ms]','Interpreter','latex'); grid on; grid minor;
+        xlim([LfpTimestampsCsdTimeWindow(1) LfpTimestampsCsdTimeWindow(end)]);
+        legend(arrayfun(@(x) num2str(x,'ch %.0f'),1:length(SelectedChannels),'unif',0),'Position',axFig4.Position*[1 0 1 0; 0 1 0 0; 0 0 0 0; 0 0 0 1]'+[0.032 0 0.01 0],'interpreter','latex');
+        title('Energy Normalized LFPs','Interpreter','latex');
+        
+        
+%% COMPUTE CSD
+    
+        CSDNumSpaceSamples=2e2;
+        [iCSDSamples,iCSDPositions]=computeiCSDCubicSplines(LfpTrialAvgCsdTimeWindow,contactPositions,conductance,diameter,CSDNumSpaceSamples);
+        
+        [iCSDSamplesNormalized,iCSDPositions]=computeiCSDCubicSplines(LfpTrialAvgCsdTimeWindowNormalized,contactPositions,conductance,diameter,CSDNumSpaceSamples);
+        
+        
+        set(0, 'currentfigure', hFig4); subplot(2,2,3);
+        surf(repmat(CSDTimestamps,length(iCSDPositions),1),repmat(iCSDPositions'*1e3,1,CSDNumTimeSamples),iCSDSamples*1e-3,'EdgeColor','none'); view(2);
+        axis([CSDTimestamps(1) CSDTimestamps(end) iCSDPositions(1)*1e3 iCSDPositions(end)*1e3]);
+        xlabel('Time aligned to stimulus onset [ms]','Interpreter','latex'); ylabel('Distance from the tip [mm]','Interpreter','Latex');
+        axFig4=gca; axFig4.YTick=(contactPositions)*1e3;
+        colormap(flipud(jet(256))); colorbar('fontsize',7,'Position',axFig4.Position*[1 0 1 0; 0 1 0 0; 0 0 0 0; 0 0 0 1]'+[0.032 0 0.01 0]);
+        axFig4.YTickLabel=arrayfun(@(x) num2str(x,'%.2f'),contactPositions*1e3,'unif',0);
+        yyaxis right; axFig4r=gca; axFig4r.YColor='k';
+        axFig4r.YLim=[iCSDPositions(1) iCSDPositions(end)]*1e3; axFig4r.YTick=(contactPositions)*1e3;
+        axFig4r.YTickLabel=arrayfun(@(x) num2str(x,'%.0f'),1:length(SelectedChannels),'unif',0);
+        aylp=axFig4r.YLabel.Position; ylabel('Channel Number','Rotation',270,'Interpreter','latex','Position',aylp);
+        axFig4r.TickLabelInterpreter='latex'; axFig4.FontSize=8;
+        yyaxis left;
+        title(['Energy Normalized iCSD [$\mu$A/mm$^3$]'] ,'Interpreter','latex');
+        
+        
+        set(0, 'currentfigure', hFig4); subplot(2,2,4);
+        surf(repmat(CSDTimestamps,length(iCSDPositions),1),repmat(iCSDPositions'*1e3,1,CSDNumTimeSamples),iCSDSamplesNormalized*1e-3,'EdgeColor','none'); view(2);
+        axis([CSDTimestamps(1) CSDTimestamps(end) iCSDPositions(1)*1e3 iCSDPositions(end)*1e3]);
+        xlabel('Time aligned to stimulus onset [ms]','Interpreter','latex'); ylabel('Distance from the tip [mm]','Interpreter','Latex');
+        axFig4=gca; axFig4.YTick=(contactPositions)*1e3;
+        colormap(flipud(jet(256))); colorbar('fontsize',7,'Position',axFig4.Position*[1 0 1 0; 0 1 0 0; 0 0 0 0; 0 0 0 1]'+[0.032 0 0.01 0]);
+        axFig4.YTickLabel=arrayfun(@(x) num2str(x,'%.2f'),contactPositions*1e3,'unif',0);
+        yyaxis right; axFig4r=gca; axFig4r.YColor='k';
+        axFig4r.YLim=[iCSDPositions(1) iCSDPositions(end)]*1e3; axFig4r.YTick=(contactPositions)*1e3;
+        axFig4r.YTickLabel=arrayfun(@(x) num2str(x,'%.0f'),1:length(SelectedChannels),'unif',0);
+        aylp=axFig4r.YLabel.Position; ylabel('Channel Number','Rotation',270,'Interpreter','latex','Position',aylp);
+        axFig4r.TickLabelInterpreter='latex'; axFig4.FontSize=8;
+        yyaxis left;
+        title(['Energy Normalized iCSD [$\mu$A/mm$^3$]'] ,'Interpreter','latex');
+        
+        ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 1],'Box','off','Visible','off','Units','normalized', 'clipping' , 'off');
+        text(0.5,.98,['LFP/CSD Cortex Visualization'],'HorizontalAlignment','center','VerticalAlignment','top','Interpreter','latex');
+        editedFF='./V1 Jones/Spline iCSD Monitor/Minus Baseline Divided StDev Energy Norm/Comparison/';
+        hFig4.InvertHardcopy='off';
+        saveas(hFig4, ['PLOT_iCSD.png']);
+        saveas(hFig4, ['PLOT_iCSD.fig']);

computeiCSDCubicSplines.m

@@ -0,0 +1,134 @@
+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