Critical_reorganization_of_adolescent_prefrontal_circuitry/code/LFP/PAC/subfunctions/UserMetaScript.m

%% Provide overview of each detected oscillation in a specified channel

clear all
close all
experiments = get_experiment_list;
% experiments=experiments(1,randperm(length(experiments)));
% experiments=fliplr(experiments);
save_data=1;
repeatCalc=0;
%% subfunction 1 - getStimProperties

getStimulationProperties(experiments,save_data,repeatCalc) % insert StimPeriods in Excel Experiment List
correctStimulationProperties(experiments,save_data,repeatCalc)

%% Calculate basic features for baseline and stimulation (spiketime, spikephase, LFPpower... )
mainGetBaselineProperties(experiments,save_data,repeatCalc)
mainGetStimulationBasicFeatures(experiments,save_data,repeatCalc)
% TO-DO
% getStimulationFiringChannels(experiments,1,0)


%% Plotting and comparisons!
% for multiple experiments
% Baseline

expType={'B1'};
plotBaselinePowerArea(experiments,expType)
plotBaselineOscProperties(experiments,expType)

% plotBaselineCoherence(experiments,1)
% plotBaselinePower(experiments,1)
% plotBaselineOscProperties(experiments,1)
% % MUA
% plotBaselineMUAPhase(experiments,1,0)


% SUA

% Stimulation
% plotStimulationPower(experiments)
% plotRaster(experiment)


expType={'B1'};
plotStimulationMUAspikeProbabiliy(experiments,expType,1) % I AM WORKING
stimulusTypes={'ramp'};
plotStimulationMUAISI(experiments,expType,stimulusTypes) % I AM WORKING
firingRateSummaryForSPSS=plotStimulationMUAFiringRate(experiments,expType,stimulusTypes); % I AM WORKING
plotStimulationMUASpikePhase(experiments,expType,stimulusTypes) % I AM WORKING
plotStimulationPowerRatio(experiments,expType,stimulusTypes) % I AM WORKING


% SUA
plotStimulationBaselineSUASpikePhase(experiments,1)
plotStimulationSUASpikePhase(experiments,1,0)
plotStimulationSUAFiringRate(experiments,1)
PlotStimulationSUAISI(experiments,1)



% for single experiment
n_experiment=101;
expType={'B1'};
plotStimulationMUAspikeProbabiliy(experiments(n_experiment),expType,1);
% [DataWindowProbability]=plotStimulationMUAspikeProbabiliy(experiment,expType,1);
CSC=23;
stimulusTypes={'ramp'};
plotStimulationRasterPlot(experiments(n_experiment), stimulusTypes, CSC) 

% get responding animals/CSCs
clear all
experiments = get_experiment_list;
expType={'B1'};
countPositveAnimal=0;
for nExp = 1:length(get_experiment_list)
    experiment=experiments(nExp);
    if ~isempty(experiment.name) && strcmp(expType{1},experiment.Exp_type) && strcmp(experiment.IUEconstruct,'CAG-ChR2ET/TC-2AtDimer2')
        [CSCcalc,~,~,~] = getStimulationRespondingCSCsRamp(experiment,1,0);
        if CSCcalc >1
            countPositveAnimal=countPositveAnimal+1;
            positiveAnimal(countPositveAnimal)=nExp;
            CSCs.(['exp' num2str(nExp)]) = CSCcalc;
        end
    end
end

%% REVISION 2 (2016-11)
close all
clear all
experiments = get_experiment_list;
save_data=0;
repeatCalc=0;
resultsdir='Q:\Personal\Sebastian\ProjectOptoPFC\Analysis\Project3_HP_PFCopto\figures\AllCSCdeepLayers\';

for n_experiment=163%[101:115 117:132 140:150 152:165 167:172 174 176 177 179:181];
    expType={'B1'};
    stimulusType='square';
    Frequency=2;
    downsampling_factor=1;
    cutFrequencies=[500 5000];

    for CSC=17:32
        stimStructure.(strcat('CSC',num2str(CSC)))=getStimulusSignal(experiments(n_experiment),CSC, stimulusType, Frequency,downsampling_factor);
    end

    Window=stimStructure.CSC17.samplingrate*2.98:stimStructure.CSC17.samplingrate*3.08;
    
    
    h=figure(n_experiment);
    hold on
    plot(stimStructure.(strcat('CSC',num2str(CSC))).time(Window),stimStructure.(strcat('CSC',num2str(CSC))).signalD(1,Window)*50-16*100,'b')
    for CSC=17:32
        plot(stimStructure.(strcat('CSC',num2str(CSC))).time(Window),stimStructure.(strcat('CSC',num2str(CSC))).signal(:,Window)-CSC*100,'k')
        for sweep=1:30
            filteredSignal(sweep,:)=ZeroPhaseFilterZeroPadding(stimStructure.(strcat('CSC',num2str(CSC))).signal(sweep,:),stimStructure.(strcat('CSC',num2str(CSC))).samplingrate,cutFrequencies);
        end
        plot(stimStructure.(strcat('CSC',num2str(CSC))).time(Window),filteredSignal(:,Window)-CSC*100,'b')

        plot(stimStructure.(strcat('CSC',num2str(CSC))).time(Window),mean(stimStructure.(strcat('CSC',num2str(CSC))).signal(:,Window))-CSC*100,'r')
        plot(stimStructure.(strcat('CSC',num2str(CSC))).time(Window),mean(filteredSignal(:,Window))-CSC*100,'g')
    axis tight
    title(num2str(n_experiment))
    end
    savefig(h,strcat(resultsdir,'allCSC',num2str(n_experiment),'.fig'))
    close all
end


%% REVISION 2 (2016-11)
close all
clear all
experiments = get_experiment_list;
save_data=0;
repeatCalc=0;
resultsdir='Q:\Personal\Sebastian\ProjectOptoPFC\Analysis\Project3_HP_PFCopto\figures\SpikeDetectionTest\';

n_experiment=163
expType={'B1'};
stimulusType='square';
Frequency=2;
downsampling_factor=1;
cutFrequencies=[500 5000];


for CSC=26%17:32%27:28%17:32
    stimStructure.(strcat('CSC',num2str(CSC)))=getStimulusSignal(experiments(n_experiment),CSC, stimulusType, Frequency,downsampling_factor);
    
    for sweep=1:30
        stimStructure.(strcat('CSC',num2str(CSC))).filteredSignal(sweep,:)=ZeroPhaseFilterZeroPadding(stimStructure.(strcat('CSC',num2str(CSC))).signal(sweep,:),stimStructure.(strcat('CSC',num2str(CSC))).samplingrate,cutFrequencies);
    end
    
    
    h=figure;
    hold on
%     Window=stimStructure.(strcat('CSC',num2str(CSC))).samplingrate*2.98:stimStructure.(strcat('CSC',num2str(CSC))).samplingrate*3.08;
    Window=1:size(stimStructure.(strcat('CSC',num2str(CSC))).time,2);
    plot(stimStructure.(strcat('CSC',num2str(CSC))).time(Window),stimStructure.(strcat('CSC',num2str(CSC))).signalD(1,Window)*50,'b')
    for sweep=1:30
        time=stimStructure.(strcat('CSC',num2str(CSC))).time;
        filteredSignal=stimStructure.(strcat('CSC',num2str(CSC))).filteredSignal(sweep,:);
        
        plot(time(Window),filteredSignal(Window)-sweep*100,'k')
        
        % amplitudeDown threshold spike detection
        thr = std(filteredSignal)*3.5;
        [peakLocAmplDown, ~] = peakfinderOpto(filteredSignal,(max(filteredSignal)-min(filteredSignal))/4 ,-thr,-1,0);
        plot([time(Window(1)) time(Window(end))],[-thr -thr]-sweep*100,'k')
        if ~isempty(peakLocAmplDown)
            plot(time(peakLocAmplDown),-sweep*100-10,'y*')
        end
        
         % amplitudeUP threshold spike detection
        thr = std(filteredSignal)*3.5;
        [peakLocAmplUp, ~] = peakfinderOpto(filteredSignal,(max(filteredSignal)-min(filteredSignal))/4  ,thr,1,0);
        plot([time(Window(1)) time(Window(end))],[thr thr]-sweep*100,'k')
        if ~isempty(peakLocAmplDown)
            plot(time(peakLocAmplUp),-sweep*100+10,'y*')
        end
        
        % slope threshold spike detection
        dv=[0 diff(filteredSignal)];
        plot(time(Window),dv(Window)-sweep*100-50,'g')
        thr = std(dv)*3.0;
        [peakLocSlope, ~] = peakfinderOpto(dv,(max(dv)-min(dv))/4  ,-thr,-1,1);
        plot([time(Window(1)) time(Window(end))],[-thr -thr]-sweep*100-50,'g')
        if ~isempty(peakLocSlope)
            plot(time(peakLocSlope),-sweep*100-50,'b*')
        end
        
        % combined slope + amplitude
        spikeWindow=round(stimStructure.(strcat('CSC',num2str(CSC))).samplingrate*0.001);
        peakLocCombined=[];
        for f1=1:length(peakLocAmplDown)
            if any(ismember(peakLocAmplDown(f1)-spikeWindow:peakLocAmplDown(f1),peakLocSlope))
                if any(ismember(peakLocAmplDown(f1)-spikeWindow:peakLocAmplDown(f1),peakLocAmplUp))
                    peakLocCombined=[peakLocCombined peakLocAmplDown(f1)];
                end
            end
        end
        if ~isempty(peakLocCombined)
            plot(time(peakLocCombined),-sweep*100-25,'r*')
        end
        
        
%         figure
%         hold on
%         plotspikeWindow=round(stimStructure.(strcat('CSC',num2str(CSC))).samplingrate*0.005);
%         for f1=1:length(peakLocCombined)
%             plot((-plotspikeWindow:plotspikeWindow)/stimStructure.(strcat('CSC',num2str(CSC))).samplingrate,filteredSignal(peakLocCombined(f1)-plotspikeWindow:peakLocCombined(f1)+plotspikeWindow))
%         end
%         axis tight




%         % Nonlinear Energy Operator (NEO) for spike detection
%         vNEO=filteredSignal(2:end-1).*filteredSignal(2:end-1)-filteredSignal(3:end).*filteredSignal(1:end-2);
%         vNEO=[0 vNEO 0];
%         plot(time(Window),vNEO(Window)-sweep*100-50,'b')
%         
%         thrNEO = std(vNEO)*10;
%         [peakLocThresholdNEO, ~] = peakfinderOpto(vNEO,0 ,thrNEO,1);
%         plot([time(Window(1)) time(Window(end))],[thrNEO thrNEO]-sweep*100-50,'b')
%         if ~isempty(peakLocThresholdNEO)
%             plot(time(peakLocThresholdNEO),-sweep*100-50,'r*')
%         end
        
    end
    axis tight
    title(strcat(num2str(n_experiment),'CSC',num2str(CSC)))
%     savefig(h,strcat(resultsdir,'exp',num2str(n_experiment),'CSC',num2str(CSC),'.fig'))
%     close all
end