%% MAGNITUDE SPECTRUM OF SPONTANEOUS ACTIVITY (Figure 4C,E)

%1) calculate magnitude spectrum of spontaneous activity. This part is
%executed separately for each neuron.
% spontdata is a vector containing spontaneous Vm fluctuations (sampling
% frequency Fs)

transforms=[];
data=spontdata-median(spontdata(1:1000));
%a) high pass filter
cutofffreq=10; %in Hz 
filterorderfactor=20;
data_hp=genfirfilt(double(data),cutofffreq, Fs,filtorder(double(data), ...
        filterorderfactor), 'hpf', 'nuttallwin');
%b) wavelet denoise
data_wd = wden(data_hp,'heursure','s','one',5,'sym8');
%c) fourier transform
F=fft(data_wd);
L=numel(data_wd);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(data_wd,NFFT);
P2 = abs(Y/NFFT);
P1 = P2(1:NFFT/2+1);
P1(2:end-1) = 2*P1(2:end-1);
absY=P1;
f = Fs*(0:(NFFT/2))/NFFT;

absY=abs(Y(1:NFFT/2+1)).*(2/NFFT); 
transforms=[...
    transforms;...
     {f} {absY}];

%2) plot
figure;
axw=10.8;
axh=1.9201;
startx=1.5;
starty=1;
interx=1;
intery=1;
positions=[...
    startx+0*(axw+interx) starty+0*(axh+intery) axw axh;...
    startx+0*(axw+interx) starty+1*(axh+intery) axw axh;...
    startx+0*(axw+interx) starty+2*(axh+intery) axw axh;...
    startx+0*(axw+interx) starty+3*(axh+intery) axw axh;...
    startx+0*(axw+interx) starty+4*(axh+intery) axw axh;...
    startx+0*(axw+interx) starty+5*(axh+intery) axw axh;...
    startx+0*(axw+interx) starty+6*(axh+intery) axw axh;...
    startx+1*(axw+interx) starty+0*(axh+intery) axw axh;...
    startx+1*(axw+interx) starty+1*(axh+intery) axw axh;...
    startx+1*(axw+interx) starty+2*(axh+intery) axw axh;...
    startx+1*(axw+interx) starty+3*(axh+intery) axw axh;...
    startx+1*(axw+interx) starty+4*(axh+intery) axw axh;...
    startx+1*(axw+interx) starty+5*(axh+intery) axw axh;...
    startx+1*(axw+interx) starty+6*(axh+intery) axw axh;...
    ];

peakfreqs=[];
freqs05=[];
for n=1:size(transforms,1)
disp(n);
f=transforms{n,1};
absY=transforms{n,2};

%calculate median power in freq bins
binwidth = 10; 
f1=cutofffreq;
binamps=[];binfreqs=[];
while f1<max(f)-binwidth
f2=f1+binwidth;

binamps=[binamps mean(absY(f>=f1 & f<f2))];
binfreqs=[binfreqs mean([f1 f2])];

f1=f2;
end

binampsN=runav(binamps,3);
binampsN=binampsN./max(binampsN);

peakfreq=binfreqs(binampsN==max(binampsN));
peakfreqs=[peakfreqs peakfreq];

p=plot(binfreqs,binampsN);
hold on;

if(n==1)
   pP=p; 
elseif(n==6)
   pNP=p;
end

if n<6
set(p,'Color',[0 0.8 0]);
else
    set(p,'Color','k');
end

%calculate upper-cutoff frequency
[ind,t0,s0] = crossing(binampsN,binfreqs,0.5,'linear');
freq05=t0(end);
freqs05=[freqs05 freq05];

end
xlim([0 1000]);
set(gca,'Box','off','TickDir','out','Units','centimeters','Position',[1.5 1.5 7 4],'XScale','log','XTick',[20 50 100 200 500 1000],'XMinorTick','off');
xlabel('Frequency (Hz)');ylabel('Amplitude (norm.)');
l=legend([pP pNP],{'Principal','Non-principal'});

%% CALCULATE MEAN AND COEFFICIENT OF VARIATION OF INTERSPIKE INTERVAL FOR SPIKE RESPONSES (Figure 5 - Figure Supplement 1)
% spikesarray: cell array where each cell corresponds to one stimulus repetition,
% and containes a vector of spike timestamps.
binwidth = 0.1; %ms
currtimes=[0:binwidth:100];
spikematrix=ones(max(cellfun(@numel,spikesarray)),numel(spikesarray)).*NaN;
for ns=1:numel(spikesarray)
   spikestemp=(spikesarray{ns}./Fs).*1e3; 
   spikematrix(1:numel(spikestemp),ns)=spikestemp;
end
nreps = size(spikematrix,2);
minno = nreps*0.2;
isivals = ones(length(currtimes),nreps).*nan;
for j = 1:nreps
    for i = 1:(numel(find(~isnan(spikematrix(:,j)))) - 1)
            isivals(currtimes >= spikematrix(i,j) & ...
                currtimes < spikematrix(i+1,j),j) ...
                = spikematrix(i+1,j) - spikematrix(i,j);
    end
end
isimeans = ones(size(currtimes)).*NaN;
isistds = isimeans;
for i = 1:length(currtimes)
    currisi = isivals(i,~isnan(isivals(i,:)));
    if length(currisi) > minno
        isimeans(i) = mean(currisi);
        isistds(i) = std(currisi,1);
    end
end
isicovars = isistds./isimeans;

%% CALCULATE PEAKS OF SUBTHRESHOLD Vm TO MONAURAL TONES FOR ONSET AND SUSTAINED RESPONSES OF THE STIMULUS (Figure 7 - Figure supplement 1) 
% data: cell array where each cell corresponds to one sound level and 
% contains a matrix of Vm data (monaural responses) where Vm rest is subtracted 
% out (so 0 mV corresponds to Vm rest). Each column of these matrices is one 
% repetition.
% spikes: cell array with the same dimensions as data. Each cell
% corresponds to one sound level, corresponding to the cell array data, and
% contains a cell array where the number of cells is the same as the number
% of repetitions. Each of those cells contains timestamps of action potentials
% obtained during one repetition.
% spls: a vector corresponding to the number of cells in data, and
% containing the sound level values.
% Fs: sampling frequency (/s)

onset=[5 15]; onset_s=round(onset.*1e-3.*Fs);
sust=[40 50];sust_s=round(sust.*1e-3.*Fs);
sust_alt=[20 25];sust_alt_s=round(sust_alt.*1e-3.*Fs);

peakssust=[];
for i=1:numel(spls)
    currdata=data{i};
    currspikes=spikes{i};
    currspl=spls(i);
            
    for j=1:size(currdata,2)
    currstrand=currdata(:,j)-median(currdata(restdur_s(1):restdur_s(2),j));
    sps=currspikes{j};
    
    if(numel(sps)>0)
    for k=1:numel(sps)
        if(sps(k)>(onset_s(1)-spikeint_s(2)) && sps(k)<(onset_s(2)+spikeint_s(1)))
            currstrand(onset_s(1):onset_s(2))=NaN;
        end
        if(burstdur<50) %for 2 cells shorter repetitions were used
         if(sps(k)>(sust_alt_s(1)-spikeint_s(2)) && sps(k)<(sust_alt_s(2)+spikeint_s(1)))
            currstrand(sust_alt_s(1):sust_alt_s(2))=NaN;
        end   
        else
        if(sps(k)>(sust_s(1)-spikeint_s(2)) && sps(k)<(sust_s(2)+spikeint_s(1)))
            currstrand(sust_s(1):sust_s(2))=NaN;
        end
        end
    end
    end
    
    currdata(:,j)=currstrand;
    
    end
    
    currmean=nanmean(currdata,2);
    
    if(burstdur<50)
    peakssust=[peakssust; splsi(i) max(abs(currmean(onset_s(1):onset_s(2)))) max(abs(currmean(sust_alt_s(1):sust_alt_s(2))))];
    else
    peakssust=[peakssust; splsi(i) max(abs(currmean(onset_s(1):onset_s(2)))) max(abs(currmean(sust_s(1):sust_s(2))))];
    end
end

%% CALCULATE METRICS OF ILD FUNCTION (Figure 9 - Figure supplement 1)
%ildvals: vector of ild values
%ratevals: vector of spike rates corresponding to ilds
%interpolate ild function
xvals=min(ildvals):0.01:max(ildvals);
yvals = pchip(ildvals,ratevals,xvals);
% find ILD values for crossings of 50%, 10% and 90 % spike rate
thr50=min(yvals)+(0.5*(max(yvals)-min(yvals)));
yvalsshift=yvals(2:end);
ind50=find((yvals(1:end-1)>thr50 & yvalsshift<=thr50) ...
    | (yvals(1:end-1)<=thr50 & yvalsshift>thr50));
if numel(ind50)>1
   ind50=ind50(end); 
end
thr90=min(yvals)+(0.9*(max(yvals)-min(yvals)));
yvalsshift=yvals(2:end);
ind90=find((yvals(1:end-1)>thr90 & yvalsshift<=thr90) ...
    | (yvals(1:end-1)<=thr90 & yvalsshift>thr90));
if(numel(ind90)>1)
   [~,mini]=min(abs(ind90-ind50));
   ind90=ind90(mini);
end
thr10=min(yvals)+(0.1*(max(yvals)-min(yvals)));
yvalsshift=yvals(2:end);
ind10=find((yvals(1:end-1)>thr10 & yvalsshift<=thr10) ...
    | (yvals(1:end-1)<=thr10 & yvalsshift>thr10));
if(numel(ind10)>1)
    [~,mini]=min(abs(ind10-ind50));
   ind10=ind10(mini);
end
%calculate ILD function metrics
ild50=xvals(ind50);val50=yvals(ind50);
ild90=xvals(ind90);val90=yvals(ind90);
ild10=xvals(ind10);val10=yvals(ind10);
range=ild10-ild90;