From 8f2645d053826fc34656ad5bd6930b6a766772db Mon Sep 17 00:00:00 2001 From: lkaraba Date: Fri, 29 May 2026 13:50:54 -0400 Subject: [PATCH] add frequency filter option to remove events in ripple band --- BARRs/BARR_PSTH.m | 4 +-- BARRs/detectBARR.m | 12 ++++++++- BARRs/pareBARRs.m | 63 +++++++++++++++++++++++++++++++++++++++++++--- 3 files changed, 73 insertions(+), 6 deletions(-) diff --git a/BARRs/BARR_PSTH.m b/BARRs/BARR_PSTH.m index a7dda0d..a8da2a0 100644 --- a/BARRs/BARR_PSTH.m +++ b/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 diff --git a/BARRs/detectBARR.m b/BARRs/detectBARR.m index 89ffcaf..8b7a464 100644 --- a/BARRs/detectBARR.m +++ b/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; diff --git a/BARRs/pareBARRs.m b/BARRs/pareBARRs.m index 23caa21..5601226 100644 --- a/BARRs/pareBARRs.m +++ b/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)