'''
Functions for computing correcting fluorescence signals for changes in
baseline activity.
Authors:
- Scott C Lowe
'''
import numpy as np
import scipy.signal
[docs]def findBaselineF0(rawF, fs, axis=0, keepdims=False):
"""Find the baseline for a fluorescence imaging trace line.
The baseline, F0, is the 5th-percentile of the 1Hz
lowpass filtered signal.
Parameters
----------
rawF : array_like
Raw fluorescence signal.
fs : float
Sampling frequency of rawF, in Hz.
axis : int, optional
Dimension which contains the time series. Default is 0.
keepdims : bool, optional
Whether to preserve the dimensionality of the input. Default is
`False`.
Returns
-------
baselineF0 : numpy.ndarray
The baseline fluorescence of each recording, as an array.
Note
----
In typical usage, the input rawF is expected to be sized
`(numROI, numTimePoints, numRecs)`
and the output will then be sized `(numROI, 1, numRecs)`
if `keepdims` is `True`.
"""
# Parameters --------------------------------------------------------------
nfilt = 30 # Number of taps to use in FIR filter
fw_base = 1 # Cut-off frequency for lowpass filter, in Hz
base_pctle = 5 # Percentile to take as baseline value
# Main --------------------------------------------------------------------
# Ensure array_like input is a numpy.ndarray
rawF = np.asarray(rawF)
# Remove the first datapoint, because it can be an erroneous sample
rawF = np.split(rawF, [1], axis)[1]
if fs <= fw_base:
# If our sampling frequency is less than our goal with the smoothing
# (sampling at less than 1Hz) we don't need to apply the filter.
filtered_f = rawF
else:
# The Nyquist rate of the signal is half the sampling frequency
nyq_rate = fs / 2.0
# Cut-off needs to be relative to the nyquist rate. For sampling
# frequencies in the range from our target lowpass filter, to
# twice our target (i.e. the 1Hz to 2Hz range) we instead filter
# at the Nyquist rate, which is the highest possible frequency to
# filter at.
cutoff = min(1.0, fw_base / nyq_rate)
# Make a set of weights to use with our taps.
# We use an FIR filter with a Hamming window.
b = scipy.signal.firwin(nfilt, cutoff=cutoff, window='hamming')
# The default padlen for filtfilt is 3 * nfilt, but in case our
# dataset is small, we need to make sure padlen is not too big
padlen = min(3 * nfilt, rawF.shape[axis] - 1)
# Use filtfilt to filter with the FIR filter, both forwards and
# backwards.
filtered_f = scipy.signal.filtfilt(b, [1.0], rawF, axis=axis,
padlen=padlen)
# Take a percentile of the filtered signal
baselineF0 = np.percentile(filtered_f, base_pctle, axis=axis,
keepdims=keepdims)
return baselineF0