123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181 |
- function [naturalMovie_output]=BC_predNaturalMovies(loadData,filename,movieSignal,STA_output,NL_output,Hz)
- %[naturalMovie_output]=BC_predNaturalMovies(loadData,filename,movieSignal,STA_output,NL_output,Hz)
- %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
- % movieSignal = the preprocessed natural movie signal the
- % output of the BC_getNaturalMovies.m. Important that the
- % movie signal has the same dimensions as the STA_output
- % 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
- % Hz = movie update rate given in excel file
- %
- %Outpus: naturalMovie_output = structure containing the output of the
- % prediction
- %
- % code by HS last modiefied March 2021
- %% Constants
- findBreak=500; %the movie is repeated multiple times. Between the movies there is
- %a break/gray screen of >0.5s. To find the start of each movie trial, you find the gray screen between the movies
- %-------------------------------------------------------------------------------
- %% 1. cut the natural movie signal to the size of the filter
- %in the BC_STA.m function, the STA was cut. Here the natural movie signal
- %has to be cut to the same dimensions.
- movieSignal_zoom=movieSignal(STA_output.rangeY_zoom,STA_output.rangeX_zoom,:);
- %make the signal 2D for prediction
- movieSignal2D_zoom=reshape(movieSignal_zoom,size(movieSignal_zoom,1)*...
- size(movieSignal_zoom,2),size(movieSignal_zoom,3));
- %-------------------------------------------------------------------------------
- %% 2. compute the prediction
- %% 2.1 convolution
- %normalize the filter
- linearFiler=normSTA(STA_output.STA2D_zoom);
- %convolve the natural movie stimulus sequence with the
- %normalized STA.
- stim2pred=filter2(linearFiler,movieSignal2D_zoom,'valid');
- %%IMPORTANT: if you have a best spatial and temporal component as described in Point 7
- %in the function BC_STA.m (see Method of Manuscript), this part hast to be done in two lines, e.g.:
- %NL_PredXaxis2=filter2(sp_bestnorm,movieSignal2D_zoom,'valid'); %first
- %convolve the spatial dimension
- %generator_signal=filter2(temp_bestnorm,NL_PredXaxis2,'valid'); %then the
- %temporal dimension
- %%Also note, to avoid noise contributions from pixels outside the receptive
- %%field,it is important to cut the STA into a smaller region around the maximum pixel, this
- %%gives a much cleaner NL-> see BC_STA.m and Methods in Schreyer&Gollisch (2021).
- %% 2.2 compute the prediction through linear interpolation
- NL_x=NL_output.bins.NL_xbin;
- NL_y=NL_output.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
- %-------------------------------------------------------------------------------
- %% 4. compute the average response to the move
- %% 4.1 load movie response
- data=load([loadData,filename]); %load the file
- %sort the different pulses for the different trials
- diffttls=diff(data.ttls_ms);
- nbrTrials=find(diffttls>findBreak); %to find out which pulse correspond to a new trial (gray break)
- nbrFrames=diff(nbrTrials);
- if isempty(nbrFrames)
- nbrFrames=nbrTrials;
- end
- %I get the membrane potential per trial & frames per trial
- reOrderFrames=nan(nbrFrames(1),length(nbrTrials));
- membranePerTrial_origRes=cell(1,length(nbrTrials));
- for kk=1:length(nbrTrials)
- if kk==1 %first trial
- reOrderFrames(:,kk)=data.ttls_ms(1:nbrTrials(kk));
- membranePerTrial_origRes{kk}=data.membPot(reOrderFrames(1,kk):reOrderFrames(end,kk)-1); %-1 because the last frame is the start of the gray screen
- else
- reOrderFrames(:,kk)=data.ttls_ms(nbrTrials(kk-1)+1:nbrTrials(kk)); %end of previous trial +1 (starts at 999)
- membranePerTrial_origRes{kk}=data.membPot(reOrderFrames(1,kk):reOrderFrames(end,kk)-1);
- end
- end
- %I only include trials here that here fully finished until the break!
- %------------------------------------------------------------------------------------------
- %% 4.2 Natural Movie Build average per frame
- %make one value for the membrane potential between two stimulus signal (two
- %frames), where the same image was shown
- memPperTrial_Movie=avgPerFrameNatMovies(data.membPot,...
- reOrderFrames,nbrFrames(1),Hz);
- memP_AVGMovie=mean(memPperTrial_Movie,2);
- %------------------------------------------------------------------------------------------
- %% 5 compute the pearson correlation coefficient
- %pearson correlation
- corrcoeff_temp=corrcoef(memP_AVGMovie(STA_output.STA_fnbr:end),predMemP_y);
- pearCorrSqrt=corrcoeff_temp(2)^2;
- %not with this method the prediction value is a bit lower than what we report in the
- %manuscript. To get the same values as in the mansucript, you have to
- %cut 3sigma around the RF and perform the SVD (see also Manuscript and
- %BC_STA, BC_NL)
- %-------------------------------------------------------------------------------
- %% 6. plot the output
- figure;
- %plot the first frame of the move
- subplot(2,1,1)
- imagesc(movieSignal(:,:,1))
- colormap('gray'); %how the sitmulus was shown
- axis equal; axis tight;
- xticks([]); yticks([]);
- %plot the prediction trace
- subplot(2,1,2)
- %only for visual comparison we normalize the predicted trace to have the
- %same mean and std as the response to the movie.
- x=predMemP_y;
- m2=mean(memP_AVGMovie);
- s2=std(memP_AVGMovie);
- s1=std(x);
- m1=mean(x);
- corrPredMemPred=m2+(x-m1)*(s2/s1);
- plot(memP_AVGMovie(STA_output.STA_fnbr:end),'k') %original membrane potential
- hold on;
- plot(corrPredMemPred,'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
- naturalMovie_output.pearCorrSqrt=pearCorrSqrt;
- naturalMovie_output.memP_AVGMovie=memP_AVGMovie;
- 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
- % %make SVD of the filter
- % [sp,sigma,temp] = svd(STA_output.STA2D_zoom,'econ');
- % sp=sp(:,1);
- % imagesc(reshape(sp,size(STA_output.rangeY_zoom,2),size(STA_output.rangeX_zoom,2)));
- % temp=temp(:,1);
- % %linear filters
- % linearFiler1=normSTA(sp);
- % linearFiler2=normSTA(temp);
- % stim2pred=filter2(linearFiler1,movieSignal2D_zoom,'valid');
- % stim2pred=filter2(linearFiler2',stim2pred,'valid');
|