%%AffVR_preproces.m is a custom written script used to detect afferent
%%neuron action potentials and ventral motor root actiity. Afferent spike
%%rates are then characterized relative to the timing of ventral motor root
%%activity

% -------------------------------------------------------------------------
% variables
% -------------------------------------------------------------------------

% analysis
    
spk_detect_mtd='MAX';    % spike detection method
spk_detect_lb=0.4;    % spike detection threshold lower bound
spk_detect_ub=1;      % spike detection threshold upper bound
spk_dur_min=0.00001;     % minimum expected spike duration (seconds)
spk_refper_min=0.0001;   % minimum expected refractory period (seconds)

vr_detect_mtd='MAX';    % ventral root detection method
vr_detect_lb=0.15;       % spike detection threshold lower bound
vr_detect_ub=1;         % spike detection threshold upper bound
vr_dur_min=0.00001;     % minimum expected spike duration (seconds)
vr_refper_min=0.0001;  % minimum expected refractory period (seconds)

bst_isi_max=0.01;            % maximum inter-spike interval within a burst (seconds)
bst_spknum_min=2;           % minimum number of spikes within a burst
bst_dur_min=0.005;            % minimum burst duration (seconds)

sb_ibi_max=0.2;            % maximum inter-burst interval within a bout (seconds)
sb_bstnum_min=3;           % minimum number of bursts within a swimming bout
sb_dur_min=0.1;            % minimum swimming bout duration (seconds)

% flags
notch_filter=1;             % remove electrical noise using notch filters
elec_stim_filter=1;              % remove electrical stimulus artifact
issue_flag=0;                   % if there is an issue while detecting afferent spikes and swimming activity

% -------------------------------------------------------------------------
% data
% -------------------------------------------------------------------------     

% select an agent file
[rcdfile,rcdpath] = uigetfile('*');

% load data
load(fullfile(rcdpath,rcdfile));

% calculate sampling rate
sr = h.dataPtsPerChan/(h.recTime(2)-h.recTime(1));

% organize data
aff = d(:,1);
vr = d(:,2);

% generate time
t = (1:length(aff))/sr;

disp(['File: ',num2str(rcdfile)]);
disp(['Data Length: ',num2str(length(aff))]);
disp(['Sampling rate: ', num2str(sr)]);


% -------------------------------------------------------------------------
% filter data
% -------------------------------------------------------------------------     

if notch_filter
    
    disp('filtering data ...');
    
    % calculate fft of raw ventral root recording
    % [fdata_vr,f_vr]=fftcal(vr,sr);
    
    % filter the data with a series of notch filters to remove electrical noise
    % and its harmonics
    [f60,f120,f180,f240,f300,f360,f420,f480,f540,f600,f660,f720,f780,f840,f900,f1020,f1140,f1260,f1380,f1500,f1620,f1740,f1860,f1980]=bpf60gen(sr);
    vr_filt = filtfilt(f60,vr); vr_filt = filtfilt(f120,vr_filt); vr_filt = filtfilt(f180,vr_filt); 
    vr_filt = filtfilt(f240,vr_filt); vr_filt = filtfilt(f300,vr_filt); vr_filt = filtfilt(f360,vr_filt);
    vr_filt = filtfilt(f420,vr_filt); vr_filt = filtfilt(f480,vr_filt); vr_filt = filtfilt(f540,vr_filt);
    vr_filt = filtfilt(f600,vr_filt); vr_filt = filtfilt(f660,vr_filt); vr_filt = filtfilt(f720,vr_filt);
    vr_filt = filtfilt(f780,vr_filt); vr_filt = filtfilt(f840,vr_filt); vr_filt = filtfilt(f900,vr_filt);
    vr_filt = filtfilt(f1020,vr_filt); vr_filt = filtfilt(f1140,vr_filt); vr_filt = filtfilt(f1260,vr_filt);
    vr_filt = filtfilt(f1380,vr_filt); vr_filt = filtfilt(f1500,vr_filt); vr_filt = filtfilt(f1620,vr_filt);
    vr_filt = filtfilt(f1740,vr_filt); vr_filt = filtfilt(f1860,vr_filt); vr_filt = filtfilt(f1980,vr_filt);
else
    
    vr_filt = vr;
    
end

vr_filt_er=vr_filt;
aff_er=aff;

% -------------------------------------------------------------------------
% detect afferent spikes
% -------------------------------------------------------------------------     

disp('detecting afferent spikes ...');

% detect afferent spikes
spktimes=spike_detect(aff_er,sr,spk_detect_mtd,spk_detect_lb,spk_detect_ub,spk_dur_min,spk_refper_min);

if size(spktimes,1)>0
    
    % generate a binary spike trace
    spkbin=zeros(length(aff_er),1);
    for k=1:size(spktimes,1)
        spkbin(spktimes(k,3))=1;
    end
    % calculate spike rate smoothing function (window size:100 ms)
    sw = 2000;
    spkrate = smooth(spkbin,sw)*sr;
    
    % calculate inter-spike interval of afferents
    isi = diff((spktimes(:,3)/sr)');
    isi = [0
        isi'];
    isi_spktimes = [isi(:,1) spktimes(:,3)];
    
else
    
    disp('no spikes were detected');
    issue_flag=1;
    
end
% -------------------------------------------------------------------------
% identify swimming bouts
% -------------------------------------------------------------------------     

disp('detecting ventral root activity ...');

% detect vr activity
vrtimes=spike_detect(abs(vr_filt_er),sr,vr_detect_mtd,vr_detect_lb,vr_detect_ub,vr_dur_min,vr_refper_min);

if size(vrtimes,1)>0

    % detect bursts in vr activity
    bst = burst_detect(vrtimes(:,3),sr,bst_isi_max,bst_spknum_min,bst_dur_min);
    % generate a binary burst activity trace
    bstbin=zeros(length(vr_filt_er),1);
    for k=1:size(bst,1)
        bstbin(bst(k,1):bst(k,2))=1;
    end
    % detect swimming bouts in vr activity
    sb = swimbout_detect(bst(:,1),bst(:,2),sr,sb_ibi_max,sb_bstnum_min,sb_dur_min);
    % generate a binary swimming bout activity trace
    sbbin=zeros(length(vr_filt_er),1);
    for k=1:size(sb,1)
        sbbin(sb(k,1):sb(k,2))=1;
    end
    
else
    
    disp('There was no swimming activity!');
    % issue_flag=1;
    
end

if issue_flag

    error('Analysis has been aborted!');

else
    
    disp('spike counting during swimming and no swimming ...');
    
    % calculate swim bout duration
    swim_dur = (sb(:,2)-sb(:,1))/sr;
    %Number of burst in a given swim
    bst_num = sb(:,3);

    % calculate swim frequency for each swim bout
    freq_sb = zeros(size(bst_num,1),1);
    for k = 1:size(bst_num,1);
        freq_sb(k) = bst_num(k)/swim_dur(k);
    end

    % calculate the duty cycle for each swim
    bst_st = bst(:,1)/sr;
    bst_end = bst(:,2)/sr;
    bst_dur=zeros(size(sb,1),1);
    for k=1:size(sb,1)
    temp = find(bst_st(:,1)>(sb(k,1)/sr) & bst_end(:,1)<(sb(k,2)/sr));
    if size(temp,1)>0
        bst_dur(k) = sum(bst_end(temp)-bst_st(temp));
    else
        bst_dur(k)=nan;
    end
    end
    
    dty_cycle=zeros(size(sb,1),1);
    for k=1:size(sb,1)
        dty_cycle(k) = bst_dur(k)/swim_dur(k);
    end
    
            
    % calculate the spike rate at times when there is a swimming bout
    spknum_insb=zeros(size(sb,1),1);        % number of spikes within a swimming bout
    spkrate_insb=zeros(size(sb,1),1);       % spike rate within a swimming bout
    for k=1:size(sb,1)
        spknum_insb(k) = sum(spktimes(:,3)>sb(k,1) & spktimes(:,3)<sb(k,2));
        spkrate_insb(k) = spknum_insb(k)*sr/(sb(k,2)-sb(k,1));
    end
    
    % calculate the isi when there is a swimming bout
    isi_insb=zeros(size(sb,1),1);
    for k=1:size(sb,1)
        temp=find(isi_spktimes(:,2)>sb(k,2));
        if size(temp,1)>0
            isi_insb(k) = isi_spktimes(temp(1));
        else
            isi_insb(k)=nan;
        end
    end
    
    % calculate the spike rate at times when there is no swimming bout
    spknum_outsb=zeros(size(sb,1)+1,1);     % number of spikes outside a swimming bout
    spkrate_outsb=zeros(size(sb,1)+1,1);    % spike rate outside a swimming bout
    spknum_outsb(1)=sum(spktimes(:,3)<sb(k,1));
    spkrate_outsb(1)=spknum_outsb(1)*sr/sb(k,1);
    for k=2:size(sb,1)
        spknum_outsb(k) = sum(spktimes(:,3)>sb(k-1,2) & spktimes(:,3)<sb(k,1));
        spkrate_outsb(k) = spknum_outsb(k)*sr/(sb(k,1)-sb(k-1,2));
    end
    spknum_outsb(end)=sum(spktimes(:,3)>sb(end,1));
    spkrate_outsb(end)=spknum_outsb(end)*sr/(length(vr_filt_er)-sb(end,2));
 
    %calculate the number of spikes & spkrate preceding a swim bout
    if size(sb,1)==1
        spknum_presb=sum(spktimes(:,3)<sb(:,1) & spktimes(:,3) >= (sb(:,1)-(sb(:,2)-sb(:,1)))); %calculates the number of spikes before a swim bout that occur within duration of the corresponding swim 
        spkrate_presb=spknum_presb*sr/(sb(:,2)-sb(:,1));
    else
    spknum_presb=zeros(size(sb,1),1);
    spkrate_presb=zeros(size(sb,1),1);
    for k=1:size(sb,1)
        spknum_presb(k) = sum(spktimes(:,3)<sb(k,1) & spktimes(:,3) >= (sb(k,1)-(sb(k,2)-sb(k,1))));
        spkrate_presb(k) = spknum_presb(k)*sr/(sb(k,2)-sb(k,1));
        end
    end
    
    %calculate the number of spikes & spkrate proceding a swim bout
    if size(sb,1)==1
        spknum_pstsb(1)=sum(spktimes(:,3)>sb(1,2) & spktimes(:,3) <= (sb(1,2)+(sb(1,2)-sb(1,1)))); %calculates the number of spikes after a swim bout that occur within duration of the corresponding swim
        spkrate_pstsb(1)=spknum_pstsb(1)*sr/(sb(1,2)-sb(1,1));
    else
        spknum_pstsb=zeros(size(sb,1),1);
        spkrate_pstsb=zeros(size(sb,1),1);
        for k=2:size(sb,1)
            spknum_pstsb(k) = sum(spktimes(:,3)>sb(k,2) & spktimes(:,3) <= (sb(k,2) + (sb(k,2)-sb(k,1))));
            spkrate_pstsb(k) = spknum_pstsb(k)*sr/(sb(k,2)-sb(k,1));        
        end
    end
   
     %isolate swim bouts where spiking occurred within a swim bouts duration before the swim
    spk_delta = [spknum_presb spknum_insb];
    
    delta_TF = zeros(size(sb,1),1);
    for k=1:size(sb,1)
    temp = find(spk_delta(k,1) > 0);  %Find indices of spknum_presb where spiking occured before a swim
        if size(temp,1)>0
            delta_TF(k) = 1;         %Index when a spike occured before the swim
        else
            delta_TF(k) = 0;         %Index when a spike did not occur before the swim
        end
    end
            
    %calculate the isi preceding a swimming bout
    isi_presb=zeros(size(sb,1),1);
    temp = find(isi_spktimes(:,2)<sb(1,1));
    if size(temp,1)>0
        isi_presb(1) = mean(isi_spktimes(temp));
    else
        isi_presb(1)=nan;
    end
    
    for k=2:size(sb,1)
    temp=find(isi_spktimes(:,2)>sb(k-1,2) & isi_spktimes(:,2)<sb(k,1));
        if size(temp,1)>0
            isi_presb(k) = mean(isi_spktimes(temp));
        else
            isi_presb(k)=nan;
        end
    end
   
    %calculate the isi when there is no swimming bout
    isi_outsb = setdiff(isi,isi_insb);

    % calcuate the time passed until the first spike after the swimming activity
    rfp_aftersb=zeros(size(sb,1),1);        % time passed until the first spike after swimming bout
    for k=1:size(sb,1)
        temp=find(spktimes(:,3)>sb(k,2));
        if size(temp,1)>0
            rfp_aftersb(k)=(spktimes(temp(1))-sb(k,2))/sr;
        else
            rfp_aftersb(k)=nan;
        end
    end
    
end

% -------------------------------------------------------------------------
% extract info from file name
% ------------------------------------------------------------------------- 
rcdfile(find(rcdfile == '.',1,'last'):end) = [];
expression = '_';
C = regexp(rcdfile,expression,'split');

f=zeros(size(sb,1),1);
dpf=zeros(size(sb,1),1);
exp=zeros(size(sb,1),1);
for k=1:size(sb,1)
    f(k) = str2num(C{1});
    dpf(k) = str2num(C{2});
    exp(k) = str2num(C{3});
end

% -------------------------------------------------------------------------
% calculate mean instantaneous spikerate using sliding window
% -------------------------------------------------------------------------     

test = [sb spknum_presb];
ind = test(:,4) > 0;
sb_spknum = test(ind,:);

for k = 2:(size(sb_spknum,1)-2);
        instspkrate1(k) = spkrate(sb_spknum(k,1)-20000);
        instspkrate2(k) = spkrate(sb_spknum(k,1)-18000);
        instspkrate3(k) = spkrate(sb_spknum(k,1)-16000);
        instspkrate4(k) = spkrate(sb_spknum(k,1)-14000);
        instspkrate5(k) = spkrate(sb_spknum(k,1)-12000);
        instspkrate6(k) = spkrate(sb_spknum(k,1)-10000);
        instspkrate7(k) = spkrate(sb_spknum(k,1)-8000);
        instspkrate8(k) = spkrate(sb_spknum(k,1)-6000);
        instspkrate9(k) = spkrate(sb_spknum(k,1)-4000);
        instspkrate10(k) = spkrate(sb_spknum(k,1)-2000);
        instspkrate11(k) = spkrate(sb_spknum(k,1));
        instspkrate12(k) = spkrate(sb_spknum(k,1)+2000);
        instspkrate13(k) = spkrate(sb_spknum(k,1)+4000);
        instspkrate14(k) = spkrate(sb_spknum(k,1)+6000);
        instspkrate15(k) = spkrate(sb_spknum(k,1)+8000);
        instspkrate16(k) = spkrate(sb_spknum(k,1)+10000);
        instspkrate17(k) = spkrate(sb_spknum(k,1)+12000);
        instspkrate18(k) = spkrate(sb_spknum(k,1)+14000);
        instspkrate19(k) = spkrate(sb_spknum(k,1)+16000);
        instspkrate20(k) = spkrate(sb_spknum(k,1)+18000);
        instspkrate21(k) = spkrate(sb_spknum(k,1)+20000);
end

m_instspkrate = [mean(instspkrate1) 0
mean(instspkrate2) 2000
mean(instspkrate3) 4000
mean(instspkrate4) 6000
mean(instspkrate5) 8000
mean(instspkrate6) 10000
mean(instspkrate7) 12000
mean(instspkrate8) 14000
mean(instspkrate9) 16000
mean(instspkrate10) 18000
mean(instspkrate11) 20000
mean(instspkrate12) 22000
mean(instspkrate13) 24000
mean(instspkrate14) 26000
mean(instspkrate15) 28000
mean(instspkrate16) 30000
mean(instspkrate17) 32000
mean(instspkrate18) 34000
mean(instspkrate19) 36000
mean(instspkrate20) 38000
mean(instspkrate21) 40000];

sem_instspkrate = [std(instspkrate1/sqrt(length(instspkrate1)))
    std(instspkrate2/sqrt(length(instspkrate2)))
    std(instspkrate3/sqrt(length(instspkrate3)))
    std(instspkrate4/sqrt(length(instspkrate4)))
    std(instspkrate5/sqrt(length(instspkrate5)))
    std(instspkrate6/sqrt(length(instspkrate6)))
    std(instspkrate7/sqrt(length(instspkrate7)))
    std(instspkrate8/sqrt(length(instspkrate8)))
    std(instspkrate9/sqrt(length(instspkrate9)))
    std(instspkrate10/sqrt(length(instspkrate10)))
    std(instspkrate11/sqrt(length(instspkrate11)))
    std(instspkrate12/sqrt(length(instspkrate12)))
    std(instspkrate13/sqrt(length(instspkrate13)))
    std(instspkrate14/sqrt(length(instspkrate14)))
    std(instspkrate15/sqrt(length(instspkrate15)))
    std(instspkrate16/sqrt(length(instspkrate16)))
    std(instspkrate17/sqrt(length(instspkrate17)))
    std(instspkrate18/sqrt(length(instspkrate18)))
    std(instspkrate19/sqrt(length(instspkrate19)))
    std(instspkrate20/sqrt(length(instspkrate20)))
    std(instspkrate21/sqrt(length(instspkrate21)))];

% -------------------------------------------------------------------------
% figures
% -------------------------------------------------------------------------     

disp('generating figures ...');

figure;
hold on;
plot(t,aff_er,'k');
plot(t(spktimes(:,3)),aff_er(spktimes(:,3)),'ro','MarkerFaceColor','r');
hold off;
xlabel('time (s)'); title('Spike detection');

figure;
hold on;
plot(t,vr_filt_er,'k');
plot(t(vrtimes(:,3)),vr_filt_er(vrtimes(:,3)),'ro','MarkerFaceColor','r');
plot(t,max(vr_filt_er)*sbbin,'g','Linewidth',2);
hold off;
xlabel('time (s)'); title('Swimming bout detection');

figure;
hold on;
plot(t,spkrate,'r','Linewidth',2);
plot(t,20*sbbin,'b','Linewidth',2);
legend('swimming bout','spike rate'); 
hold off;
xlabel('time (s)'); title('Swimming activity versus spike rate');

figure
hold on
for k = 2:(size(sb_spknum,1)-2);
    plot(spkrate((sb_spknum(k,1)-20000):(sb_spknum(k,1)+20000)),'LineWidth',1);
end
xlabel('time (sec)'); ylabel('instantaneous spike rate (hz)'); title('Spike rate dynamics during each swim');
xticks([0 10000 20000 30000 40000])
xticklabels([-1 -0.5 0 0.5 1])
xlim([0 40000])
ylim([0 25]);
hold off

figure
hold on
plot(m_instspkrate(:,2),m_instspkrate(:,1),'LineWidth',3)
errorbar(m_instspkrate(:,2),m_instspkrate(:,1),sem_instspkrate(:,1),'LineWidth',2)
xlabel('time (sec)'); ylabel('instantaneous spike rate (hz)'); title('Mean spike rate dynamics during swimming');
xticks([0 10000 20000 30000 40000])
xticklabels([-1 -0.5 0 0.5 1])
ylim([0 20]);
hold off

% -------------------------------------------------------------------------
% save results
% -------------------------------------------------------------------------   

% info
res.info.rcdfile=rcdfile;
res.info.rcdpath=rcdpath;
res.info.sr=sr;
res.info.fish=f(1,1);
res.info.age=dpf(1,1);

% variables to analyze data   
res.var.spk_detect_mtd=spk_detect_mtd;    % spike detection method
res.var.spk_detect_lb=spk_detect_lb;    % spike detection threshold lower bound
res.var.spk_detect_ub=spk_detect_ub;      % spike detection threshold upper bound
res.var.spk_dur_min=spk_dur_min;     % minimum expected spike duration (seconds)
res.var.spk_refper_min=spk_refper_min;   % minimum expected refractory period (seconds) 

res.var.vr_detect_mtd=vr_detect_mtd;    % ventral root detection method
res.var.vr_detect_lb=vr_detect_lb;       % spike detection threshold lower bound
res.var.vr_detect_ub=vr_detect_ub;         % spike detection threshold upper bound
res.var.vr_dur_min=vr_dur_min;     % minimum expected spike duration (seconds)
res.var.vr_refper_min=vr_refper_min;  % minimum expected refractory period (seconds)

res.var.bst_isi_max=bst_isi_max;            % maximum inter-spike interval within a burst (seconds)
res.var.bst_spknum_min=bst_spknum_min;           % minimum number of spikes within a burst
res.var.bst_dur_min=bst_dur_min;            % minimummean burst duration (seconds)

% flags
%res.flags.notch_filter=notch_filter;
%res.issue_flag=issue_flag;

% preliminary analysis
res.preanal.spktimes=spktimes;
res.preanal.isi=isi;
res.preanal.vrtimes=vrtimes;
res.preanal.bst=bst;
res.preanal.sb=sb;
res.preanal.sb_freq=freq_sb;
res.preanal.dty_cycle=dty_cycle;
res.preanal.spknum_delta=spk_delta;
res.preanal.subset_delta=delta_TF;
res.preanal.spknum_insb=spknum_insb;
res.preanal.spkrate_insb=spkrate_insb;
res.preanal.isi_presb=isi_presb;
res.preanal.isi_insb=isi_insb;
res.preanal.spknum_outsb=spknum_outsb;
res.preanal.spkrate_outsb=spkrate_outsb;
res.preanal.spkrate_presb=spkrate_presb;
res.preanal.spkrate_pstsb=spkrate_pstsb;
res.preanal.isi_outsb=isi_outsb;
res.preanal.rfp_aftersb=rfp_aftersb;


prompt = 'Do you want to save the results ? Y/N [Y]: ';
str = input(prompt,'s');
if isempty(str)
    str = 'Y';
end

if str=='Y'
    save(['./preprocess_output/',rcdfile,'.mat'],'res');
end

% -------------------------------------------------------------------------
% save to workspace
% -------------------------------------------------------------------------  

fish = [f];
    %fish];
    
age = [dpf];
   % age];

burst = [bst];  %Start_End_#Spikes
   % burst];

swimbout_st = [(sb(:,1)/sr)];  %Start_End_#Bursts
    % swimbout_st];
swimbout_end = [(sb(:,2)/sr)];
    % swimbout_end];
sb_dur = [swimbout_end-swimbout_st];
swimbout_bst = [sb(:,3)];
   % swimbout_bst];

sb_dty = [dty_cycle];
   % sb_dty];
sb_freq = [freq_sb];
   % sb_freq];

spknum_delta = [spk_delta]; %Pre-Swim Spike#_During Swim Spike#
   % spknum_delta];
    
affspkrate_insb = [spkrate_insb];
   % affspkrate_insb];
affspkrate_outsb = [spkrate_outsb];
   % affspkrate_outsb];
affspkrate_presb = [spkrate_presb];
   % affspkrate_presb];
affspkrate_pstsb = [spkrate_pstsb];
   % affspkrate_pstsb];
    
%total_isi = [isi];
    %total_isi];
%affisi_presb = [isi_presb];
    %affisi_presb];
%affisi_insb = [isi_insb];
    %affisi_insb];
%affisi_outsb = [isi_outsb];
    %affisi_outsb];

rfp = [rfp_aftersb];
    %rfp];
    
%Swimbout Start Time(sec)-Swimbout End Time(sec)-#Burst w/in Swimbout-Swimbout Frequency(Hz)-Aff Freq Swim(Hz)-Aff Freq No Swim(Hz)-Aff Interspike Interval Swim(sec)-Aff Interspike Interval No Swim(sec)-Refractory Period(sec)-#Spikes PreSwim-#Spikes Swim-Subset Indicator (1,0)
temp = [exp, fish, age, sb_dur, sb_freq, sb_dty, affspkrate_presb, affspkrate_insb, affspkrate_pstsb,spknum_delta(:,1), spknum_delta(:,2)];

ind0 = temp(:,7) > 0;
temp = temp(ind0,:);

rel_inhib = zeros(size(temp,1),1);
for k = 1:size(temp,1)
    rel_inhib(k) = 1-(temp(k,8)/temp(k,7));
        if rel_inhib(k) < 0
            rel_inhib(k) = 0;
        end
end

data_out = [temp rel_inhib zeros(size(rel_inhib,1),1)];
header = {'Date','Fish','Age','SwimDuration','SwimFrequency','SwimDutyCycle','PreSwimSpikeRate','SwimSpikeRate','PostSwimSpikeRate','PreSwimSpikeNumber','SwimSpikeNumber','RelativeInhibition','Ablation'};
data_out = [header;num2cell(data_out)];

xlswrite('data_out.xls',data_out);

%clearvars -except fish age sb_dur sb_freq sb_dty affspkrate_presb affspkrate_insb affspkrate_pstsb spknum_delta(:,1) spknum_delta(:,2) rfp rel_inhib data_out
      
