22-May-2014 12:54:30 function Theta = thetaFeatures(Filebase,Channels, Epochs, varargin) % THETAFEATURES This code takes in a set of channels from an .egg % or .lfp file and outputs the main features of % hippocampal theta (5-12Hz by default). These % include: peak and trough times, phase, and % amplitude. % % 1) uses default parameters % Theta = thetaFeatures( Filebase, Channels, Epochs ) % % 2) supply parameters as a structure with fields an values % Theta = thetaFeatures( ..., options ) % % 3) supply parameters as Property/Value pairs % Theta = thetaFeatures( ..., 'swBP', [2 50], 'WinSize', 300, etc. ) % author: John D. Long II, PhD contact: jlong29@gmail.com % wavelet analysis routines adapted from 'getWavelet' by Eran Stark and % 'wavelet' and 'wave_bases', Torrence and Compo % %%%%%%%%%%%%%% %%% INPUTS %%% %%%%%%%%%%%%%% % Filebase: can be: -full or relative path with or without .xml % extension % e.g. /mnt/Data/buddy140_060813_reo/buddy140_060813_reo.xml % buddy140_060813_reo.xml % buddy140_060813_reo % % Channels: an array of channel ids following LoadBinary format(base 1) % % Epochs: a list of time stamps [start stop] in seconds demarcating % epochs in which to estimate theta features. If this % argument is empty, the entire .eeg/.lfp file will be used. % %%%%%%%%%%%%%%% %%% OUTPUTS %%% %%%%%%%%%%%%%%% % ThetaFeatures: structure of size N, one for each channel with fields % % cell, size M = # of epochs % peaks: per cycle theta peaks (seconds @ phase = 0) % troughs: per cyle theta troughs (seconds @ phase = -/+ pi) % power: z-scored and saved at floor(SR/DSfactor) sampling rate % % The analysis also writes out log file (.txt) containing the code used. % %%%%%%%%%%%%%%%%%%%%%%% %%% FREE PARAMETERS %%% %%%%%%%%%%%%%%%%%%%%%%% % DSfactor: an integer > 1 for downsampling the field potential % thetaBP: [low high], Hz, passband for filtering theta activity % nWaves: the number of wavelets in log spacing within thetaBP % SAVE: boolean flag indicating whether or not to save the output % to file within the Filebase directory % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% SET DEFAULT FREE PARAMETERS %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Defaults.DSfactor = 10; Defaults.thetaBP = [5 12]; Defaults.nWaves = 16; Defaults.SAVE = true; %%%%%%%%%%%%%%%%%%%% %%% INPUT CHECKS %%% %%%%%%%%%%%%%%%%%%%% mfname = mfilename; % 1) fileabase % user might have input an absolute or relative path [pathname, filename, extname] = fileparts(Filebase); if isempty(pathname) pathname = pwd; end % check whether user provided a file extension if isempty(extname) filechk = dir([Filebase '.xml']); if isempty(filechk) error(['%s: incompatible format for 1st argument (Filebase)\n', ... '\tIf user supplied a relative path, the .xml and all\n', ... '\trelevant files must be in the current directory.\n'],mfname); end else filechk = dir(Filebase); if isempty(filechk) error(['%s: incompatible format for 1st argument (Filebase)\n', ... '\tIf user supplied a relative path, the .xml and all\n', ... '\trelevant files must be in the current directory.\n'],mfname); end end % Read in xml of file parameters fprintf(1, 'Loading XML file...\n'); par = LoadXml( Filebase ); % pull parameters from .xml file SR = par.lfpSampleRate; % lfp sampling rate nChan = par.nChannels; % number of channels in the recording % 2) Channels if ~all(Channels > 0 & Channels < nChan) error(['%s: incompatible input for 2nd argument (Channels)\n', ... '\tAccording to the .xml file, the user did not supply\n', ... '\ta valid list of channels.\n'],mfname); end Nchan = length(Channels); % Access lfp looking for an .eeg file first and a .lfp file second if ~isempty(dir([Filebase '.eeg'])) lfp_file = [Filebase '.eeg']; lfp_info = dir(lfp_file); elseif ~isempty(dir([Filebase '.lfp'])) lfp_file = [Filebase '.lfp']; lfp_info = dir(lfp_file); else error(['%s: Field potential file could not be found\n', ... '\tThere is neither a .eeg nor a .lfp file matching the\n', ... '\tbase name of the directory supplied by the user.\n'],mfname); end % 3) Epochs if nargin < 3 || isempty(Epochs) % use the whole recording if no restrictions supplied by user fprintf( 1, '%s: Using Default Epochs:\n', mfname ); fprintf( 1, '\tNo Epochs supplied by user. Using the whole recording.\n'); Epochs = [1/SR lfp_info.bytes/(2*nChan*SR)]; elseif ~all(Epochs(:,1) < Epochs(:,2)) || ... ~all(Epochs(:) > 0 & Epochs(:) < lfp_info.bytes/(2*nChan*SR)) || ... size(Epochs,2) > 2 fprintf( 1, '%s: Epochs are not formatted properly:\n', mfname ); fprintf( 1, '\tEpochs must be formatted as a list of [start stop] times in seconnds.\n'); end Nepochs = size(Epochs,1); %%%%%%%%%%%%%%%%%%%%%%% %%% PARSE ARGUMENTS %%% %%%%%%%%%%%%%%%%%%%%%%% if ~isempty(varargin) % There is no argument validation herein. Defaults = parseArgs(Defaults,varargin); end % return parameters stored in Defaults v2struct(Defaults); % store parameter settings for this run params.filedir = pathname; params.filename = lfp_info.name; params.SReeg = SR; params.DSfactor = DSfactor; params.thetaBP = thetaBP; params.nWaves = nWaves; %%%%%%%%%%%%%%%%%%%%%%% %Access LFP T0 = tic; fprintf(1,'Loading LFP channels...'); lfp = LoadBinary(lfp_file,'frequency',SR,'nChannels',nChan,'channels',Channels); fprintf(1,'in %6.2f seconds\n',toc(T0)); timeline = (1:size(lfp,1))/SR; DSsr = floor(SR/DSfactor); % factor for downsampling % Initialize Output Structure Theta = struct('peaks',cell(size(lfp,2),1),'troughs',cell(size(lfp,2),1), ... 'power',cell(size(lfp,2),1)); % calculate wavelet analysis of each channel for each epoch for cc = 1:Nchan fprintf(1,'Processing Channel %d, %d out of %d...\n',Channels(cc), cc, Nchan); % Initialize per channel variables Peaks = cell(size(Epochs,1),1); Troughs = cell(size(Epochs,1),1); Powers = cell(size(Epochs,1),1); for tt = 1:Nepochs fprintf(1,'Channel %d, Epoch %d out of %d...\n',Channels(cc),tt, Nepochs); startInd = find(timeline >= Epochs(tt,1),1,'first'); endInd = find(timeline <= Epochs(tt,2),1,'last'); segment = lfp(startInd:endInd,cc); timelineEpoch = timeline(startInd:endInd); Nsamples = length(timelineEpoch); % Resample LFP by DSfactor segmentRS = resample(segment,1,DSfactor); timelineRS = (1:DSfactor:Nsamples)/SR + (startInd - 1)/SR; % Run wavelet analysis [w, ~, ~, ~, phases, raw, ~, ~, ~, ~ ] = ... getWavelet(segmentRS,DSsr, thetaBP(1), thetaBP(2), nWaves, 'var',0); % reconstruct signal as the max power across wavelets at each samp [~,maxWInds] = max(w); Nseg = size(raw,2); reconst = zeros(1,Nseg); phase = zeros(1,Nseg); power = zeros(1,Nseg); for ii = 1:Nseg reconst(ii) = raw(maxWInds(ii),ii); phase(ii) = phases(maxWInds(ii),ii); power(ii) = w(maxWInds(ii),ii); end % Now find the zeros (peaks) of theta interpolated up to original sampling % rate % THETA PEAKS: find the zeros of the theta cycle as a 3 point min upon the % rectified signal peaks = zeros(length(phase),1); counter = 0; for ii = 2:length(phase)-1 [~,ind] = min(abs(phase(ii-1:ii+1))); if ind == 2 counter = counter + 1; peaks(counter) = ii; end end peaks = peaks(1:counter); % cycle through all 3 point mins found and interpolate min up to sampling % rate of the original signal ThetaPeaksTs = zeros(counter,1); for ii = 1:counter T0 = phase(peaks(ii)); T1 = phase(peaks(ii)+1); % test whether we've changed signs, phase [-pi, pi] if sign(T0) > 0 && sign(T1) > 0 T1 = T0; T0 = phase(peaks(ii)-1); baseTs = find(timeline <= timelineRS(1+peaks(ii)-1),1,'last'); else baseTs = find(timeline <= timelineRS(1+peaks(ii)),1,'last'); end % interpolate up the minimum [~,minInd] = min(abs(linspace(T0,T1,DSfactor))); ThetaPeaksTs(ii) = baseTs + minInd; end % THETA TROUGHS: find the zeros of the theta cycle as a 3 point min upon % the rectified signal troughs = zeros(length(phase),1); counter = 0; for ii = 2:length(phase)-1 [~,ind] = min(abs(phase(ii-1:ii+1) - pi)); if ind == 2 counter = counter + 1; troughs(counter) = ii; end end troughs = troughs(1:counter); % cycle through all 3 point mins found and interpolate min up to sampling % rate of the original signal ThetaTroughsTs = zeros(counter,1); for ii = 1:counter T0 = phase(troughs(ii)); T1 = phase(troughs(ii)+1); % test whether we've changed signs, phase [-pi, pi] if sign(T0) > 0 && sign(T1) > 0 T1 = T0; T0 = phase(troughs(ii)-1); baseTs = find(timeline <= timelineRS(1+troughs(ii)-1),1,'last'); else baseTs = find(timeline <= timelineRS(1+troughs(ii)),1,'last'); end % interpolate up the minimum [~,minInd] = min(abs(linspace(T0,2*pi+T1,DSfactor)-pi)); ThetaTroughsTs(ii) = baseTs + minInd; end % Now let's use the peaks and troughs found using the reconstructed signal % to do a local max/min search on the smoothed lfp signal smoothlfp = spline(timelineRS,segmentRS,timelineEpoch); ThetaPeaks = zeros(size(ThetaPeaksTs)); for ii = 1:length(ThetaPeaksTs) % find nearest trough on both sides of this peak lefttrough = find(ThetaTroughsTs < ThetaPeaksTs(ii),1,'last'); righttrough = find(ThetaTroughsTs > ThetaPeaksTs(ii),1,'first'); % We want the number of samples on each side of the peak that % approximate +/- pi/2 % check ends if isempty(lefttrough) samples = floor((ThetaTroughsTs(righttrough) - ... ThetaPeaksTs(ii))/2); elseif isempty(righttrough) samples = floor((ThetaPeaksTs(ii) - ... ThetaTroughsTs(lefttrough))/2); else samples = floor((ThetaTroughsTs(righttrough) - ... ThetaTroughsTs(lefttrough))/4); end % Access local segment around peak, find max sample, and align % subtract off start index of Epoch from ThetaPeaksTs Peak = ThetaPeaksTs(ii) - startInd + 1; samplesL = samples; if Peak-samples < 1 samplesL = Peak - 1; end samplesR = samples; if Peak+samples > Nsamples samplesR = Nsamples - Peak; end localsegment = smoothlfp(Peak-samplesL:Peak+samplesR); [~,maxInd] = max(localsegment); ThetaPeaks(ii) = ThetaPeaksTs(ii) + maxInd - samples - 1; end ThetaTroughs = zeros(size(ThetaTroughsTs)); for ii = 1:length(ThetaTroughsTs) % find nearest peak on both sides of this peak leftpeak = find(ThetaPeaksTs < ThetaTroughsTs(ii),1,'last'); rightpeak = find(ThetaPeaksTs > ThetaTroughsTs(ii),1,'first'); % We want the number of samples on each side of the peak that % approximate +/- pi/2 % check ends if isempty(leftpeak) samples = floor((ThetaPeaksTs(rightpeak) - ... ThetaTroughsTs(ii))/2); elseif isempty(rightpeak) samples = floor((ThetaTroughsTs(ii) - ... ThetaPeaksTs(leftpeak))/2); else samples = floor((ThetaPeaksTs(rightpeak) - ... ThetaPeaksTs(leftpeak))/4); end % Access local segment around peak, find max sample, and align % subtract off start index of Epoch from ThetaPeaksTs Trough = ThetaTroughsTs(ii) - startInd + 1; samplesL = samples; if Trough-samples < 1 samplesL = Trough - 1; end samplesR = samples; if Trough+samples > Nsamples samplesR = Nsamples - Trough; end localsegment = smoothlfp(Trough-samplesL:Trough+samplesR); [~,minInd] = min(localsegment); ThetaTroughs(ii) = ThetaTroughsTs(ii) + minInd - samplesL - 1; end % Log per trial data Peaks{tt} = ThetaPeaks/SR; Troughs{tt} = ThetaTroughs/SR; % power is zscored and then shifted such that the minimum is zero temp = zscore(power); Powers{tt} = temp - min(temp); end % Log per channel data Theta(cc).peaks = Peaks; Theta(cc).troughs = Troughs; Theta(cc).power = Powers; end % attach parameter log Theta.params = params; % write output to file if requested if SAVE save(sprintf('%s%s%s_theta.mat',pathname,filesep,filename),'Theta') end % write out log file logAnalysisFile(mfilename('fullpath')); end %EOF