function [NL_output]=BC_NL(STA_output)
%[NL_output]=BC_NL(STA_output)
%this function computes the output nonlinearity (NL) from the linear-nonlinear model (LN-model)
%for continous voltage signals of bipolar cells (BC).
%Inputs:
%               STA_output = output of the BC_STA function; containing the
%               STA and stimulus signal
%
%Outpus:        NL_output = structure containing the output nonlinearity
%
%
% Note: the function is built only for cells recorded without frozen frames.
% But additional comments are made for cells with frozen frames.
% code by HS, last modified March 2021

%-------------------------------------------------------------------------------
%% 1. convolve zoom STA with stimulus (to get the generator signal)
%get the linear filter (normalized, see Methods of Schreyer&Gollisch,2021)
linearFiler=normSTA(STA_output.STA2D_zoom);
%get stimulus signal
stimulus=STA_output.stimulus2D_zoom; 
%%if you have frozen frames, only the stimulus signal of the running frames
%%is used for the convolution. Also the stimulus signal has to be convolved
%%for each trial of the running noise separately (because it is not a
%%continous signal -> it was interupted by the frozen noise).

%do the convolution with the stimulus signal
generator_signal=filter2(linearFiler,stimulus,'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,stimulus,'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. get the y axis of the NL
%get the voltage responses
y_axisNL_unsorted=STA_output.membranPotAVG(STA_output.STA_fnbr:end); %sthe first value is at the length of the filter

%-------------------------------------------------------------------------------
%% 3. sort the axis and do the binning
%sort the x axis
[NL_xAxis,indx]=sort(generator_signal);
%sort the corresponding y axis
NL_yAxis=y_axisNL_unsorted(indx);
%do a binning with percentile to have same amount of data in each bin
%make 40 bins
nbrbins=40;
bins= doBin(NL_xAxis,NL_yAxis,nbrbins);

%-------------------------------------------------------------------------------
%% 4. plot the nonlinearity
figure;
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: filename: ' STA_output.filename])
axis tight
xlabel('generator signal')
ylabel('voltage (mV)')

%-------------------------------------------------------------------------------
%% 5. add variables into the output
NL_output.NL_xAxis=NL_xAxis;
% STA_output.STA=STA; %not  normalized
NL_output.NL_yAxis=NL_yAxis;
NL_output.bins=bins; %for reshapeing to 3D structure if needed

end