add frequency filter option to remove events in ripple band

BARRs/BARR_PSTH.m

@@ -241,7 +241,7 @@ function BARR_PSTH(plotSave,useState)
         end
 
     else
-        warning(['No pyramidal cells in ' br]);
+        warning(strcat('No pyramidal cells in ', br));
     end
     subNum = subNum+1; % should max at 4
 end
@@ -312,7 +312,7 @@ function BARR_PSTH(plotSave,useState)
         end
 
     else
-        warning(['No interneurons in ' br]);
+        warning(strcat('No interneurons in ', br));
         intEmpt = intEmpt+1;
     end
     subNum = subNum+1; % should max at 4

BARRs/detectBARR.m

@@ -53,6 +53,12 @@
 %               Default: 0
 % pareDur:      Minimum duration of BARRs kept, in seconds. This should
 %               generally be kept at or below 0.2. Default: 0.2
+% freqFilt:     Frequency band to restrict. Recommended to restrict
+%               ~100-200 Hz to exclude ripple events. If an event has this
+%               instantaneous frequency at the detected peak, it will be
+%               removed. Recommended [100 200]. Default: [] (skips)
+% freqChan:     1-based channel to use for restricting events with an
+%               instaneous peak frequency in the restricted band. 
 % zeroRip:      Logical option to remove BARRs which overlap with ripples.
 %               This seems to be a better option than remRip. 
 %               Default: false
@@ -99,6 +105,8 @@
 addParameter(p, 'spkHz', 100, @isnumeric);
 addParameter(p, 'unMax', 0, @isnumeric);
 addParameter(p, 'pareDur', 0.2, @isnumeric);
+addParameter(p, 'freqFilt',[],@isnumeric);
+addParameter(p, 'freqChan', 0, @isnumeric);
 addParameter(p, 'zeroRip', true, @islogical);
 addParameter(p, 'remRip', false, @islogical);
 
@@ -116,6 +124,8 @@
 spkHz = p.Results.spkHz;
 unMax = p.Results.unMax;
 pareDur = p.Results.pareDur;
+freqFilt = p.Results.freqFilt;
+freqChan = p.Results.freqChan;
 zeroRip = p.Results.zeroRip;
 remRip = p.Results.remRip;
 
@@ -179,7 +189,7 @@
                 'Notes',note_all,'sstd',-1*(nSigma-0.5),'estd',(nSigma-0.5),...
                 'recordMetrics',true,'remRip',remRip);
 
-HSE = pareBARRs(basepath, HSE, spikes, savePath, unMin, spkNum, pareDur, spkHz, zeroRip, unMax);
+HSE = pareBARRs(basepath, HSE, spikes, savePath, unMin, spkNum, pareDur, spkHz, freqFilt, freqChan, zeroRip, unMax);
 %add unitsForDetection parameters
 HSE.detectorinfo.Hz = Hz;
 HSE.detectorinfo.ft = ft;

BARRs/pareBARRs.m

@@ -1,5 +1,5 @@
 %% Pare events
-function HSE = pareBARRs(basepath, HSE, spikes, savePath, numCell, numSpk, pareDur, singleHz, stim, maxCell)
+function HSE = pareBARRs(basepath, HSE, spikes, savePath, numCell, numSpk, pareDur, singleHz, freqFilt, freqChan, stim, maxCell)
 %% Refine BARR detection
 %
 % This script pares back BARR detections based on restrictions provided as
@@ -25,6 +25,13 @@
 % singleHz:     Firing rate needed for a single BARR unit to be allowed to
 %               drive a BARR. Note that this does not mean it's the only 
 %               cell firing, just that it's the only flagged cell firing.
+% freqFilt:     Frequency band to restrict. Recommended to restrict
+%               ~100-200 Hz to exclude ripple events. If an event has this
+%               instantaneous frequency at the detected peak, it will be
+%               removed. Recommended [100 200]. Default: [] (skips)
+% freqChan:     1-based channel to use for restricting events with an
+%               instaneous peak frequency in the restricted band. Default
+%               will try to pull from ripples file, requires lfp file. 
 %% OPTIONAL %%
 % stim:         Logical option to indicate whether or not the session
 %               includes stimulation (ie ripple generation). If so, 
@@ -60,10 +67,14 @@
 %        HSE.timestamps(HSE.keep(HSE.NREM),:)
 %%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-if nargin < 8
+if nargin < 9 
+    freqChan = 0;
     stim = 0;
     maxCell = 0;
-elseif nargin <9
+elseif nargin < 10 %options not implemented into main script
+    stim = 0;
+    maxCell = 0;
+elseif nargin < 11
     maxCell = 0;
 end
 
@@ -113,6 +124,52 @@
     HSE.note = "HSE.keep has removed barrages overlapping with ripples";
 end
 
+if ~isempty(freqFilt)
+    if freqChan==0
+        load([basepath '\' basename '.ripples.events.mat']);
+        try
+            freqChan = ripples.detectorinfo.detectionchannel1;
+        catch
+            try
+                freqChan = ripples.detectorinfo.detectionparms.Channels(1);
+            catch
+                disp('Check for this specific case');
+            end
+        end
+    end
+    disp(['Using channel ' num2str(freqChan) ' for frequency filter']);
+    load([basepath '\' basename '.session.mat']);
+    disp('Loading LFP...');
+    SR = 1250;
+    % Compute instantaneous phase and amplitude
+    lfp = LoadBinary([basepath '\' basename '.lfp'],'frequency',SR, 'nChannels', session.extracellular.nChannels, 'channels',freqChan);
+    %     lfpfilt = bz_Filter(lfp(:, 1), 'passband', ripBP, 'filter', 'fir1', 'nyquist', SR/2);
+    h = hilbert(lfp);
+    phase = angle(h);
+    amplitude = abs(h);
+    unwrapped = unwrap(phase);
+    % Compute instantaneous frequency
+    frequency = bz_Diff(medfilt1(unwrapped, 12), (0:length(lfp) - 1)'/SR, 'smooth', 0);
+    frequency = frequency / (2 * pi);
+
+    t = linspace(0,session.epochs{end}.stopTime,length(lfp));
+%t = linspace(0,(1/SR)*length(lfp), length(lfp));
+tempKeep4 = [];
+toCheckFreq = [];
+for i = 1:length(HSE.keep)
+    [~,keep4_t] = min(abs(HSE.peaks(HSE.keep(i))-t));
+    toCheckFreq(i) = frequency(keep4_t);
+    if ~((frequency(keep4_t) < freqFilt(2)) && (frequency(keep4_t) > freqFilt(1)))
+        tempKeep4 = [tempKeep4 HSE.keep(i)];
+    end
+end
+HSE.keep = tempKeep4;
+if isfield(HSE, 'note')
+    HSE.note = "Frequency filter used for HSE.keep";
+else
+    HSE.note = [HSE.note "Frequency filter used"];
+end
+end
 
 load(strcat(basepath,filesep,basename,'.SleepState.states.mat'));
 if ~isempty(SleepState.ints.NREMstate)