|
@@ -0,0 +1,115 @@
|
|
|
+function [naturalMovie_output]=BC_predNaturalMovies(loadData,filename,Hz,STA_output,NL_output)
|
|
|
+%[naturalMovie_output]=BC_predNaturalMovies(loadData,filename,pixelSize,Hz,STA_output,NL_output)
|
|
|
+%this function shows the prinicpal calculation for the prediction of the
|
|
|
+%response to the natural movie, from the linear-nonlinear model (LN-model)
|
|
|
+%computed with white noise stimulus.
|
|
|
+%Inputs:
|
|
|
+% loadData = the folder path of the natural movie files
|
|
|
+% filename = the filename of the BC of interest
|
|
|
+% movieData
|
|
|
+% Hz = stimulus update rate
|
|
|
+% STA_output = output of the BC_STA function; containing the
|
|
|
+% STA (filter) of the white noise stimulus
|
|
|
+% NL_output = output of the BC_NL function; containing the
|
|
|
+% output nonlinearity of the white noise stimulus
|
|
|
+%
|
|
|
+%Outpus: naturalMovie_output = structure containing the output of the
|
|
|
+% prediction
|
|
|
+%
|
|
|
+%
|
|
|
+
|
|
|
+%-------------------------------------------------------------------------------
|
|
|
+%% 1. Get the movie contrast values
|
|
|
+%the movies can be downloaded under the following
|
|
|
+%link:https://zenodo.org/record/46481. See documentation.data.pdf for
|
|
|
+%further instructions, also which frames to take.
|
|
|
+%they have to be prepared as follow:
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+%-------------------------------------------------------------------------------
|
|
|
+%% 2. compute the prediction
|
|
|
+%% 2.1 convolution
|
|
|
+%convolve the new stimulus sequence (not used for prediciton) with the
|
|
|
+%normalized STA
|
|
|
+stim2pred=filter2(STA_norm,stimPred,'valid');
|
|
|
+
|
|
|
+%% 2.2 compute the prediction through linear interpolation
|
|
|
+NL_x=bins.NL_xbin;
|
|
|
+NL_y=bins.NL_ybin;
|
|
|
+%example of how to make the prediction through linear interpolation
|
|
|
+predMemP_y=nan(1,size(stim2pred,2));
|
|
|
+for kk=1:size(stim2pred,2)
|
|
|
+ nbrTemp=stim2pred(kk); %first x-axis value to predict
|
|
|
+ [xIndxS,xIndxL]= getIndex(NL_x,nbrTemp,1); %get those values that are on left and rigth
|
|
|
+ xValueS=NL_x(xIndxS);
|
|
|
+ xValueL=NL_x(xIndxL);
|
|
|
+ %do the linear interpolation
|
|
|
+ [polyCoeff] = polyfit([xValueS;xValueL],[NL_y(xIndxS);...
|
|
|
+ NL_y(xIndxL)],1);
|
|
|
+ predMemP_y(kk)=polyCoeff(1)*nbrTemp+polyCoeff(2);
|
|
|
+end
|
|
|
+
|
|
|
+%-------------------------------------------------------------------------------
|
|
|
+%% 3. evaluate the prediction performance
|
|
|
+%pearson correlation
|
|
|
+corrcoeff_temp=corrcoef(voltPred(STA_output.STA_fnbr:end),predMemP_y); %the first value is at the length of the filter
|
|
|
+pearCorrSqrt=corrcoeff_temp(2)^2;
|
|
|
+
|
|
|
+%-------------------------------------------------------------------------------
|
|
|
+%% 4. plot the output
|
|
|
+figure;
|
|
|
+%STA
|
|
|
+subplot(2,2,1)
|
|
|
+[~,indx_Zoom]=max(abs(STA3D_norm(:))); %max intensity value of STA
|
|
|
+[~,~,posZ]=ind2sub(size(STA3D_norm),indx_Zoom); %get the frame position
|
|
|
+imagesc(STA3D_norm(:,:,posZ))
|
|
|
+colormap('gray'); %how the sitmulus was shown
|
|
|
+axis equal; axis tight;
|
|
|
+title('STA best frame')
|
|
|
+%NL
|
|
|
+subplot(2,2,2)
|
|
|
+plot(NL_xAxis,NL_yAxis,'.','color','k') %plot non-binned signal
|
|
|
+hold on;
|
|
|
+plot(bins.NL_xbin,bins.NL_ybin,'o','markerfacecolor','r','markeredgecolor','r') %plot binned signal
|
|
|
+title('NL')
|
|
|
+axis tight
|
|
|
+xlabel('generator signal')
|
|
|
+ylabel('voltage (mV)')
|
|
|
+%prediction
|
|
|
+subplot(2,2,3:4)
|
|
|
+plot(voltPred(STA_output.STA_fnbr:end),'k') %original membrane potential
|
|
|
+hold on;
|
|
|
+plot(predMemP_y,'r')%predicted membrane potential
|
|
|
+title(['pCorr:' mat2str(round(pearCorrSqrt,2)) ' Filename: ' STA_output.filename])
|
|
|
+xlabel('stimulus frames')
|
|
|
+ylabel('voltage (ms)')
|
|
|
+legend('original response','predicted response')
|
|
|
+
|
|
|
+%-------------------------------------------------------------------------------
|
|
|
+%% 5. add variables into the output
|
|
|
+pred_output.pearCorrSqrt=pearCorrSqrt;
|
|
|
+
|
|
|
+end
|
|
|
+
|
|
|
+function [IndS,IndL]= getIndex(x,valConv,Points)
|
|
|
+IndS=find(x<valConv,Points,'last'); %find the last 5 closest smaller values
|
|
|
+IndL=find(x>=valConv,Points);
|
|
|
+if isempty(IndS)|| length(IndS)<Points %if you are on the left border
|
|
|
+ if isempty(IndS) %take value from 1:10
|
|
|
+ IndL=find(x>=valConv,Points*2);
|
|
|
+ IndS=IndL(1);
|
|
|
+ IndL=IndL(2:end);
|
|
|
+ else
|
|
|
+ IndL=find(x>=valConv,Points*2-length(IndS));
|
|
|
+ end
|
|
|
+elseif isempty(IndL)|| length(IndL)<Points %if you are on the right border
|
|
|
+ if isempty(IndL) %take value from 1:10
|
|
|
+ IndS=find(x<valConv,Points*2,'last');
|
|
|
+ IndL=IndS(end);
|
|
|
+ IndS=IndS(1:end-1);
|
|
|
+ else
|
|
|
+ IndS=find(x<valConv,Points*2-length(IndL));
|
|
|
+ end
|
|
|
+end
|
|
|
+end
|