MorletWaveletFilter timeseries dtype bug
Closed this issue · 3 comments
esolomon commented
The following code doesn't work with the latest version of cmlreaders (0.7.2) and ptsa (2.0.0):
import numpy as np
from cmlreaders import CMLReader, get_data_index
df = get_data_index("r1")
s = 'R1111M'
exp = 'FR1'
sess = 0
sessions = df[np.logical_and(df["subject"] == s, df['experiment']==exp)]['session'].unique()
pairs = CMLReader(s, exp).load("pairs")
reader = CMLReader(s, exp, sess)
evs = reader.load("events")
evs = evs[evs.type=='COUNTDOWN_START']
eeg = reader.load_eeg(events=evs, rel_start=0, rel_stop=10000, scheme=pairs)
eeg = eeg.to_ptsa()
#Extract spectral power
from ptsa.data.filters import MorletWaveletFilter
wf = MorletWaveletFilter(timeseries=eeg, freqs=np.array(range(70, 171, 10)), width=5, output='power')
powers = wf.filter()
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-30-67a550edab27> in <module>()
21 from ptsa.data.filters import MorletWaveletFilter
22 wf = MorletWaveletFilter(timeseries=eeg, freqs=np.array(range(70, 171, 10)), width=5, output='power')
---> 23 powers = wf.filter()
~/anaconda3/envs/CML/lib/python3.6/site-packages/ptsa/data/filters/morlet.py in filter(self)
116 mt.set_output_type(morlet.COMPLEX)
117
--> 118 mt.set_signal_array(timeseries_reshaped)
119 mt.set_wavelet_pow_array(powers_reshaped)
120 mt.set_wavelet_phase_array(phases_reshaped)
~/anaconda3/envs/CML/lib/python3.6/site-packages/ptsa/extensions/morlet/morlet.py in set_signal_array(self, signal_array)
259
260 def set_signal_array(self, signal_array):
--> 261 return _morlet.MorletWaveletTransformMP_set_signal_array(self, signal_array)
262
263 def set_wavelet_pow_array(self, wavelet_pow_array):
TypeError: Array of type 'double' required. Array of type 'short' given
Looks like the timeseries dtype either needs to be output by cmlreaders to_ptsa() as float64 or reset to float64 by PTSA.
mivade commented
Workaround for the time being:
eeg = eeg.to_ptsa().astype(np.double)
mivade commented
The easiest solution is to call .astype(np.double)
here when setting the timeseries on the filter. Otherwise, the underlying C++ code would have to be reworked significantly to explicitly handle typecasting or use templates to provide for different types.
mivade commented
On second thought, it looks like the old PTSA EEG readers always converted to double
, so maybe it would make more sense for to_ptsa
to do this to avoid other issues.