Skip to content

Commit

Permalink
fix formatting
Browse files Browse the repository at this point in the history
  • Loading branch information
@ryanharvey1
ryanharvey1 committed Nov 20, 2024
1 parent fa32263 commit b0cc816
Showing 1 changed file with 95 additions and 95 deletions.
190 changes: 95 additions & 95 deletions lfp/getEMGFromLFP.m
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,15 @@
function [EMGFromLFP] = getEMGFromLFP(basepath,varargin)
function [EMGFromLFP] = getEMGFromLFP(basepath, varargin)
% Based on Erik Schomburg's work and code. Grabs channels and calculates
% their correlations in the 300-600Hz band over sliding windows of 0.5sec.

% Channels are automatically selected and are a combination of first and last channels
% on each shank. This is based on the xml formatting standard that channel ordering goes
% from superficial to deep for each channel/spike group.
% on each shank. This is based on the xml formatting standard that channel ordering goes
% from superficial to deep for each channel/spike group.
%
% Special channels should be 0-indexed, per neuroscope convention
% Requires .lfp/lfp and .xml. Assumes each spikegroup in the .xml
% represents a "shank"
%
%
% Mean pairwise correlations are calculated for each time point.

% =========================================================================
Expand All @@ -34,10 +34,10 @@
% 'samplingFrequency' - desired sampling rate for EMG output. default:2
% 'noPrompts' (default: false) prevents prompts about saving/adding metadata
% 'fromDat' -uses the .dat file instead of .lfp (default:false)
%
%
%
% OUTPUTS
%
%
% EMGFromLFP - struct of the LFP datatype. saved at
% basePath/baseName.EMGFromLFP.LFP.mat
% .timestamps - timestamps (in seconds) that match .data samples
Expand All @@ -53,23 +53,23 @@

%% Buzcode name of the EMGCorr.LFP.mat file
recordingname = basenameFromBasepath(basepath);
matfilename = fullfile(basepath,[recordingname,'.EMGFromLFP.LFP.mat']);
matfilename = fullfile(basepath, [recordingname, '.EMGFromLFP.LFP.mat']);

%% xmlPameters
p = inputParser;
addParameter(p,'restrict',[0 inf],@isnumeric)
addParameter(p,'specialChannels',[],@isnumeric)
addParameter(p,'rejectChannels',[],@isnumeric)
addParameter(p,'restrictChannels',[],@isnumeric)
addParameter(p,'saveMat',true,@islogical)
addParameter(p,'saveLocation','',@isstr)
addParameter(p,'overwrite',false,@islogical)
addParameter(p,'samplingFrequency',2,@isnumeric)
addParameter(p,'noPrompts',false,@islogical);
addParameter(p,'fromDat',false,@islogical);
addParameter(p,'chInfo',[],@isstruct);
parse(p,varargin{:})
addParameter(p, 'restrict', [0, inf], @isnumeric)
addParameter(p, 'specialChannels', [], @isnumeric)
addParameter(p, 'rejectChannels', [], @isnumeric)
addParameter(p, 'restrictChannels', [], @isnumeric)
addParameter(p, 'saveMat', true, @islogical)
addParameter(p, 'saveLocation', '', @isstr)
addParameter(p, 'overwrite', false, @islogical)
addParameter(p, 'samplingFrequency', 2, @isnumeric)
addParameter(p, 'noPrompts', false, @islogical);
addParameter(p, 'fromDat', false, @islogical);
addParameter(p, 'chInfo', [], @isstruct);
parse(p, varargin{:})

restrict = p.Results.restrict;
specialChannels = p.Results.specialChannels;
rejectChannels = p.Results.rejectChannels;
Expand All @@ -82,177 +82,178 @@
chInfo = p.Results.chInfo;

if ~isempty(p.Results.saveLocation)
matfilename = fullfile(p.Results.saveLocation,[recordingname,'.EMGFromLFP.LFP.mat']);
matfilename = fullfile(p.Results.saveLocation, [recordingname, '.EMGFromLFP.LFP.mat']);
end

%% Check if EMGCorr has already been calculated for this recording
%If the EMGCorr file already exists, load and return with EMGCorr in hand
if exist(matfilename,'file') && ~overwrite
if exist(matfilename, 'file') && ~overwrite
display('EMGFromLFP Correlation already calculated - loading from EMGFromLFP.LFP.mat')
load(matfilename)
if exist('EMGCorr','var')%for backcompatability
EMGFromLFP = EMGCorr;
end
if ~exist('EMGFromLFP','var')
display([matfilename,' does not contain a variable called EMGFromLFP'])
if exist('EMGCorr', 'var') %for backcompatability
EMGFromLFP = EMGCorr;
end
if ~exist('EMGFromLFP', 'var')
display([matfilename, ' does not contain a variable called EMGFromLFP'])
end
return
end
display('Calculating EMGFromLFP from High Frequency LFP Correlation')

load(fullfile(basepath,[recordingname,'.session.mat']))
load(fullfile(basepath, [recordingname, '.session.mat']))
nChannels = session.extracellular.nChannels;
SpkGrps = session.extracellular.spikeGroups.channels;
Fs = session.extracellular.srLfp;
if exist(fullfile(basepath,[recordingname,'.lfp']),'file')
lfpFile = checkFile('basepath',basepath,'fileTypes',{'.lfp','.eeg'});
lfpFile = [basepath filesep lfpFile(1).name];
elseif exist(fullfile(basepath,[recordingname,'.eeg']),'file')
lfpFile = [basepath filesep recordingname '.eeg'];
if exist(fullfile(basepath, [recordingname, '.lfp']), 'file')
lfpFile = checkFile('basepath', basepath, 'fileTypes', {'.lfp', '.eeg'});
lfpFile = [basepath, filesep, lfpFile(1).name];
elseif exist(fullfile(basepath, [recordingname, '.eeg']), 'file')
lfpFile = [basepath, filesep, recordingname, '.eeg'];
end
if fromDat
datFile = checkFile('basepath',basepath,'fileTypes',{'.dat'});
datFile = [basepath filesep datFile(1).name];
datFile = checkFile('basepath', basepath, 'fileTypes', {'.dat'});
datFile = [basepath, filesep, datFile(1).name];
datFs = session.extracellular.sr;
end

%%

binScootS = 1 ./ samplingFrequency;
binScootSamps = round(Fs*binScootS); % must be integer, or error on line 190
corrChunkSz = 20; %for batch-processed correlations


%% Pick channels and to analyze
% get spike groups,
% pick every other one... unless specialshanks, in which case pick non-adjacent
%This is potentially dangerous in combination with rejectChannels... i.e.
%what if you pick every other shank but then the ones you pick are all
%reject because noisy shank.

% xcorrs_chs is a list of channels that will be loaded
% spkgrpstouse is a list of spike groups to find channels from
% xcorrs_chs is a list of channels that will be loaded
% spkgrpstouse is a list of spike groups to find channels from

if ~isempty(restrictChannels) % If restrict channel case:
if ~isempty(restrictChannels) % If restrict channel case:
xcorr_chs = restrictChannels;
else
% get list of spike groups (aka shanks) that should be used
usablechannels = [];
spkgrpstouse = [];
if length(SpkGrps)>1, n=1; else n=5; end
if length(SpkGrps) > 1, n = 1;
else n = 5;
end
for gidx = 1:length(SpkGrps)
usableshankchannels{gidx} = setdiff(SpkGrps{gidx},rejectChannels);
usablechannels = cat(2,usablechannels,usableshankchannels{gidx});
usableshankchannels{gidx} = setdiff(SpkGrps{gidx}, rejectChannels);
usablechannels = cat(2, usablechannels, usableshankchannels{gidx});
if ~isempty(usableshankchannels{gidx})
spkgrpstouse = cat(2,spkgrpstouse,gidx);
spkgrpstouse = cat(2, spkgrpstouse, gidx);
end
end

% okay, lets limit usableshankchannels with spkgrpstouse. Thus, empty shanks will be excluded.
% okay, lets limit usableshankchannels with spkgrpstouse. Thus, empty shanks will be excluded.
usableshankchannels = usableshankchannels(spkgrpstouse);

% get list of channels (1 from each good spike group)
xcorr_chs = [];
for gidx=1:length(usableshankchannels)
for gidx = 1:length(usableshankchannels)
%Remove rejectChannels
% usableshankchannels = setdiff(SpkGrps(spkgrpstouse(i)).Channels,rejectChannels);
% usableshankchannels = setdiff(SpkGrps(spkgrpstouse(i)).Channels,rejectChannels);

%grab random channel on each shank
if ~isempty(usableshankchannels{gidx})
randChfromShank = usableshankchannels{gidx}(randi(length(usableshankchannels{gidx}),n,1));
xcorr_chs = [xcorr_chs,randChfromShank];
%grab random channel on each shank
if ~isempty(usableshankchannels{gidx})
randChfromShank = usableshankchannels{gidx}(randi(length(usableshankchannels{gidx}), n, 1));
xcorr_chs = [xcorr_chs, randChfromShank];

end
end
end
xcorr_chs = unique([xcorr_chs,specialChannels]);
xcorr_chs = unique([xcorr_chs, specialChannels]);
end

%% Read and filter channel
% read channels
% old: bz_sessionInfo is 0 indexed (neuroscope) channels,
% but Loadbinary.m takes 1-indexed channel #'s
% new: should use bz_GetLFP
% old: bz_sessionInfo is 0 indexed (neuroscope) channels,
% but Loadbinary.m takes 1-indexed channel #'s
% new: should use bz_GetLFP

switch fromDat
case false
lfp = LoadBinary(lfpFile ,'nChannels',nChannels,'channels',xcorr_chs,...
'start',restrict(1),'duration',diff(restrict),'frequency',Fs); %read and convert to mV
lfp = LoadBinary(lfpFile, 'nChannels', nChannels, 'channels', xcorr_chs, ...
'start', restrict(1), 'duration', diff(restrict), 'frequency', Fs); %read and convert to mV

if isempty(lfp)
error('LFP file empty. Must have valid binary file to compute EMG.')
error('LFP file empty. Must have valid binary file to compute EMG.')
end

case true
lfp = LoadBinary(datFile ,'nChannels',nChannels,'channels',xcorr_chs,...
'start',restrict(1),'duration',diff(restrict),'frequency',datFs,...
'downsample',datFs./Fs); %read and convert to mV
lfp = LoadBinary(datFile, 'nChannels', nChannels, 'channels', xcorr_chs, ...
'start', restrict(1), 'duration', diff(restrict), 'frequency', datFs, ...
'downsample', datFs./Fs); %read and convert to mV
if isempty(lfp)
error('DAT file empty. Must have valid binary file to compute EMG.')
error('DAT file empty. Must have valid binary file to compute EMG.')
end
end

% Filter first in high frequency band to remove low-freq physiologically
% correlated LFPs (e.g., theta, delta, SPWs, etc.)

maxfreqband = floor(max([625 Fs/2]));
maxfreqband = floor(max([625, Fs / 2]));
% xcorr_freqband = [275 300 600 625]; % Hz
xcorr_freqband = [275 300 maxfreqband-25 maxfreqband]; % Hz
xcorr_freqband = [275, 300, maxfreqband - 25, maxfreqband]; % Hz
lfp = filtsig_in(lfp, Fs, xcorr_freqband);

%% xcorr 'strength' is the summed correlation coefficients between channel
% pairs for a sliding window of 25 ms
xcorr_window_samps = round(binScootS*Fs);
xcorr_window_inds = -xcorr_window_samps:xcorr_window_samps;%+- that number of ms in samples
xcorr_window_inds = -xcorr_window_samps:xcorr_window_samps; %+- that number of ms in samples

% new version... batches of correlation calculated at once
timestamps = (1+xcorr_window_inds(end)):binScootSamps:(size(lfp,1)-xcorr_window_inds(end));
timestamps = (1 + xcorr_window_inds(end)):binScootSamps:(size(lfp, 1) - xcorr_window_inds(end));
numbins = length(timestamps);
EMGCorr = zeros(numbins, 1);
% tic
counter = 1;
for j=1:(length(xcorr_chs))
for k=(j+1):length(xcorr_chs)
for j = 1:(length(xcorr_chs))
for k = (j + 1):length(xcorr_chs)
%disp([num2str(counter*2 ./ (length(xcorr_chs)*length(xcorr_chs)*length(timestamps)))])
bz_Counter(counter,(length(xcorr_chs)*(length(xcorr_chs)-1))./2,'Channel Pair')
bz_Counter(counter, (length(xcorr_chs) * (length(xcorr_chs) - 1))./2, 'Channel Pair')
c1 = [];
c2 = [];
binind = 0;
binindstart = 1;
for i = timestamps
binind = binind+1;
s1 =lfp(i + xcorr_window_inds, j);
s2 =lfp(i + xcorr_window_inds, k);
c1 = cat(2,c1,s1);
c2 = cat(2,c2,s2);
if size(c1,2) == corrChunkSz || i == timestamps(end)
binind = binind + 1;
s1 = lfp(i + xcorr_window_inds, j);
s2 = lfp(i + xcorr_window_inds, k);
c1 = cat(2, c1, s1);
c2 = cat(2, c2, s2);
if size(c1, 2) == corrChunkSz || i == timestamps(end)
binindend = binind;
tmp = corr(c1,c2);
tmp = corr(c1, c2);
tmp = diag(tmp);
EMGCorr(binindstart:binindend) = EMGCorr(binindstart:binindend) + tmp;
c1 = [];
c2 = [];
binindstart = binind+1;
end
binindstart = binind + 1;
end
end
counter = counter+1;
counter = counter + 1;
end
end
% toc

EMGCorr = EMGCorr/(length(xcorr_chs)*(length(xcorr_chs)-1)/2); % normalize
EMGCorr = EMGCorr / (length(xcorr_chs) * (length(xcorr_chs) - 1) / 2); % normalize

EMGFromLFP.timestamps = timestamps'./Fs;
EMGFromLFP.timestamps = timestamps' ./ Fs;
EMGFromLFP.data = EMGCorr;
EMGFromLFP.channels = xcorr_chs;
EMGFromLFP.detectorName = 'getEMGFromLFP';
EMGFromLFP.samplingFrequency = samplingFrequency;
EMGFromLFP.samplingFrequency = samplingFrequency;

if ~any(~isnan(EMGFromLFP.data)),
keyboard
end
if saveMat
%Save in buzcodeformat
save(matfilename,'EMGFromLFP');
save(matfilename, 'EMGFromLFP');
end


Expand All @@ -262,7 +263,7 @@
% Created by: Erik Schomburg, 2011

if isnumeric(filtband_or_Filt)
h = fdesign.bandpass(filtband_or_Filt(1), filtband_or_Filt(2), filtband_or_Filt(3), filtband_or_Filt(4), ...
h = fdesign.bandpass(filtband_or_Filt(1), filtband_or_Filt(2), filtband_or_Filt(3), filtband_or_Filt(4), ...
60, 1, 60, Fs);
Filt = design(h, 'butter', 'MatchExactly', 'passband');
else
Expand All @@ -272,17 +273,17 @@
if ~isempty(sig)
if iscell(sig)
filt_sig = cell(size(sig));
for i=1:length(sig(:))
for i = 1:length(sig(:))
filt_sig{i} = filter(Filt, sig{i});
filt_sig{i} = filter(Filt, filt_sig{i}(end:-1:1));
filt_sig{i} = filt_sig{i}(end:-1:1);
end
elseif ((size(sig,1) > 1) && (size(sig,2) > 1))
elseif ((size(sig, 1) > 1) && (size(sig, 2) > 1))
filt_sig = zeros(size(sig));
for i=1:size(filt_sig,2)
filt_sig(:,i) = filter(Filt, sig(:,i));
filt_sig(:,i) = filter(Filt, filt_sig(end:-1:1,i));
filt_sig(:,i) = filt_sig(end:-1:1,i);
for i = 1:size(filt_sig, 2)
filt_sig(:, i) = filter(Filt, sig(:, i));
filt_sig(:, i) = filter(Filt, filt_sig(end:-1:1, i));
filt_sig(:, i) = filt_sig(end:-1:1, i);
end
else
filt_sig = filter(Filt, sig);
Expand All @@ -292,4 +293,3 @@
else
filt_sig = [];
end

0 comments on commit b0cc816

Please sign in to comment.