% DC-ERG Analysis program 
% JM 4-28-2020
%
% LOAD DATA
% PERFORM LINEAR DRIFT SUBTRACTION
% FIND DC-ERG COMPONENTS
% SAVES VALUES TO STRUCTURE: Parameters
% EXPORTS DATA TO EXCEL  
%
% Note:  This program was tested with 10cd/m2 (1 log cd/m2) light stimulation.
% Therefore the off response is defined as a minimum. 
% See lines 344 and 359. This light level
% is submaximal and well below the threshold where the 
% off response reverses polarity (> 2.4 log cd/m2 for WT).
%___________________________________________________________

%%%%%%%%%%%%%%%%%%%%%% LOAD DATA %%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all
clc

endshift=100;
ends=0;
n=1;

reply = 'y';
[c,pathc] = uigetfile({'*.csv'},'Select the data','V:\data','MultiSelect', 'on');

if reply == 'N'||'n';
    c=cellstr(c);
elseif reply == 'Y'||'y';
    display('. . . . calculating');
end

clear i;

%create Excel filename for export - append current month, day, year, time
Filename = sprintf('DCERG_Results_%s.xlsx', datestr(now,'mm-dd-yyyy HH-MM'));
    
for i=1:length(c)
    file=[pathc c{i}];
    fid = fopen(file);
    [pathstr, name, ext] = fileparts(file);
    
    % Load a data file even if it has comments.
    
    format long
    
    %# find and open file
    
    fseek(fid, 0, -1);
    
    Header1 = textscan(fid, '%s %s %s %s %s', 1);
    
    Header2 = textscan(fid, '%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s', 1);
    
    Header3 = textscan(fid, '%s %s %s %s %s %s %s %s %s %s %s %s', 1);
    
    Skip1 = textscan(fid, '%s %s %s %s %s %s %s %s %s %s %s %s', 1);
    
    Skip2 = textscan(fid, '%s %s %s %s %s %s %s %s %s %s %s %s', 1);
    
    Skip3 = textscan(fid, '%s %s %s %s %s %s %s %s %s %s %s %s', 1);
    
    Skip4 = textscan(fid, '%s %s %s %s %s %s %s %s %s %s %s %s', 1);
    
    %# find empty cells
    
    clear i
    
    for i=1:length(Header3)
        emptyCells(i) = cellfun(@isempty,Header3{i});
    end
    
    %# remove empty cells
    
    Header3(emptyCells) = [];
    
    if length(Header3)==10
        
        temp=textscan(fid, '%n %n %n %n %n' ,-1); % NEW FILES HAVE 3 columns
        
    elseif length(Header3)==6
        
        temp=textscan(fid, '%n %n %n' ,-1); % NEW FILES HAVE 3 columns
        
    elseif length(Header3)==12
        
        temp=textscan(fid, '%n %n %n %n %n %n' ,-1); % OLD FILES HAVE 6 columns
        
        %# close the file again
        
        fclose(fid);
        
    end
    
    clear data
    
    % catenate, because textscan returns a column of cells for each column in the data
    
    if length(Header3)==10
        
        data = cat(2,temp{1:5});  % NEW DATA 3 columns
        
        data(:,5)=[];
        
        data(:,2)=[];
        
    elseif length(Header3)==6
        
        data = cat(2,temp{1:3});  % NEW DATA 3 columns
        
    elseif length(Header3)==12
        
        data = cat(2,temp{4:6}); % OLD DATA 6 columns
        
    end
    
    data=data(1:end,:);
    
    %%%%%%%%%%%%%%%%%%%%%%%% PLOT RAW DATA  %%%%%%%%%%%%%%%%%%%%%%%%%%
    
    time=data(:,1)./60000;  % time is 540,000 ms (convert it to 9 min)
    
    data(:,1)=data(:,1)./60000; % convert time to min
    data(:,2)=data(:,2)./10^6; % convert nV to mV
    data(:,3)=data(:,3)./10^6; % convert nV to mV
    
    figure
    subplot(2,1,1)
    Lrawdata=line(data(:,1),data(:,2));
    set(Lrawdata,'Color',[0 0 0]);
    set(Lrawdata,'LineWidth',1.5);
    
    subplot(2,1,2)
    Rrawdata=line(data(:,1),data(:,3));
    set(Rrawdata,'Color',[1 0 0]);
    set(Rrawdata,'LineWidth',1.5);
    
    set(0,'DefaultLegendAutoUpdate','off')
    
    %str2='.eps'
    %nameRaw = strcat(str1, str2)
    %print(nameRaw, '-depsc2')
    
    %%%%%%%%%%%%%%%%%%%%%%%% PLOT LINEAR FIT  %%%%%%%%%%%%%%%%%%%%%%%%%%
    
    clear x1
    clear y1
    clear p1
    clear f1
    
    x1=data(10:260,1);
    
    y1=data(10:260,2);
    
    subplot(2,1,1)
    hold on
    plot(x1,y1,'g.');
    
    [c2,gof2,output2] = fit(x1,y1,'poly1');
    
    hold on
    
    plot(c2,'m');
    
    coeff1 = coeffvalues(c2)';
    
    confinterval1 = confint(c2)';
    
    %%%%%%%%%%%%%%%%%%%  CALCULATE LINEAR FIT %%%%%%%%%%%%%%%
    b1=coeff1(1);
    
    b2=coeff1(2);
    
    yhat = @(coeff1,x)(b1*x+b2);
    
    yfit=yhat(coeff1,time);
    
    hold on
    
    % show results
    
    plot(time,yfit,'c');
    
    hXlabel=xlabel('Time (min)');
    hYlabel=ylabel('Amplitude (mV)')
    hTitle= title(num2str(name),'Interpreter','none');
    
    set( gca                       , ...
        'FontName'   , 'Arial' );
    set([hTitle, hXlabel, hYlabel], ...
        'FontName'   , 'Arial');
    set([hXlabel, hYlabel]  , ...
        'FontSize'   , 15          );
    set( hTitle                    , ...
        'FontSize'   , 18          , ...
        'FontWeight' , 'bold'      );
    
    set(gca, ...
        'Box'         , 'off'     , ...
        'TickDir'     , 'out'     , ...
        'TickLength'  , [.02 .02] , ...
        'XMinorTick'  , 'off'      , ...
        'YMinorTick'  , 'off'      , ...
        'YGrid'       , 'off'      , ...
        'XColor'      , [.3 .3 .3], ...
        'YColor'      , [.3 .3 .3], ...
        'LineWidth'   , 1.0        );
    
    hLegend = legend( ...
        [Lrawdata], ...
        'Left eye (Ch.1)' , ...
        'location', 'bestoutside' );
    
    set([hLegend, gca]             , ...
        'FontSize'   , 8           );
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    clear x2
    clear y2
    clear p2
    clear f2
    
    x2=data(10:260,1);
    
    y2=data(10:260,3);
    
    subplot(2,1,2)
    hold on
    plot(x2,y2,'g.')
    
    [c3,gof3,output3] = fit(x2,y2,'poly1')
    
    hold on
    
    plot(c3,'m')
    
    coeff2 = coeffvalues(c3)';
    
    confinterval2 = confint(c3)';
    
    b01=coeff2(1);
    
    b02=coeff2(2);
    
    yhat2 = @(coeff2,x)(b01*x+b02);
    
    yfit2=yhat2(coeff2,time);
    
    
    hold on
    
    % show results
    
    plot(time,yfit2,'c');
     
    hXlabel=xlabel('Time (min)');
    hYlabel=ylabel('Amplitude (mV)')
    hTitle= title(num2str(name),'Interpreter','none');
    
    set( gca                       , ...
        'FontName'   , 'Arial' );
    set([hTitle, hXlabel, hYlabel], ...
        'FontName'   , 'Arial');
    set([hXlabel, hYlabel]  , ...
        'FontSize'   , 15          );
    set( hTitle                    , ...
        'FontSize'   , 18          , ...
        'FontWeight' , 'bold'      );
    
    set(gca, ...
        'Box'         , 'off'     , ...
        'TickDir'     , 'out'     , ...
        'TickLength'  , [.02 .02] , ...
        'XMinorTick'  , 'off'      , ...
        'YMinorTick'  , 'off'      , ...
        'YGrid'       , 'off'      , ...
        'XColor'      , [.3 .3 .3], ...
        'YColor'      , [.3 .3 .3], ...
        'LineWidth'   , 1.0         );
    
    hLegend = legend( ...
        [Rrawdata], ...
        'Right eye (Ch.2)' , ...
        'location', 'bestoutside' );
    
    set([hLegend, gca]             , ...
        'FontSize'   , 8           );
    
    set(gcf, 'PaperPositionMode', 'auto');
    
    %testfit=['FLinearfit ' name];
    %saveas(gcf,testfit,'tiffn');
    
    %%%%%%%%%%%%%%%%%%%%%%    APPLY DRIFT SUBRACTION   %%%%%%%%%%%%%%%
    
    %drift subtract and 50-point moving average (50 pts =  5.292 s)
    Eye1=smooth(data(:,2)-yfit,50); 
    Eye2=smooth(data(:,3)-yfit2,50);
    
    %plot left eye (ch1)
    figure
    subplot(2,1,1)
    LEye1=line(data(:,1),Eye1);
    set(LEye1,'Color',[0 0 0]);
    set(LEye1,'LineWidth',1.5);
    
    %plot right eye (ch2)
    subplot(2,1,2)
    LEye2=line(data(:,1),Eye2);
    set(LEye2,'Color',[1 0 0]);
    set(LEye2,'LineWidth',1.5);
    
    
    %%%%%%%%%%%%%%%%%%%%%%    FIND DC-ERG PARAMETERS   %%%%%%%%%%%%%%%
    
    [cwave1,cindex1]=max(Eye1(300:600));
    
    [FOmin1,FOindex1]=min(Eye1(600:1200));
    
    [cwave2,cindex2]=max(Eye2(300:600));
    
    [FOmin2,FOindex2]=min(Eye2(600:1200));
    
    %plot left eye (ch1)
    subplot(2,1,1)
    hold on
    plot(data(cindex1+300,1),cwave1,'co');
    plot(data(FOindex1+600,1),FOmin1,'ms');
    
    %plot right eye (ch2)
    subplot(2,1,2)
    hold on
    plot(data(cindex2+300,1),cwave2,'co');
    plot(data(FOindex2+600,1),FOmin2,'ms');
    
    FO1=cwave1-FOmin1;
    FO2=cwave2-FOmin2;
    
    [LPmax1,LPindex1]=max(Eye1(4480:4600));
    
    [Offmin1,Offindex1]=min(Eye1(4480:4900));
    
    subplot(2,1,1)
    
    hold on
    
    plot(data(LPindex1+4480,1),LPmax1,'bd');
    
    hold on
    
    plot(data(Offindex1+4480,1),Offmin1,'gp');
    
    [LPmax2,LPindex2]=max(Eye2(4480:4600));
    
    [Offmin2,Offindex2]=min(Eye2(4480:4900));
    
    subplot(2,1,2)
    
    hold on
    
    plot(data(LPindex2+4480,1),LPmax2,'bd');
    
    hold on
    
    plot(data(Offindex2+4480,1),Offmin2,'gp');
    
    if FOmin1>0
        LP1=LPmax1-FOmin1;
    else
        LP1=LPmax1+abs(FOmin1);
    end
    
    if FOmin2>0
        LP2=LPmax2-FOmin2;
    else
        LP2=LPmax2+abs(FOmin2);
    end
    
    Off1=abs(LPmax1-Offmin1);
    Off2=abs(LPmax2-Offmin2);
    
    %round the amplitude values to the nearest hundredth
    cwave1=round(cwave1,2);
    cwave2=round(cwave2,2);
    FO1=round(FO1,2);
    FO2=round(FO2,2);
    LP1=round(LP1,2);
    LP2=round(LP2,2);
    Off1=round(Off1,2);
    Off2=round(Off2,2);
    
    %create structure Parameters to store all data values
    Parameters.Name=name;
    Parameters.Data=data;
    Parameters.SmoothData(:,1)=time;
    Parameters.SmoothData(:,2)=Eye1;
    Parameters.SmoothData(:,3)=Eye2;
    Parameters.Cwave(1,1)=round(data(cindex1+300,1),3)-0.5; %Subtract 0.5 min (30s) prestimulus
    Parameters.Cwave(1,2)=cwave1;
    Parameters.Cwave(2,1)=round(data(cindex2+300,1),3)-0.5; %Subtract 0.5 min (30s) prestimulus
    Parameters.Cwave(2,2)=cwave2;
    Parameters.FO(1,1)=round(data(FOindex1+600,1),3)-0.5; %Subtract 0.5 min (30s) prestimulus
    Parameters.FO(1,2)=FO1;
    Parameters.FO(2,1)=round(data(FOindex2+600,1),3)-0.5; %Subtract 0.5 min (30s) prestimulus 
    Parameters.FO(2,2)=FO2;
    Parameters.LP(1,1)=round(data(LPindex1+4480,1),3)-0.5; %Subtract 0.5 min (30s) prestimulus
    Parameters.LP(1,2)=LP1;
    Parameters.LP(2,1)=round(data(LPindex2+4480,1),3)-0.5; %Subtract 0.5 min (30s) prestimulus
    Parameters.LP(2,2)=LP2;
    Parameters.Off(1,1)=round(data(Offindex1+4480,1),3)-0.5; %Subtract 0.5 min (30s) prestimulus
    Parameters.Off(1,2)=Off1;
    Parameters.Off(2,1)=round(data(Offindex2+4480,1),3)-0.5; %Subtract 0.5 min (30s) prestimulus
    Parameters.Off(2,2)=Off2;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %uncomment these two lines if you wish to store each Parameters file
    %after every iteration of the loop (for every mouse)
    
    %%test2=['Linearfit ' name];
    %%save(test2, 'Parameters');
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    axis auto
    ax1=subplot(2,1,1)
    h2legend=legend('Left Eye (Ch.1)',['cwave1:  ' num2str(cwave1)],['FO1:  ' num2str(FO1)],['LP1:  ' num2str(LP1)],['Off1:  ' num2str(Off1)], 'Location','Bestoutside')
    h2Xlabel=xlabel('Time (min)');
    h2Ylabel=ylabel('Amplitude (mV)')
    h2Title= title(num2str(name),'Interpreter','none');
    
    set( ax1                       , ...
        'FontName'   , 'Arial' );
    set([h2Title, h2Xlabel, h2Ylabel], ...
        'FontName'   , 'Arial');
    set([h2Xlabel, h2Ylabel]  , ...
        'FontSize'   , 15          );
    set( h2Title                    , ...
        'FontSize'   , 18          , ...
        'FontWeight' , 'bold'      );
    
    set(ax1, ...
        'Box'         , 'off'     , ...
        'TickDir'     , 'out'     , ...
        'TickLength'  , [.02 .02] , ...
        'XMinorTick'  , 'off'      , ...
        'YMinorTick'  , 'off'      , ...
        'YGrid'       , 'off'      , ...
        'XColor'      , [.3 .3 .3], ...
        'YColor'      , [.3 .3 .3], ...
        'LineWidth'   , 1.0        );
    
    set([h2legend, gca]             , ...
        'FontSize'   , 8           );
    
    axis auto
    ax2=subplot(2,1,2)
    
    h3legend=legend('Right Eye (Ch.2)',['cwave2:  ' num2str(cwave2)],['FO2:  ' num2str(FO2)], ['LP2:  ' num2str(LP2)],['Off2:  ' num2str(Off2)],'Location','BestOutside')
    
    h3Xlabel=xlabel('Time (min)');
    h3Ylabel=ylabel('Amplitude (mV)')
    h3Title= title(num2str(name),'Interpreter','none');
    
    set( ax2                       , ...
        'FontName'   , 'Arial' );
    set([h3Title, h3Xlabel, h3Ylabel], ...
        'FontName'   , 'Arial');
    set([h3Xlabel, h3Ylabel]  , ...
        'FontSize'   , 15          );
    set( h3Title                    , ...
        'FontSize'   , 18          , ...
        'FontWeight' , 'bold'      );
    
    set(ax2, ...
        'Box'         , 'off'     , ...
        'TickDir'     , 'out'     , ...
        'TickLength'  , [.02 .02] , ...
        'XMinorTick'  , 'off'      , ...
        'YMinorTick'  , 'off'      , ...
        'YGrid'       , 'off'      , ...
        'XColor'      , [.3 .3 .3], ...
        'YColor'      , [.3 .3 .3], ...
        'LineWidth'   , 1.0        );
    
    set([h3legend, gca]             , ...
        'FontSize'   , 8           );
        
    
    %%%%%% WRITE TO EXCEL %%%%%%%%%%%%
    format short
    sheet=n;
    
    %Amplitude (mV)
    cwaveVal=[cwave1;cwave2];
    foVal=[FO1;FO2];
    lpVal=[LP1;LP2];
    offVal=[Off1;Off2];
            
    %Time (min)
    cwavetme=[round(data(cindex1+300,1),3)-0.5; round(data(cindex2+300,1),3)-0.5];
    fotme=[round(data(FOindex1+600,1),3)-0.5; round(data(FOindex2+600,1),3)-0.5];
    lptme=[round(data(LPindex1+4480,1),3)-0.5; round(data(LPindex2+4480,1),3)-0.5];
    offtme=[round(data(Offindex1+4480,1),3)-0.5; round(data(Offindex2+4480,1),3)-0.5];
    
    %create a temp variable to be overwritten with the Mouse ID (name)
    blank=[n];
        
    %create column names for table
    varNames={'Cwave_mV','FO_mV','LP_mV','OFF_mV'}; 
    varNames2={'Cwave_min','FO_min','LP_min','OFF_min'};
    varNames3={name};
    varNames6={'Time_min','Eye1_Left','Eye2_Right'};
    
    %define table values
    T= table(cwaveVal,foVal,lpVal,offVal,'VariableNames',varNames);
    T2= table(cwavetme,fotme,lptme,offtme, 'VariableNames',varNames2);
    T3= table(blank,'VariableNames',varNames3);
    T6= table(time,Eye1,Eye2, 'VariableNames',varNames6);
    
    %write table to excel
    writetable(T,Filename,'Sheet',name,'Range','A3:H5');
    writetable(T2,Filename,'Sheet',name,'Range','A7:D9');
    writetable(T3,Filename,'Sheet',name,'Range','A1:B1');
    writetable(T6,Filename,'Sheet',name,'Range','F1:H5406');
   
    %Store amplitude values for summary data
    cstor(:,n)=cwaveVal;
    fstor(:,n)=foVal;
    lstor(:,n)=lpVal;
    ostor(:,n)=offVal;
    
    %Store time values for summary data
    cstor2(:,n)=cwavetme;
    fstor2(:,n)=fotme;
    lstor2(:,n)=lptme;
    ostor2(:,n)=offtme;
    
    %convert the output amplitude values (mV) to single column
    cstorout=reshape(cstor,[],1);
    fstorout=reshape(fstor,[],1);
    lstorout=reshape(lstor,[],1);
    ostorout=reshape(ostor,[],1);
        
    %convert the output time values (min) to single column
    cstorout2=reshape(cstor2,[],1);
    fstorout2=reshape(fstor2,[],1);
    lstorout2=reshape(lstor2,[],1);
    ostorout2=reshape(ostor2,[],1);
    
    %duplicate mouse ID (name) for left (ch1) and right eye (ch2)
    tempname(n,1)=varNames3;
    tempname2 = repelem(tempname,2);
                            
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    n=n+1;
    
end

%%%%%%%%%%%%%%%%% Outside of the loop

%write to Excel DC-ERG Components Amplitude (mV) Summary  

varNames4=['Mouse_ID' varNames];
T4 = table(tempname2,cstorout,fstorout,lstorout,ostorout,'VariableNames',varNames4);
writetable(T4,Filename,'Sheet','Summary_mV');

%write to Excel DC-ERG Components Time-to-Peak (min) Summary
%Time-to-peak is defined from time of light onset (t=0 min)

varNames5=['Mouse_ID' varNames2];
T5 = table(tempname2,cstorout2,fstorout2,lstorout2,ostorout2,'VariableNames',varNames5);
writetable(T5,Filename,'Sheet','Summary2_min');

%____________________________________________________
% REMOVE EXCESS EXCEL SHEETS

excelFileName =Filename; 
excelFilePath = pwd; % Current working directory.
sheetName = 'Sheet'; % EN: Sheet, DE: Tabelle, etc. (Lang. dependent)

% Open Excel file.
objExcel = actxserver('Excel.Application');
objExcel.Workbooks.Open(fullfile(excelFilePath, excelFileName)); % Full path is necessary!

% Delete sheets.
try
% Throws an error if the sheets do not exist.
objExcel.ActiveWorkbook.Worksheets.Item([sheetName '1']).Delete;
objExcel.ActiveWorkbook.Worksheets.Item([sheetName '2']).Delete;
objExcel.ActiveWorkbook.Worksheets.Item([sheetName '3']).Delete;
catch
; % Do nothing.
end

% Save, close and clean up.
objExcel.ActiveWorkbook.Save;
objExcel.ActiveWorkbook.Close;
objExcel.Quit;
objExcel.delete;
%____________________________________________________

