Skip to content

Commit

Permalink
adjusted beer lambert, freq filter, stim to xarray
Browse files Browse the repository at this point in the history
  • Loading branch information
thomasfischer11 committed Dec 3, 2023
1 parent 43abb10 commit 8178f2a
Show file tree
Hide file tree
Showing 2 changed files with 20 additions and 2 deletions.
16 changes: 15 additions & 1 deletion src/cedalion/accessors.py
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,10 @@ def freq_filter(self, fmin, fmax, butter_order=4):
array = self._obj

fny = array.cd.sampling_rate / 2
b, a = scipy.signal.butter(butter_order, (fmin / fny, fmax / fny), "bandpass")
if fmin == 0:
b, a = scipy.signal.butter(butter_order, fmax / fny, "lowpass")
else:
b, a = scipy.signal.butter(butter_order, (fmin / fny, fmax / fny), "bandpass")

if (units := array.pint.units) is not None:
array = array.pint.dequantify()
Expand Down Expand Up @@ -86,3 +89,14 @@ def rename_events(self, rename_dict):
stim = self._obj
for old_trial_type, new_trial_type in rename_dict.items():
stim.loc[stim.trial_type == old_trial_type, "trial_type"] = new_trial_type

def conditions(self):
return self._obj.trial_type.unique()

def to_xarray(self, time : xr.DataArray):
stim = self._obj
conds = self.conditions()
stim_arr = xr.DataArray(np.zeros((time.shape[0], len(conds))), dims=["time", "condition"], coords={"time" : time, "condition" : conds})
for index, row in stim.iterrows():
stim_arr.loc[row.onset, row.trial_type] = 1
return stim_arr
6 changes: 5 additions & 1 deletion src/cedalion/nirs.py
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,7 @@ def beer_lambert(
geo3d: xr.DataArray,
dpf: xr.DataArray,
spectrum: str = "prahl",
calc_od: bool = True
):
"""Calculate concentration changes from amplitude data."""
validators.has_channel(amplitudes)
Expand All @@ -69,7 +70,10 @@ def beer_lambert(
dists = channel_distances(amplitudes, geo3d)
dists = dists.pint.to("mm")

optical_density = -np.log(amplitudes / amplitudes.mean("time"))
if calc_od:
optical_density = -np.log(amplitudes / amplitudes.mean("time"))
else:
optical_density = amplitudes
# conc = Einv @ (optical_density / ( dists * dpf))
conc = xr.dot(Einv, optical_density / (dists * dpf), dims=["wavelength"])
conc = conc.pint.to("micromolar")
Expand Down

0 comments on commit 8178f2a

Please sign in to comment.