Upload files to 'code/LFP'

code/LFP/PAC/MainSinglePAC.m

@@ -0,0 +1,40 @@
+function Analysis= MainSinglePAC(rangePhase,rangeAmplitude,measure,Merge_num,Signal,fs)
+
+% This is a helperfunction to calculate a matrix of PAC values using either the ESC, MI(MVL) 
+% or CFC measure. 
+
+% INPUTS:
+%  fs - sampling frequency;
+%  Signal - vector of 0.5-100 Hz filtered LFP signal
+%  measure - measure to be used - it should be: 'esc', 'mi' or 'cfc';
+%  rangePhase - range of frequencies for low signal; for example 1:1:20
+%  rangeAmplitude - range of frequencies for high signal; for example 30:2:200
+%  Merge_num - used 50
+
+% Output
+%PAC results for each analyzed Signal window, including:
+    %pac.pacmat=pacmat;
+    %pac.freqvec_ph---the phase frequency band;
+    %pac.freqvec_amp---the amplitude frequency band;
+    %pac.shf_data_mean---shuffle mean;
+    %pac.shf_data_std----shuffle std;
+    %pac.relat_mi-----Zscore of the pac.pacmat
+
+% Author: Angela Onslow, May 2010
+%% calc PAC
+segs=size(Signal,1);
+Merge_segs=floor(2*segs/Merge_num-1);
+for Merge_seg=1:Merge_segs
+    Tem=floor(Merge_seg-1)*Merge_num/2+1:(Merge_seg+1)*Merge_num/2;
+    X1=Signal(Tem,:);
+    X2=Signal(Tem,:);
+    
+    [PAC]=find_pac_shf_several_CutGlue_segs(X2,fs,measure,X1,rangePhase,rangeAmplitude);
+    
+    
+    pacwin{Merge_seg}=PAC;
+    clear PAC
+end
+Analysis.pacwin=pacwin;
+clear Signal pacwin
+

code/LFP/PAC/subfunctions/PACmeanALL.m

@@ -0,0 +1,111 @@
+function []=PACmeanALL(Param,Path,Experiment,exp_num)
+
+for inum= 1: length(exp_num)
+    iexperiment=exp_num(inum);
+    
+filename1=Experiment(iexperiment).name;
+filename2=Experiment(iexperiment).name2;
+
+
+    
+    
+load( strcat(Path.output,filesep,'Mi_PAC OB Phase LEC Amplitude 1s no zeropad filter 1',filesep,filename1,'.mat'));
+Phases(inum,:)=pacwin{1}.freqvec_ph;
+Amplitudes(inum,:)=pacwin{1}.freqvec_amp;
+for iseg=1:length(pacwin)
+Miseg(iseg,:,:)=pacwin{iseg}.relat_mi;
+end
+Mi(inum,:,:)=mean(Miseg,1);
+
+clear pacwin
+load( strcat(Path.output,filesep,'Mi_PAC OB Phase LEC Amplitude 1s no zeropad filter 1',filesep,filename2,'.mat'));
+Phases_nO(inum,:)=pacwin{1}.freqvec_ph;
+Amplitudes_nO(inum,:)=pacwin{1}.freqvec_amp;
+for iseg=1:length(pacwin)
+Miseg_nO(iseg,:,:)=pacwin{iseg}.relat_mi;
+end
+size(Miseg_nO)
+Mi_nO(inum,:,:)=mean(Miseg_nO,1);
+size(Mi_nO)
+
+Mi_sub(inum,:,:)=Mi(inum,:,:)-Mi_nO(inum,:,:);
+Mi_div(inum,:,:)=Mi(inum,:,:)./Mi_nO(inum,:,:)
+
+figure
+subplot(2,2,1)
+imagesc(Phases(inum,:),Amplitudes(inum,:),squeeze(Mi(inum,:,:)));axis xy;  
+title('Baseline PAC OB-LEC')
+xlabel('OB phase frequency [Hz]')
+ylabel('LEC amplitude frequency [Hz]')
+subplot(2,2,2)
+imagesc(Phases_nO(inum,:),Amplitudes_nO(inum,:),squeeze(Mi_nO(inum,:,:)));axis xy;  
+title('naris occlusion PAC OB-LEC')
+xlabel('OB phase frequency [Hz]')
+ylabel('LEC amplitude frequency [Hz]')
+colormap(jet)
+subplot(2,2,3)
+imagesc(Phases(inum,:),Amplitudes(inum,:),squeeze(Mi_sub(inum,:,:)));axis xy;  
+title(' PAC OB-LEC baseline - nO')
+xlabel('OB phase frequency [Hz]')
+ylabel('LEC amplitude frequency [Hz]')
+subplot(2,2,4)
+imagesc(Phases_nO(inum,:),Amplitudes_nO(inum,:),squeeze(Mi_div(inum,:,:)));axis xy;  
+title('PAC OB-LEC baseline/nO')
+xlabel('OB phase frequency [Hz]')
+ylabel('LEC amplitude frequency [Hz]')
+colormap(jet)
+
+
+end
+
+MeanMi=mean(Mi,1);
+MeanMi_nO=mean(Mi_nO,1);
+MeanMi_sub=mean(Mi_sub,1);
+MeanMi_div=mean(Mi_div,1);
+
+%permutation testing
+c = Mi;
+d = Mi_nO;
+
+
+npermutes = 1000;
+
+diff = mean(c-d);
+diff_null = zeros(npermutes,1);
+for permi = 1:npermutes
+    anull = [c; d];
+    ind = randperm(length(anull),length(anull)/2);
+    bnull = anull(ind);
+    anull(ind) = [];
+    
+    diff_null(permi) = mean(anull-bnull);
+end
+
+pMi = 1-sum(abs(diff)>=abs(diff_null))/npermutes
+
+
+
+figure
+subplot(2,2,1)
+imagesc(Phases(1,:),Amplitudes(1,:),squeeze(MeanMi));axis xy;  
+title('mean Baseline PAC OB-LEC')
+xlabel('OB phase frequency [Hz]')
+ylabel('LEC amplitude frequency [Hz]')
+subplot(2,2,2)
+imagesc(Phases_nO(1,:),Amplitudes_nO(1,:),squeeze(MeanMi_nO));axis xy; 
+title('naris occlusion PAC OB-LEC')
+xlabel('OB phase frequency [Hz]')
+ylabel('LEC amplitude frequency [Hz]')
+colormap(jet)
+subplot(2,2,3)
+imagesc(Phases(1,:),Amplitudes(1,:),squeeze(MeanMi_sub));axis xy;  
+title('mean PAC OB-LEC baselin-nO')
+xlabel('OB phase frequency [Hz]')
+ylabel('LEC amplitude frequency [Hz]')
+subplot(2,2,4)
+imagesc(Phases_nO(1,:),Amplitudes_nO(1,:),squeeze(MeanMi_div));axis xy; 
+title('mean PAC OB-LEC baseline/nO')
+xlabel('OB phase frequency [Hz]')
+ylabel('LEC amplitude frequency [Hz]')
+colormap(jet)
+end

code/LFP/PAC/subfunctions/UMA_PAC.m

@@ -0,0 +1,25 @@
+%% UMA PAC
+close all
+clear all
+exp_num=[2,3,4,5,8,9,11,12,14,15];
+Experiment=get_experiment_list_OB_LEC_nose_closed;
+Path=get_path;
+Param=get_parameters;
+
+
+%PAC_TORTetal_main(Param,Path,Experiment,exp_num)  %Maltes Method
+%meanPAC(Param,Path,Experiment,exp_num)
+
+
+
+range1=[2:1:12];
+range2=[6:2:80];
+
+cross='ba'; 
+measure='mi';
+Merge_num=40;
+
+
+main_function_calculate_PAC_several_GutGlue_segs(Experiment,exp_num, Path,range1,range2,cross,measure,Merge_num);
+
+PACmeanALL(Param,Path,Experiment,exp_num)

code/LFP/PAC/subfunctions/UserMetaScript.m

@@ -0,0 +1,235 @@
+%% 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
+
+

code/LFP/PAC/subfunctions/ampvec.m

@@ -0,0 +1,30 @@
+function y = ampvec(f,s,Fs,width)
+% function y = ampvec(f,s,Fs,width)
+%
+% Returns a vector 'y' containing the instantaneous amplitude values of 
+% signal 's' filtered for frequency 'f' (via convolution with a complex 
+% Morlet wavelet described by 'width')
+%
+% INPUTS:
+% f - frequency to filter at 
+% s - signal to be filtered
+% Fs - sampling frequency of s
+% width - parameter which defines the mother wavelet (this is then 
+% scaled and translated in order to filter for different frequencies, 
+% >= 5 is suggested, see Tallon-Baudry et al., J. Neurosci. 15, 722-734 
+% (1997))
+%
+% NOTE: This function is a modification of code written Ole Jensen, August
+% 1998, see ENERGYVEC
+%
+% Author: Angela Onslow, May 2010
+
+dt = 1/Fs;
+sf = f/width;
+st = 1/(2*pi*sf); 
+
+t=-(width/2)*st:dt:(width/2)*st; 
+m = morlet_wav(f,t,width);
+y = conv(s,m');
+y = abs(y);
+y = y(ceil(length(m)/2):length(y)-floor(length(m)/2));

code/LFP/PAC/subfunctions/bpvec.m

@@ -0,0 +1,30 @@
+function y = bpvec(f,s,Fs,width)
+% function y = bpvec(f,s,Fs,width)
+%
+% Returns a vector 'y' containing signal 's' filtered for frequency 'f'
+% (via convolution with a complex Morlet wavelet described by 'width')
+%
+% INPUTS:
+% f - frequency to filter at 
+% s - signal to be filtered
+% Fs - sampling frequency of s
+% width - parameter which defines the mother wavelet (this is then 
+% scaled and translated in order to filter for different frequencies, 
+% >= 5 is suggested, see Tallon-Baudry et al., J. Neurosci. 15, 722-734 
+% (1997))
+%
+% NOTE: This function is a modification of code written Ole Jensen, August
+% 1998, see ENERGYVEC
+%
+% Author: Angela Onslow, May 2010
+
+dt = 1/Fs;
+sf = f/width;
+st = 1/(2*pi*sf); 
+
+t=-(width/2)*st:dt:(width/2)*st; 
+m = morlet_wav(f,t,width);
+y = conv(s,m);
+y = real(y);
+y = y(ceil(length(m)/2):length(y)-floor(length(m)/2));
+

code/LFP/PAC/subfunctions/esc_measure.m

@@ -0,0 +1,50 @@
+function escval = esc_measure(ph_sig, amp_sig, avg)
+% function escval = esc_measure(ph_sig, amp_sig, avg)
+%
+% Returns a value (or vector of values) for the ESC measure calculated
+% between two signals. Signals may contain one of more trials. Multiple
+% trials may be averaged so as to return one ESC value or a vector of
+% ESC values calculated for each trial may be returned, depending on the 
+% 'avg' argument. Signals should be passed as column vectors, multiple 
+% trials stored as multiple columns.
+%
+% INPUTS:
+% ph_sig - signal filtered for a lower, modulating frequency (e.g. theta
+% band oscillations)
+%
+% amp_sig - signal filtered for a higher, modulated frequency (e.g. gamma
+% band oscillations)
+%
+% avg - string, either 'y' or 'n', determines whether ESC values are
+% averaged over trials or returned as a vector
+%
+% Author: Angela Onslow, May 2010
+
+escsum = 0;
+if size(ph_sig, 2) ~= size(amp_sig, 2)
+    sprintf('Error - Signals must have the same number of trials')
+    return
+end
+num_trials = size(ph_sig, 2);
+
+if strcmp(avg, 'y')
+    
+    %Average over trials using the Fisher transform
+    for c = 1:num_trials
+            
+            r = corrcoef(ph_sig(:,c), amp_sig(:,c));
+            escsum = escsum + atanh(r(1,2));
+            
+    end
+
+    escsum = escsum/num_trials;
+    escval = tanh(escsum);
+    
+else
+    escval = zeros(num_trials,1);
+
+    for i = 1:num_trials
+        r = corrcoef(ph_sig(:,i), amp_sig(:,i));
+        escval(i,1) = r(1,2);
+    end
+end

code/LFP/PAC/subfunctions/filt_signalsWAV_several_CutGlue_segs.m

@@ -0,0 +1,101 @@
+function [ph_filt_signals, amp_filt_signals] = filt_signalsWAV_several_CutGlue_segs(sig1,sig2,...
+    Fs, ph_freq_vec, amp_freq_vec, measure, width)
+% function [ph_filt_signals, amp_filt_signals] = filt_signalsWAV(sig1,sig2,...
+% Fs, ph_freq_vec, amp_freq_vec, measure, width)
+%
+% Returns cell arrays 'ph_filt_signals' and 'amp_filt_signals'. These 
+% contain 'sig1' and 'sig2' correctly filtered in order to calculate the
+% PAC measure, either ESC or MI given by 'measure'. The frequencies
+% to filter at are determined by the vectors 'ph_freq_vec' and
+% 'amp_freq_vec'. The frequency range is defined by the minimum and the
+% maximum value and the interval, or bandwith, between values. For example
+% if 'ph_freq_vec' is defined as 1:5:101 ([1, 6, 11, 16, 21...101]) then
+% the centre frequencies that will be filtered for are [3, 8, 13, 18...98]
+% i.e. the centre of the 5 Hz bins from 1 Hz to 101 Hz. Filtering is
+% achieved via convolution with a complex Morlet wavelet (number of cycles 
+% defined by 'width' parameter). 
+%
+% INPUTS:
+% sig1 - signal which will be analysed as containing the higher frequency,
+% modulated PAC signal
+% sig2 - signal which will be analysed as containing the lower frequency,
+% modulating signal
+% Fs - sampling frequency, in Hz
+% ph_freq_vec - vector of frequencies to filter sig2 for, in Hz
+% amp_freq_vec - vector of frequencies to filter sig1 for, in Hz
+% measure - PAC measure which will be calculated using these filtered
+% signals
+% width - number of cycles which will define the mother Morlet wavelet
+%
+% OUTPUTS:
+% ph_filt_signals - 1 x (number of bins determined by ph_freq_vec) cell 
+% array, each element has the same dimensions as the original signals
+% amp_filt_signals - 1 x (number of bins determined by amp_freq_vec) cell 
+% array, each element has the same dimensions as the original signals
+%
+% Author: Angela Onslow, May 2010
+
+total_num_dp = size(sig1,1);  %  length of data
+num_trials = 1:size(sig1,2);  %  1:sample numbers
+
+phfreq_low = min(ph_freq_vec);
+phfreq_high = max(ph_freq_vec);
+phfreq_bw = diff(ph_freq_vec(1:2));
+ampfreq_low = min(amp_freq_vec);
+ampfreq_high = max(amp_freq_vec);
+ampfreq_bw = diff(amp_freq_vec(1:2));
+
+xbins = ceil((phfreq_high - phfreq_low)/phfreq_bw);  % number of frequency bands of sig2
+ybins = ceil((ampfreq_high - ampfreq_low)/ampfreq_bw);  % number of frequency bands of sig1
+
+
+% Create structures to store filtered signals
+ph_filt_signals = cell(1,xbins);
+amp_filt_signals = cell(1,ybins);
+
+for i2=1:xbins
+    ph_filt_signals{1,i2} = zeros(total_num_dp*max(num_trials),1);
+end
+
+
+for i2=1:ybins
+    amp_filt_signals{1,i2} = zeros(total_num_dp*max(num_trials),1);
+end
+
+% Filter and store filtered signals %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ====  Filter sig2 for phase freq  ====
+for i2=1:xbins
+    
+    upper_bin = phfreq_low+(i2*phfreq_bw);
+    lower_bin = upper_bin-phfreq_bw;
+ 
+    if strcmp(measure, 'esc')
+        for i3=1:max(num_trials)
+            tem=(i3-1)*total_num_dp+1:(i3)*total_num_dp;
+            ph_filt_signals{1,i2}(tem,1) = bpvec((lower_bin + floor((upper_bin- lower_bin)/2)), sig2(:,i3), Fs, width);
+        end
+    else
+        for i3=1:max(num_trials)
+            tem=(i3-1)*total_num_dp+1:(i3)*total_num_dp;
+            aa=lower_bin + floor((upper_bin- lower_bin)/2);
+            ph_filt_signals{1,i2}(tem,1) = phasevec((lower_bin + floor((upper_bin- lower_bin)/2)),sig2(:,i3), Fs,width);
+        end
+    end
+    
+end
+
+% ====  Filter sig1 for amplitude freq  ====
+for i2=1:ybins
+    
+    upper_bin = ampfreq_low+(i2*ampfreq_bw);
+    lower_bin = upper_bin-ampfreq_bw;
+        
+    for i3=1:max(num_trials)
+        tem=(i3-1)*total_num_dp+1:(i3)*total_num_dp;
+        amp_filt_signals{1,i2}(tem,1) = ampvec((lower_bin + floor((upper_bin- lower_bin)/2)),sig1(:,i3), Fs, width);
+
+    end
+    
+  
+end

code/LFP/PAC/subfunctions/find_pac_nofilt.m

@@ -0,0 +1,175 @@
+function [pacmat, freqvec_ph, freqvec_amp] = find_pac_nofilt (sig_pac, Fs,...
+    measure, sig_mod, ph_freq_vec, amp_freq_vec, plt, waitbar, width, nfft,...
+        dataname, sig_pac_name, sig_mod_name)
+%function [pacmat, freqvec_ph, freqvec_amp] = find_pac_nofilt (sig_pac, Fs,...
+%   measure, sig_mod, ph_freq_vec, amp_freq_vec, plt, waitbar, width, nfft,...
+%       dataname, sig_pac_name, sig_mod_name)
+%
+% This function calculates a matrix of PAC values using either the ESC, MI 
+% or CFC measure. 
+% It assumes that the input is a prefiltered signal since this function 
+% does not include any filtering. The output is a matrix of PAC values and
+% a plot of this matrix (depending on the value of the 'plt' argument).
+% This signal can only take single vector signals as input but the
+% provision for multiple trials will be added (partly implemented for the
+% CFC measure)
+%
+% Basic function call:
+% function pacmat = find_pac_nofilt(sig_pac, Fs, measure)
+%
+% REQUIRED INPUTS:
+%  sig_pac - signal suspected of containing PAC
+%  Fs - sampling frequency
+%  measure - measure to be used - it should be: 'esc', 'mi' or 'cfc'
+%
+% The function can be executed with many optional inputs:
+%  function pacmat = find_pac_nofilt(sig_pac, Fs, measure, ...
+%    sig_mod, ph_freq_vec, amp_freq_vec, plt, width, nfft, dataname,...
+%       sig_pac_name, sig_mod_name)
+%
+% OPTIONAL INPUTS:
+%  sig_mod - signal containing modulating frequency; default = sig_pac
+%  ph_freq_vec - range of frequencies for modulating signal; default = 1:5:101
+%  amp_freq_vec - range of frequencies for PAC signal; default = 1:5:101
+%  plt - flag indicating whether the output should be plotted - it should
+%        be 'y' or 'n'; default = 'y'
+%  width - width of the wavelet filter; default = 7
+%  nfft - the number of points in fft; default = 200
+%  dataname - the name of the dataset to be included as a graph title;
+%             default = ''
+%  sig_pac_name - the name of sig_pac to be printed on the y-axis; default = ''
+%  sig_mod_name - the name of sig_mod to be printed on the x-axis; default = ''
+%
+% Author: Angela Onslow, May 2010 
+
+% Checks of input variables %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if nargin < 4
+    sig_mod = sig_pac;
+end
+if nargin < 5
+    ph_freq_vec = 1:5:101;
+end
+if nargin < 6
+    amp_freq_vec = 1:5:101;
+end
+if nargin < 7
+    plt = 'y';
+end
+if nargin < 8 
+    waitbar = 0;
+end
+if nargin < 9
+    width = 7;
+end
+if nargin < 10
+    nfft = ceil(Fs/(diff(ph_freq_vec(1:2))));
+end
+if nargin < 11
+    dataname = '';
+end
+if nargin < 12
+    sig_pac_name = '';
+end
+if nargin < 13
+    sig_mod_name = '';
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Set up some parameters for clarity/storage
+xbins = ceil((max(ph_freq_vec) - min(ph_freq_vec))/(diff(ph_freq_vec(1:2))));
+ybins = ceil((max(amp_freq_vec) - min(amp_freq_vec))/(diff(amp_freq_vec(1:2))));
+cent_freq_vec = zeros(xbins,1);
+cent_freq_vec2 = zeros(ybins, 1);
+
+
+if isa(sig_pac, 'cell')
+    % Data represents filtered signals
+    num_trials = size(sig_pac{1,1},2); 
+else
+    % Data represents signals 
+    num_trials = size(sig_pac,2);
+    
+end
+
+if waitbar == 1
+    counter = 0;
+    countermax = xbins*ybins;
+    fprintf('\nCalculating PAC values\n');
+end
+
+for i =1:xbins
+    upper_bin = min(ph_freq_vec)+i*(diff(ph_freq_vec(1:2)));
+    lower_bin = upper_bin-(diff(ph_freq_vec(1:2)));
+    cent_freq_vec(i) =  lower_bin + floor((upper_bin- lower_bin)/2);
+end
+
+
+for i =1:ybins
+    upper_bin = min(amp_freq_vec)+i*(diff(amp_freq_vec(1:2)));
+    lower_bin = upper_bin-(diff(amp_freq_vec(1:2)));
+    cent_freq_vec2(i) =  lower_bin + floor((upper_bin- lower_bin)/2);
+end
+
+
+freqvec_amp = cent_freq_vec2;
+freqvec_ph = cent_freq_vec;
+
+
+% Calculate PAC measures %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if (strcmp(measure, 'esc') || strcmp(measure, 'mi'))
+pacmat = zeros(ybins, xbins);    
+    
+    for i = 1:ybins
+        for j = 1:xbins
+            % Calculate matrix of PAC values
+            if strcmp(measure, 'esc')
+                pacmat(i,j) = esc_measure(sig_mod{1,j}, ...
+                    sig_pac{1,i}, 'y');
+            end
+            
+            if strcmp(measure, 'mi')
+                % Pacmat full of raw mi values, not yet normalized
+                pacmat(i,j) = mi_measure(sig_mod{1,j}, sig_pac{1,i});
+            end
+            
+            % Display current computational step to user
+            if waitbar ==1
+                counter = counter+1;
+                if counter == 1
+                    fprintf('%03i%% ', floor((counter/countermax)*100));
+                else
+                    fprintf('\b\b\b\b\b%03i%% ', floor((counter/countermax)*100));
+                end
+                if counter == countermax
+                    fprintf('\n');
+                end
+            end
+        end
+    end
+    
+end
+
+
+
+if strcmp(measure, 'cfc')
+    
+    % Calculate matrix of PAC values (these have been averaged over trials)
+    [pacmat, freqvec_ph, freqvec_amp] = cfc_measure(sig_pac, ...
+             sig_mod, 'y',ph_freq_vec, freqvec_amp, Fs, nfft, width, waitbar);
+              
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if strcmp(plt, 'y')
+    pac_plot_fun(pacmat, freqvec_ph, freqvec_amp, measure, sig_pac_name,...
+                 sig_mod_name, dataname, sig_pac, Fs, sig_mod);
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+     
+ 
+

code/LFP/PAC/subfunctions/find_pac_shf_several_CutGlue_segs.m

@@ -0,0 +1,142 @@
+function pac = find_pac_shf_several_CutGlue_segs (SIG_pac, Fs, measure, ...
+    SIG_mod, ph_freq_vec, amp_freq_vec)
+% Xiaxia 
+% This function calculates a matrix of PAC values using either the ESC, MI(MVL) 
+% or CFC measure. 
+% It uses shuffled datasets to conduct a significance analysis of the PAC
+% values found
+
+% INPUTS:
+%  sig_pac - high frequency,several cut-glue row signals; 
+%  Fs - sampling frequency;
+%  measure - measure to be used - it should be: 'esc', 'mi' or 'cfc';
+%  SIG_mod - low frequency,several cut-glue row signal;
+%  ph_freq_vec - range of frequencies for low signal; for example 1:1:20
+%  amp_freq_vec - range of frequencies for high signal; for example
+%  30:2:200
+
+% Output
+%pac.pacmat=pacmat;
+%pac.freqvec_ph---the phase frequency band;
+%pac.freqvec_amp---the amplitude frequency band;
+%pac.shf_data_mean---shuffle mean;
+%pac.shf_data_std----shuffle std;
+%pac.relat_mi-----Zscore of the pac.pacmat
+
+% Author: Angela Onslow, May 2010
+
+% set some paremeters: 
+   waitbar = 1;%  waitbar - display progress in the command window; suggest 1 (Xiaxia)
+   width = 7;%  width - width of the wavelet filter; sugguest 7 (Xiaxia)
+   nfft = ceil(Fs/(diff(ph_freq_vec(1:2)))); %  nfft - the number of points in fft; 
+   num_shf = 50; %  num_shf - the number of shuffled data sets to use during significancetesting; suggest 50 (Xiaxia)
+   alpha = 0.05; %  alpha - significance value to use; default = 0.05 
+    
+% Set up some parameters for clarity
+xbins = ceil((max(ph_freq_vec) - min(ph_freq_vec))/(diff(ph_freq_vec(1:2))));
+ybins = ceil((max(amp_freq_vec) - min(amp_freq_vec))/(diff(amp_freq_vec(1:2))));
+alpha = alpha/(xbins*ybins); % Uncomment to use Bonferonni Correction
+    
+   sig_pac=SIG_pac';
+   sig_mod=SIG_mod';
+    
+% Filtering phase and amplitude and glue together to obtain a long dataset
+if (strcmp(measure, 'esc')) ||(strcmp(measure, 'mi')) 
+[filt_sig_mod, filt_sig_pac] =filt_signalsWAV_several_CutGlue_segs(sig_pac, sig_mod, Fs, ...
+    ph_freq_vec, amp_freq_vec, measure, width);
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+% Create shuffled datasets and distribution of PAC values %%%%%%%%%%%%%%%%%
+if num_shf ~= 0
+if waitbar == 1
+    fprintf('\nCreating shuffled data sets\n');
+end
+for s = 1:num_shf
+    
+    if strcmp(measure, 'esc')
+           
+        shuffled_sig_amp = shuffle_esc(filt_sig_pac, Fs);
+        shf_pacmat_final(s,:,:) = find_pac_nofilt(shuffled_sig_amp, Fs, measure, filt_sig_mod, ph_freq_vec, amp_freq_vec,'n');
+        
+    end
+     
+
+    if strcmp(measure, 'mi')
+
+        shuffled_sig_amp = shuffle_esc(filt_sig_pac, Fs);
+        shf_pacmat_final(s,:,:) = find_pac_nofilt(shuffled_sig_amp, Fs, measure, filt_sig_mod, ph_freq_vec, amp_freq_vec,'n');
+        
+    end
+    
+    if strcmp(measure, 'cfc')
+          
+        shuffled_sig1 = shuffle_esc(sig_pac, Fs);
+        shf_pacmat_final(s,:,:) = find_pac_nofilt(shuffled_sig1, Fs,measure, sig_mod, ph_freq_vec, amp_freq_vec,'n', 0, width, nfft);
+        
+    end
+    
+    % Display current computational step to user
+    if waitbar == 1
+        if s == 1
+            fprintf('%03i%% ', floor((s/num_shf)*100));
+        else
+            fprintf('\b\b\b\b\b%03i%% ', floor((s/num_shf)*100));
+        end
+        if s == num_shf
+            fprintf('\n');
+        end
+    end
+end
+
+%Find mean and standard deviation of shuffled data sets
+if strcmp(measure, 'mi')
+    for i =1:ybins
+        for j=1:xbins
+            [shf_data_mean(i,j), shf_data_std(i,j)] = normfit(shf_pacmat_final(:,i,j));
+        end
+    end
+    
+else
+    shf_data_mean = squeeze (mean (shf_pacmat_final, 1));
+    shf_data_std = squeeze (std (shf_pacmat_final, 1));
+end
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Calculate PAC measures %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if strcmp(measure, 'esc')
+    [pacmat, freqvec_ph, freqvec_amp] = find_pac_nofilt(filt_sig_pac, Fs, measure, filt_sig_mod, ph_freq_vec, amp_freq_vec, 'n', 1);
+end
+
+if strcmp(measure, 'mi')
+    [pacmat, freqvec_ph, freqvec_amp] = find_pac_nofilt(filt_sig_pac, Fs, measure, filt_sig_mod, ph_freq_vec, amp_freq_vec, 'n', 1);
+end
+
+if strcmp(measure, 'cfc')
+    [pacmat, freqvec_ph, freqvec_amp] = find_pac_nofilt(sig_pac, Fs, measure, sig_mod, ph_freq_vec, amp_freq_vec, 'n', 1, width, nfft);
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Compute significance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if num_shf ~= 0 
+for i = 1:size(pacmat,1)
+    for j = 1:size(pacmat,2)
+        [h, p] = my_sig_test(pacmat(i,j), squeeze(shf_pacmat_final(:,i,j)), alpha);
+        if h == 0
+            pacmat(i,j) = 0;
+        end
+    end
+end
+pac.shf_data_mean=shf_data_mean;
+pac.shf_data_std=shf_data_std;
+end
+
+
+pac.pacmat=pacmat;
+pac.freqvec_ph=freqvec_ph;
+pac.freqvec_amp=freqvec_amp;
+pac=relat_shaf(pac);
+

code/LFP/PAC/subfunctions/main_function_calculate_PAC_several_GutGlue_segs.m

@@ -0,0 +1,93 @@
+function pacwin= main_function_calculate_PAC_several_GutGlue_segs(Experiment,exp_num, Path,rangePhase,rangeAmplitude,cross,measure,Merge_num)
+%written by Xiaxia
+% cross --- 'aa' , ab, bb, ba , region
+% measure---- the method name, 'mi', 'esc', 'cfc'
+% rangePhase - range of frequencies for low signal; for example 1:1:20
+% rangeAmplitude - range of frequencies for high signal; for example
+%  30:2:200
+% Merge_num, the number of segments to merge
+   cd(Path.output)
+     dirName=strcat('Mi_PAC OB Phase LEC Amplitude 2s no zeropad filter 1')
+     mkdir(dirName)
+
+
+
+th_r=2;flag_signal='xndetrend';
+for inum=1:length(exp_num)
+    iExperiment=exp_num(inum);
+    nameBaseline=Experiment(iExperiment).name1;
+    name_nO=Experiment(iExperiment).name2;
+    Name={nameBaseline,name_nO};
+    Channel1=Experiment(iExperiment).OBchannels; %OB
+    Channel2=Experiment(iExperiment).LECchannels; %LEC
+    
+    %%%%%%%%% load filtered and cut glued signal to calculate PE
+
+for group=1:length(Name)
+filename=Name{group};
+    
+clear signal_LFP
+clear Signal1
+load( strcat(Path.temp,filesep,'nlx_load_LFP_LP100',filesep,filename,'\CSC',num2str(Channel1),'.mat'));
+fs=fsOutput;
+Signal1=ZeroPhaseFilter(signal_LFP,fs,[1 100]);
+
+clear signal_LFP
+clear Signal2
+load( strcat(Path.temp,filesep,'nlx_load_LFP_LP100',filesep,filename,'\CSC',num2str(Channel2),'.mat'));
+Signal2=ZeroPhaseFilter(signal_LFP,fs,[1 100]);
+
+
+%load overlapping Event timestamps
+clear oscStartInner
+clear oscEndInner
+clear timestamps1  
+load( strcat(Path.output,filesep,'SymOsc50',filesep,filename,'.mat'));
+timestamps1(1,:)=oscStartInner/fs;
+timestamps1(2,:)=oscEndInner/fs;
+        
+[ Sig1,Nwindows1]=cutandglue_nozeropadding(Signal1,timestamps1',fs);
+[ Sig2,Nwindows2]=cutandglue_nozeropadding(Signal2,timestamps1',fs);
+Signal1=Sig1.xndetrend;
+size(Signal1)
+Signal2=Sig2.xndetrend;
+    
+    segs=size(Signal2,1)
+    Merge_segs=floor(2*segs/Merge_num-1)
+    for Merge_seg=1:Merge_segs
+        %Tem=(Merge_seg-1)*Merge_num+1:Merge_seg*Merge_num;
+        Tem=floor(Merge_seg-1)*Merge_num/2+1:(Merge_seg+1)*Merge_num/2;
+    if strcmp(cross,'aa')
+        X1=Signal1(Tem,:);
+        X2=Signal1(Tem,:);
+    end
+    if strcmp(cross,'ab')
+        X1=Signal1(Tem,:);
+        X2=Signal2(Tem,:);
+    end
+    if strcmp(cross,'ba')
+        X2=Signal2(Tem,:);
+        X1=Signal1(Tem,:);
+    end
+    if strcmp(cross,'bb')
+        X1=Signal2(Tem,:);
+        X2=Signal2(Tem,:);
+    end
+    
+    [PAC_2to12_15to80]=find_pac_shf_several_CutGlue_segs(X2,fs,measure,X1,rangePhase,rangeAmplitude);
+     
+
+    
+    pacwin{Merge_seg}=PAC_2to12_15to80
+        clear PAC_2to12_15to80;
+    end
+    %pacwin=relat_shaf_win(pacwin{iExperiment});
+    
+     dirName
+    cd(Path.output)
+    cd(dirName)
+    save(filename,'pacwin')
+end   
+end
+end
+

code/LFP/PAC/subfunctions/mi_measure.m

@@ -0,0 +1,32 @@
+ function mival = mi_measure(phase_sig, amp_sig)
+% function mival = mi_measure(phase_sig, amp_sig)
+%
+% Returns a value for the MI measure calculated between two signals.
+% (Functionality to deal with multiple trials will be added soon)
+%
+% INPUTS:
+%
+% phase_sig - the instantaneous phase values for a signal which has been
+% filtered for a lower, modulating frequency, passed as a column vector
+%
+% amp_sig - the amplitude values for a signal which has been filtered for a
+% higher, modulated frequency, passed as a column vector 
+%
+% Author: Angela Onslow, May 2010
+
+num_trials = size(phase_sig, 2);
+
+for count = 1:num_trials
+    
+    %Create composite signal
+    z = amp_sig(:,count).*exp(1i*phase_sig(:,count));
+    m_raw(count) = mean(z);  %Compute the mean length of composite signal.
+        
+
+    mival(count,1) = abs((m_raw(count)));
+end
+
+if num_trials > 1
+    mival = mean(mival);
+end
+

code/LFP/PAC/subfunctions/morlet_wav.m

@@ -0,0 +1,25 @@
+function y = morlet_wav(f,t,width)
+% function y = morlet_wav(f,t,width)
+% 
+% Create a Morlet wavelet 'y' with frequency resolution 'f' and temporal 
+% resolution 't'. The wavelet will be normalized so the total energy is 1.
+% The 'width' defines the temporal and frequency resolution for the given
+% centre frequency 'f' by determining the number of cycles of the wavelet
+% itself (see Tallon-Baudry et al., J. Neurosci. 15, 722-734 (1997) or 
+% Event-Related Potentials: A Methods Handbook, Handy (editor), MIT Press, 
+% (2005))
+%
+% INPUTS:
+% f - frequency to filter at 
+% s - signal to be filtered
+% Fs - sampling frequency of s
+% width - parameter which defines the mother wavelet (>= 5 is suggested)
+%
+% Author: Ole Jensen, August 1998 
+
+
+sf = f/width;
+st = 1/(2*pi*sf);
+A = 1/sqrt((st*sqrt(pi)));
+y = A*exp(-t.^2/(2*st^2)).*exp(1i*2*pi*f.*t);
+

code/LFP/PAC/subfunctions/my_sig_test.m

@@ -0,0 +1,20 @@
+function [h,p] = my_sig_test(pacval, pac_shf_values, alpha)
+% function [h,p] = my_sig_test(pacval, pac_shf_values, alpha)
+%
+% Returns a binary value for significance 'h' and a p-value 'p', given a
+% value under analysis 'pacval', a distribution of values for comparison
+% 'pac_shf_values' and a significance level 'alpha'
+%
+% INPUTS:
+% pacval - PAC value 
+% pac_shf_values - distribution of PAC values obtained from shuffled data
+% sets 
+% alpha - significance level, i.e. if alpha = 0.05 then 'pacval' must fall
+% in the top 5% of the 'pac_shf_values' to be deemed significant
+%
+% Author: Angela Onslow, May 2010
+  
+
+p = mean (pac_shf_values > pacval);
+
+h = (p <= alpha);

code/LFP/PAC/subfunctions/my_sig_test_fdr.m

@@ -0,0 +1,48 @@
+function [h,p] = my_sig_test_fdr(pacval, pac_shf_values, alpha, measure)
+% function [h,p] = my_sig_test_fdr(pacval, pac_shf_values, alpha, measure)
+%
+% Returns a binary value for significance 'h' and a p-value 'p', given a
+% value under analysis 'pacval', a distribution of values for comparison
+% 'pac_shf_values' and a significance level 'alpha'
+%
+% INPUTS:
+% pacval - PAC value 
+% pac_shf_values - distribution of PAC values obtained from shuffled data
+% sets 
+% alpha - significance level, i.e. if alpha = 0.05 then 'pacval' must fall
+% in the top 5% of the 'pac_shf_values' to be deemed significant
+%
+% Author: Angela Onslow, May 2010
+
+if strcmp(measure, 'esc')
+    
+    p_tail1 = mean (pac_shf_values > pacval);
+    p_tail2 = mean (pac_shf_values < pacval);
+    
+    if (p_tail1 <= alpha/2)
+        h = 1;
+        p = 2*(p_tail1);
+    elseif (p_tail2 <= alpha/2)
+        h = 1;
+        p = 2*(p_tail2);
+    elseif (p_tail1 > alpha/2)
+        h = 0;
+        p = 2*(p_tail1);
+    elseif (p_tail2 > alpha/2)
+        h = 0;
+        p= 2*(p_tail2);
+    end
+ 
+    
+else
+    
+    % percentage of times pacval is greater than a pac_shf_value
+    p = mean (pac_shf_values > pacval);
+    
+    h = (p <= alpha);
+    
+end
+
+if p == 0
+    p =  (1/length(pac_shf_values))-(1/(10*length(pac_shf_values)));
+end

code/LFP/PAC/subfunctions/phasevec.m

@@ -0,0 +1,30 @@
+function y = phasevec(f,s,Fs,width)
+% function y = phasevec(f,s,Fs,width)
+%
+% Returns a vector 'y' containing the instantaneous phase values of signal 
+% 's' filtered for frequency 'f' (via convolution with a complex Morlet 
+% wavelet described by 'width')
+%
+% INPUTS:
+% f - frequency to filter at 
+% s - signal to be filtered
+% Fs - sampling frequency of s
+% width - parameter which defines the mother wavelet (this is then 
+% scaled and translated in order to filter for different frequencies, 
+% >= 5 is suggested, see Tallon-Baudry et al., J. Neurosci. 15, 722-734 
+% (1997))
+%
+% NOTE: This function is a modification of code written Ole Jensen, August
+% 1998, see ENERGYVEC
+%
+% Author: Angela Onslow, May 2010
+
+dt = 1/Fs;
+sf = f/width;
+st = 1/(2*pi*sf); 
+
+t=-3.5*st:dt:3.5*st; 
+m = morlet_wav(f,t,width);
+y = conv(s,m');
+y = angle(y);
+y = y(ceil(length(m)/2):length(y)-floor(length(m)/2));

code/LFP/PAC/subfunctions/relat_shaf.m

@@ -0,0 +1,24 @@
+function PAC0=relat_shaf(PAC)
+% 将原始数据PAC值与shaffle数据PAC做相对值。PAC为结构数据
+[m,n]=size(PAC.pacmat);
+relat=PAC.pacmat-PAC.shf_data_mean;
+relatperc=(PAC.pacmat-PAC.shf_data_mean)./PAC.shf_data_mean;
+relat_mi=(PAC.pacmat-PAC.shf_data_mean)./PAC.shf_data_std;
+for i=1:n
+    for j=1:m
+        if relat(j,i)<0
+            relat(j,i)=0;
+        end
+        if relatperc(j,i)<0
+            relatperc(j,i)=0;
+        end
+        if relat_mi(j,i)<0
+            relat_mi(j,i)=0;
+        end
+    end
+end
+PAC.pac_relat=relat;
+PAC.pac_relatperc=relatperc;
+PAC.relat_mi=relat_mi;
+PAC0=PAC;
+end

code/LFP/PAC/subfunctions/shuffle_esc.m

@@ -0,0 +1,96 @@
+function shf_sig = shuffle_esc(sig, Fs)
+% function shf_sig = shuffle_esc(sig, Fs)
+%
+% This function shuffles a signal using random insertion type shuffling. 
+% The signal 'sig' is divided into a number of sections, equal to the
+% number of seconds of data or 1000, which ever is smaller. These sections 
+% are of random lengths, chosen from a uniform distribution. The sections 
+% are then randomly re-ordered. 
+% The input signal 'sig' must contain several seconds on data (>2) and must
+% be a vector (i.e. number of trials = 1) or a cell array containing 
+% multiple signals
+%
+% INPUTS:
+% sig - input signal which is to be shuffled, passed as a column vector
+% Fs - sampling frequency of 'sig', in Hz
+%
+% OUTPUTS:
+% shf_sig - either a vector or a cell array depending on the input 'sig'.
+% When used within find_pac_shf.m this output is of the same dimension as 
+% filt_sig_mod and filt_sig_pac: number of cells - number of frequency bins
+% and each cell element is a matrix(num_data_points, num_trials)
+%
+% Author: Rafal Bogacz, Angela Onslow, May 2010
+
+sig_type = class(sig);
+
+if iscell(sig)
+    ybins = size(sig,1);
+    xbins = size(sig,2);
+    shf_sig_cell = cell(size(sig));
+    num_sec = ceil((size(sig{1,1},1))/Fs);
+else
+    shf_sig = [];
+    num_sec = ceil((size(sig,1))/Fs);
+end
+
+if num_sec > 1000
+    num_sec = 1000;
+end
+
+
+switch sig_type
+    
+    case 'double'
+        
+        % Choose num_sec random 'cut' positions
+        dpsplit = ceil(size(sig,1).*rand(num_sec,1));
+        % Arrange these in ascending order
+        dpsplit = sort (dpsplit);
+
+        start(1)=1;
+        start(2:num_sec)=dpsplit(1:num_sec-1);
+        ending(1:num_sec-1)=dpsplit(1:num_sec-1)-1;
+        ending(num_sec) =  size(sig,1);
+
+        order = randperm(num_sec);  % Function randperm: Random permutation
+    
+        for c = 1:num_sec
+        
+            %shuffle the signal
+            shf_sig = [shf_sig; sig(start(order(c)):ending(order(c)),:)];
+        
+        end
+
+
+    case 'cell'
+        for i = 1:ybins
+            for j = 1:xbins
+                
+                current_sig = sig{i,j};
+                
+                % Choose num_sec random 'cut' positions
+                dpsplit = ceil(size(sig{1,1},1).*rand(num_sec,1));
+                % Arrange these in ascending order
+                dpsplit = sort (dpsplit);
+
+                start(1)=1;
+                start(2:num_sec)=dpsplit(1:num_sec-1);
+                ending(1:num_sec-1)=dpsplit(1:num_sec-1)-1;
+                ending(num_sec) =  size(sig{1,1},1);
+
+                order = randperm(num_sec);
+    
+                for c = 1:num_sec
+        
+                    %shuffle the signal
+                    shf_sig_cell{i,j} = [shf_sig_cell{i,j}; current_sig(start(order(c)):ending(order(c)),:)];
+        
+                end
+                
+            end
+        end
+        
+        shf_sig = shf_sig_cell;
+end
+