update to phase-amplitude coupling. time analysis should now work. exports aveStats. UI reworked to work better with low resolution displays.

Author: caccavanoap

analyzeLFPBatch.m

@@ -60,6 +60,7 @@ function analyzeLFPBatch(param, dataFolder, saveFolder, expEvFolder, expDataFold
 %     param.xFreqBin         = Frequency bin size for n x n PAC analysis (Default = 5 Hz)
 %     param.xFreqLow         = cell: low frequency band for x-freq (Theta, Alpha, Beta, SW)
 %     param.morlWidth        = width/number of cycles of the morlet wavelet filter, default = 7
+%     param.nShuffle         = # shuffles to calculate Z-value for total PAC - does not due for nxn or time PAC, very computationally expensive. (default = 200)
 %     param.winLength        = time binning for phase-amplitude analysis (s). Dictates min low freq (=1/winLength), so default = 0.5s results in min freq. of 2Hz
 %     param.winOverlap       = Amount to overlap time bins (default = 0.2s)
 %     param.importStimOption = option to import stim file from pClamp (default = 0)
@@ -140,6 +141,7 @@ function analyzeLFPBatch(param, dataFolder, saveFolder, expEvFolder, expDataFold
 if ~isfield(param,'xFreqBin');         param.xFreqBin          = 5;    end  % [Hz]
 if ~isfield(param,'xFreqLow');         param.xFreqLow          = 'Theta'; end
 if ~isfield(param,'morlWidth');        param.morlWidth         = 7;    end
+if ~isfield(param,'nShuffle');         param.nShuffle          = 200;  end 
 if ~isfield(param,'winLength');        param.winLength         = 0.5;  end
 if ~isfield(param,'winOverlap');       param.winOverlap        = 0.2;  end
 if ~isfield(param,'importStimOption'); param.importStimOption  = 0;    end

analyzeLFPFile.m

@@ -67,6 +67,7 @@
 %     param.xFreqOption      = boolean flag to perform cross-frequency analysis
 %     param.xFreqBin         = Frequency bin size for n x n PAC analysis (Default = 5 Hz)
 %     param.xFreqLow         = cell: low frequency band for x-freq (Theta, Alpha, Beta, SW)
+%     param.nShuffle         = # shuffles to calculate Z-value for total PAC - does not due for nxn or time PAC, very computationally expensive. (default = 200)
 %     param.morlWidth        = width/number of cycles of the morlet wavelet filter, default = 7
 %     param.winLength        = time binning for phase-amplitude analysis (s). Dictates min low freq (=1/winLength), so default = 0.5s results in min freq. of 2Hz
 %     param.winOverlap       = Amount to overlap time bins (default = 0.2s)
@@ -157,6 +158,7 @@
 if ~isfield(param,'xFreqBin');         param.xFreqBin          = 5;    end  % [Hz]
 if ~isfield(param,'xFreqLow');         param.xFreqLow          = 'Theta'; end
 if ~isfield(param,'morlWidth');        param.morlWidth         = 7;    end
+if ~isfield(param,'nShuffle');         param.nShuffle          = 200;  end 
 if ~isfield(param,'winLength');        param.winLength         = 0.5;  end  % [s]
 if ~isfield(param,'winOverlap');       param.winOverlap        = 0.2;  end  % [s]
 if ~isfield(param,'importStimOption'); param.importStimOption  = 0;    end
@@ -853,18 +855,17 @@
 
 
 %% Cross-Frequency Phase-Amplitude Coupling (PAC)
+%  Adapted from Canolty et al 2006 and Onslow et al 2011
 if param.xFreqOption
-  % Below is an adaptation/simplification of code written by Author: Angela Onslow, May 2010
-  fprintf(['calculating phase-amplitude coupling for file ' dataFileName '... ']);
-  
-  %% Total PAC Analysis for n x n matrix:
   % Assign a temp LFP tSeries (it will be trimmed later)
   tSeries  = data.LFP.tSeries;
   nSample  = length(data.LFP.tSeries);
-  nShuffle = 200; % Add to UI?
 
+  %% Total PAC Analysis for n x n matrix:
   if isfield(param, 'spectLim1') && isfield(param, 'spectLim2')
 
+    fprintf(['calculating n x n phase-amplitude coupling (file ' dataFileName ')... ']);
+    
     % Initialize data structures:
     if ~isfield(data.LFP,'xFreq'); data.LFP.xFreq = struct; end
 
@@ -876,29 +877,32 @@
     nFreq = length(data.LFP.xFreq.morlFreq);
 
     % Initialize phase, amplitude, and PAC arrays
-    data.LFP.xFreq.modPhase = zeros(nSample, nFreq);
-    data.LFP.xFreq.pacAmp   = zeros(nSample, nFreq);
-    data.LFP.xFreq.pacMI    = zeros(nFreq, nFreq);
+    data.LFP.xFreq.phsPAC = zeros(nSample, nFreq);
+    data.LFP.xFreq.ampPAC = zeros(nSample, nFreq);
+    data.LFP.xFreq.pacMI  = zeros(nFreq, nFreq);
 
     % Determine phase and amplitude via Morlet wavelet:
     for i = 1 : nFreq
-      data.LFP.xFreq.modPhase(:, i) = morletPhase(data.LFP.xFreq.morlFreq(i), tSeries, param.Fs, param.morlWidth);
-      data.LFP.xFreq.pacAmp(:, i)   = morletAmp(data.LFP.xFreq.morlFreq(i), tSeries, param.Fs, param.morlWidth);
+      data.LFP.xFreq.phsPAC(:,i) = morletPhase(data.LFP.xFreq.morlFreq(i), tSeries, param.Fs, param.morlWidth);
+      data.LFP.xFreq.ampPAC(:,i) = morletAmp(data.LFP.xFreq.morlFreq(i), tSeries, param.Fs, param.morlWidth);
     end
 
     % Calculate PAC modulation index (MI):
-    for i = 1 : nFreq
-      for j = 1 : nFreq
-        data.LFP.xFreq = calcPACMI(data.LFP.xFreq.modPhase(:,j), data.LFP.xFreq.pacAmp(:,i), nSample, param.Fs, nShuffle);
-      end
-    end
+    data.LFP.xFreq = calcPACMI(data.LFP.xFreq, param.Fs, 0); % nShuffle > 1 gets very computationally expensive!
+    
+    % Order structure:
+    data.LFP.xFreq = orderStruct(data.LFP.xFreq);
+    
+    fprintf('done\n');
   end
 
   %% Total PAC Analysis for each higher frequency band:
   % Calculate modulating phase (lower frequency) via Morlet wavelet:
   if isfield(data, param.xFreqLow)
+    fprintf(['calculating Z-corrected phase-amplitude coupling for frequency bands of interest (file ' dataFileName ')... ']);
+    
     morlFreqP = data.(param.xFreqLow).lim1 + floor((data.(param.xFreqLow).lim2 - data.(param.xFreqLow).lim1)/2);
-    modPhase  = morletPhase(morlFreqP, tSeries, param.Fs, param.morlWidth);
+    phsPAC    = morletPhase(morlFreqP, tSeries, param.Fs, param.morlWidth);
 
     % Determine available higher phase-modulated frequency(ies) available:
     nHi = 0;
@@ -920,63 +924,71 @@
     end
 
     if nHi > 0
-      
-      % Initialize higher frequency arrays:
-      morlFreqA = zeros(1, nHi);
-      pacAmp    = zeros(nSample, nHi);
-      
       for i = 1:nHi
-        
         if ~isfield(data.(xFreqHi{i}),'xFreq'); data.(xFreqHi{i}).xFreq = struct; end
 
         % Calculate amplitude of higher frequency(ies) via Morlet wavelet:
-        morlFreqA(i) = data.(xFreqHi{i}).lim1 + floor((data.(xFreqHi{i}).lim2 - data.(xFreqHi{i}).lim1)/2);
-        pacAmp(:,i)  = morletAmp(morlFreqA(i), tSeries, param.Fs, param.morlWidth);
-        
-        % Calculate total PAC measure:
-        data.(xFreqHi{i}).xFreq = calcPACMI(modPhase, pacAmp(:,i), nSample, Fs, nShuffle);
+        morlFreqA = data.(xFreqHi{i}).lim1 + floor((data.(xFreqHi{i}).lim2 - data.(xFreqHi{i}).lim1)/2);
+        ampPAC    = morletAmp(morlFreqA, tSeries, param.Fs, param.morlWidth);
 
         % Update data structure:
         data.(xFreqHi{i}).xFreq.xFreqLow  = param.xFreqLow;
         data.(xFreqHi{i}).xFreq.xFreqHi   = xFreqHi{i};
         data.(xFreqHi{i}).xFreq.morlFreqP = morlFreqP;
-        data.(xFreqHi{i}).xFreq.morlFreqA = morlFreqA(i);
-        data.(xFreqHi{i}).xFreq.modPhase  = modPhase;
-        data.(xFreqHi{i}).xFreq.pacAmp    = pacAmp(:,i);
+        data.(xFreqHi{i}).xFreq.morlFreqA = morlFreqA;
+        data.(xFreqHi{i}).xFreq.phsPAC    = phsPAC;
+        data.(xFreqHi{i}).xFreq.ampPAC    = ampPAC;
+
+        % Calculate total PAC measure:
+        data.(xFreqHi{i}).xFreq = calcPACMI(data.(xFreqHi{i}).xFreq, param.Fs, param.nShuffle);
       end
 
       %% Time PAC Analysis
       % Truncate signals to get integer number of time windows
       nSampWn   = ceil(param.winLength * param.Fs);
       nSampOl   = ceil(param.winOverlap * param.Fs);
       remSample = mod(nSample, nSampWn);
-      modPhase  = modPhase(1 : nSample - remSample);
-      pacAmp    = pacAmp(1 : nSample - remSample, :);
+      phsPAC    = phsPAC(1 : nSample - remSample);
+      ampPAC    = ampPAC(1 : nSample - remSample, :);
 
       % Update nSample
-      nSample = length(modPhase);
+      nSample = length(phsPAC);
       idx     = bsxfun(@plus, (1:nSampWn)', 1+(0:(fix((nSample - nSampOl)/(nSampWn - nSampOl)) - 1))*(nSampWn - nSampOl)) - 1;
       nWin    = size(idx,2);
 
       % Determine average time of windows:
       timingWin = zeros(nWin, 1);
-      for j = 1:size(idx, 2)
-        timingWin(j) = mean(data.LFP.timing(idx(:, j)));
+      for k = 1:size(idx, 2)
+        timingWin(k) = mean(data.LFP.timing(idx(:, k)));
       end
-
+      
       % Calculate windowed time-series PAC:
-      pacMIWin  = zeros(nWin, 1);
+      pacMIWin       = zeros(nWin, 1);
+      pacMIWin_Len   = zeros(nWin, 1);
+      pacMIWin_Phase = zeros(nWin, 1);
+      
       for i = 1:nHi
-        for j = 1:size(idx, 2)
+        for k = 1:size(idx, 2)
+          
+          z = ampPAC(idx(:,k), i) .* exp(1i * phsPAC(idx(:,k))); % Create composite signal
+          pacMIWin(k)       = mean(z);  % Compute the mean length of composite signal
+          pacMIWin_Len(k)   = abs(pacMIWin(k));
+          pacMIWin_Phase(k) = angle(pacMIWin(k));
 
-          z = pacAmp(idx(:,j), i) .* exp(1i * modPhase(idx(:,j))); % Create composite signal
-          pacRaw = mean(z);  % Compute the mean length of composite signal
-          pacMIWin(j) = abs(pacRaw);
-
-          % Update data structure:
-          data.(xFreqHi{i}).xFreq.pacMIWin  = pacMIWin;
-          data.(xFreqHi{i}).xFreq.timingWin = timingWin;
         end
+        
+        % Update data structure:
+        data.(xFreqHi{i}).xFreq.pacMIWin        = pacMIWin;
+        data.(xFreqHi{i}).xFreq.pacMIWin_Len    = pacMIWin_Len;
+        data.(xFreqHi{i}).xFreq.pacMIWin_Phase  = pacMIWin_Phase;
+        data.(xFreqHi{i}).xFreq.timingWin       = timingWin;
+        
+        % Linear interpolation of pacMIWin_Len to data samplingInt for later correlations
+        data.(xFreqHi{i}).xFreq.timingI = data.LFP.timing(1:nSample);
+        data.(xFreqHi{i}).xFreq.pacMIWin_LenI = interp1(timingWin, pacMIWin_Len, data.(xFreqHi{i}).xFreq.timingI);
+        
+        % Order structure:
+        data.(xFreqHi{i}).xFreq = orderStruct(data.(xFreqHi{i}).xFreq);
       end
     end
   end

calcPACMI.m

@@ -1,24 +1,46 @@
-function xFreq = calcPACMI(modPhase, pacAmp, nSample, Fs, nShuffle)
+function xFreq = calcPACMI(xFreq, Fs, nShuffle)
 % Computes normalized phase-amplitude coupled modulation index
-% Adapted from Canolty et al Science 2006
-% See also Onslow et al 2011
+% Adapted from Canolty et al 2006 and Onslow et al 2011
 
-z = pacAmp .* exp(1i * modPhase); % Create composite signal
-pacRaw = mean(z);  % Compute the mean length of composite signal
-xFreq.pacMI_Raw = abs(pacRaw);
-pacShuffle = zeros(nShuffle, 1);
+nSample = size(xFreq.phsPAC, 1);
+nPhs    = size(xFreq.phsPAC, 2);
+nAmp    = size(xFreq.ampPAC, 2);
 
-% Compute shuffled values:
-for i = 1 : nShuffle
-  skip = shuffleInd(nSample, Fs, nShuffle);
-  pacAmpShuffle = [pacAmp(skip(i):end) pacAmp(1:skip(s)-1)];
-  pacShuffle(i) = abs(mean(pacAmpShuffle .* exp(1i * modPhase)));
-end
+% Initialize output arrays:
+xFreq.pacMI       = zeros(nAmp, nPhs);
+xFreq.pacMI_Len   = zeros(nAmp, nPhs);
+xFreq.pacMI_Phase = zeros(nAmp, nPhs);
 
-% Fit gaussian to surrogate data, uses normfit.m from MATLAB Statistics toolbox 
-[xFreq.shuffleMean, xFreq.shuffleSTD] = normfit(pacShuffle); 
+if nShuffle > 1
+  xFreq.pacMIShf_LenAve = zeros(nAmp, nPhs);
+  xFreq.pacMIShf_LenSTD = zeros(nAmp, nPhs);
+  xFreq.pacMI_Z         = zeros(nAmp, nPhs);
+  xFreq.pacMI_LenZ      = zeros(nAmp, nPhs);
+end
 
-% Normalize length using shuffled data (z-score) 
-xFreq.pacMI_Length = (xFreq.pacMI_Raw - xFreq.shuffleMean) / xFreq.shuffleSTD; 
-xFreq.pacMI_Phase = angle(pacRaw);
-xFreq.pacMI_Norm = pacMI_Length * exp(1i * pacMI_Phase);
+for i = 1 : nAmp
+  for j = 1 : nPhs
+    z = xFreq.ampPAC(:,i) .* exp(1i * xFreq.phsPAC(:,j)); % Create composite signal
+    xFreq.pacMI(i,j)       = mean(z);  % Compute the mean length of composite signal
+    xFreq.pacMI_Len(i,j)   = abs(xFreq.pacMI(i,j));
+    xFreq.pacMI_Phase(i,j) = angle(xFreq.pacMI(i,j));
+    
+    if nShuffle > 1
+      pacMIShf_Len = zeros(nShuffle, 1);
+      
+      % Compute shuffled values:
+      for k = 1 : nShuffle
+        skip = shuffleInd(nSample, Fs, nShuffle);
+        ampPacShf = vertcat(xFreq.ampPAC(skip(k):end,j), xFreq.ampPAC(1:skip(k)-1,j));
+        pacMIShf_Len(k) = abs(mean(ampPacShf .* exp(1i * xFreq.phsPAC(:,i))));
+      end
+      
+      % Fit gaussian to shuffled data, uses normfit.m from MATLAB Statistics toolbox
+      [xFreq.pacMIShf_LenAve(i,j), xFreq.pacMIShf_LenSTD(i,j)] = normfit(pacMIShf_Len);
+      
+      % Normalize length using shuffled data (z-score)
+      xFreq.pacMI_LenZ(i,j) = (xFreq.pacMI_Len(i,j) - xFreq.pacMIShf_LenAve(i,j)) / xFreq.pacMIShf_LenSTD(i,j);
+      xFreq.pacMI_Z(i,j) = xFreq.pacMI_LenZ(i,j) * exp(1i * xFreq.pacMI_Phase(i,j));
+    end
+  end
+end

exportAveStats.m

@@ -78,15 +78,22 @@ function exportAveStats(data, saveFile, exportFile)
 if isfield(data, 'gamma')
   % Total Stats:
   if isfield(data.gamma, 'tPower'); outTable = [outTable table(10^6 * data.gamma.tPower, 'VariableNames', {'Gamma_Tot_Power_uV2'})]; end
-  if isfield(data.gamma, 'FFT') 
+  if isfield(data.gamma, 'FFT')
     if isfield(data.gamma.FFT, 'pkFreq');  outTable = [outTable table(data.gamma.FFT.pkFreq, 'VariableNames', {'Gamma_FFT_pkFreq_Hz'})]; end
-  	if isfield(data.gamma.FFT, 'fitMean'); outTable = [outTable table(data.gamma.FFT.fitMean, 'VariableNames', {'Gamma_FFT_fitMean_Hz'})]; end
+    if isfield(data.gamma.FFT, 'fitMean'); outTable = [outTable table(data.gamma.FFT.fitMean, 'VariableNames', {'Gamma_FFT_fitMean_Hz'})]; end
     if isfield(data.gamma.FFT, 'fitSD');   outTable = [outTable table(data.gamma.FFT.fitSD, 'VariableNames', {'Gamma_FFT_fitSD_Hz'})]; end
-  	if isfield(data.gamma.FFT, 'fitFWHM'); outTable = [outTable table(data.gamma.FFT.fitFWHM, 'VariableNames', {'Gamma_FFT_fitFWHM_Hz'})]; end
+    if isfield(data.gamma.FFT, 'fitFWHM'); outTable = [outTable table(data.gamma.FFT.fitFWHM, 'VariableNames', {'Gamma_FFT_fitFWHM_Hz'})]; end
   end
-  if isfield(data.gamma, 'phase') 
+  if isfield(data.gamma, 'phase')
     if isfield(data.gamma.phase, 'nCycle'); outTable = [outTable table(data.gamma.phase.nCycle, 'VariableNames', {'Gamma_Tot_nCycle'})]; end
-  	if isfield(data.gamma.phase, 'phFreq'); outTable = [outTable table(data.gamma.phase.phFreq, 'VariableNames', {'Gamma_Tot_phFreq_Hz'})]; end
+    if isfield(data.gamma.phase, 'phFreq'); outTable = [outTable table(data.gamma.phase.phFreq, 'VariableNames', {'Gamma_Tot_phFreq_Hz'})]; end
+  end
+  if isfield(data.gamma, 'xFreq')
+    if isfield(data.gamma.xFreq, 'pacMI_Len');       outTable = [outTable table(data.gamma.xFreq.pacMI_Len, 'VariableNames', {'Gamma_pacMI_Len'})]; end
+    if isfield(data.gamma.xFreq, 'pacMI_Phase');     outTable = [outTable table(data.gamma.xFreq.pacMI_Phase, 'VariableNames', {'Gamma_pacMI_Phase'})]; end
+    if isfield(data.gamma.xFreq, 'pacMI_LenZ');      outTable = [outTable table(data.gamma.xFreq.pacMI_LenZ, 'VariableNames', {'Gamma_pacMI_LenZ'})]; end
+    if isfield(data.gamma.xFreq, 'pacMIShf_LenAve'); outTable = [outTable table(data.gamma.xFreq.pacMIShf_LenAve, 'VariableNames', {'Gamma_pacMIShf_LenAve'})]; end
+    if isfield(data.gamma.xFreq, 'pacMIShf_LenSTD'); outTable = [outTable table(data.gamma.xFreq.pacMIShf_LenSTD, 'VariableNames', {'Gamma_pacMIShf_LenSTD'})]; end
   end
   % SWR Stats:
   if isfield(data.gamma, 'SWR')
@@ -115,6 +122,13 @@ function exportAveStats(data, saveFile, exportFile)
     if isfield(data.hgamma.phase, 'nCycle'); outTable = [outTable table(data.hgamma.phase.nCycle, 'VariableNames', {'HGamma_Tot_nCycle'})]; end
   	if isfield(data.hgamma.phase, 'phFreq'); outTable = [outTable table(data.hgamma.phase.phFreq, 'VariableNames', {'HGamma_Tot_phFreq_Hz'})]; end
   end
+  if isfield(data.hgamma, 'xFreq')
+    if isfield(data.hgamma.xFreq, 'pacMI_Len');       outTable = [outTable table(data.hgamma.xFreq.pacMI_Len, 'VariableNames', {'HGamma_pacMI_Len'})]; end
+    if isfield(data.hgamma.xFreq, 'pacMI_Phase');     outTable = [outTable table(data.hgamma.xFreq.pacMI_Phase, 'VariableNames', {'HGamma_pacMI_Phase'})]; end
+    if isfield(data.hgamma.xFreq, 'pacMI_LenZ');      outTable = [outTable table(data.hgamma.xFreq.pacMI_LenZ, 'VariableNames', {'HGamma_pacMI_LenZ'})]; end
+    if isfield(data.hgamma.xFreq, 'pacMIShf_LenAve'); outTable = [outTable table(data.hgamma.xFreq.pacMIShf_LenAve, 'VariableNames', {'HGamma_pacMIShf_LenAve'})]; end
+    if isfield(data.hgamma.xFreq, 'pacMIShf_LenSTD'); outTable = [outTable table(data.hgamma.xFreq.pacMIShf_LenSTD, 'VariableNames', {'HGamma_pacMIShf_LenSTD'})]; end
+  end
   % SWR Stats:
   if isfield(data.hgamma, 'SWR')
     if isfield(data.hgamma.SWR, 'power'); outTable = [outTable table(10^6 * mean(data.hgamma.SWR.power,'omitnan'), 'VariableNames', {'HGamma_SWR_Power_uV2'})]; end

morletAmp.m

@@ -33,6 +33,6 @@
 t  = -(width/2)*st : dt : (width/2)*st; 
 A  = 1 / sqrt((st*sqrt(pi)));
 m  = A * exp(-t.^2/(2*st^2)) .* exp(1i*2*pi*f.*t);
-y = conv(s, m');
-y = abs(y);
-y = y(ceil(length(m)/2) : length(y) - floor(length(m)/2));
+y  = conv(s, m');
+y  = abs(y);
+y  = y(ceil(length(m)/2) : length(y) - floor(length(m)/2));

morletPhase.m

@@ -33,6 +33,6 @@
 t  = -(width/2)*st : dt : (width/2)*st; 
 A  = 1 / sqrt((st*sqrt(pi)));
 m  = A * exp(-t.^2/(2*st^2)) .* exp(1i*2*pi*f.*t);
-y = conv(s, m');
-y = angle(y);
-y = y(ceil(length(m)/2) : length(y) - floor(length(m)/2));
+y  = conv(s, m');
+y  = angle(y);
+y  = y(ceil(length(m)/2) : length(y) - floor(length(m)/2));