Repository: esi-neuroscience/syncopy

Issue:
Wavelet Cross-Spectra
We now have ~~one~~ two interested parties in Wavelet Cross-Spectra (and time-dependent Granger I assume), and our `CrossSpectralData` already has a time axis so this feature addition should be rather painless.

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index ba408ebf..1a4549cf 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,7 @@ and this project follows [Semantic Versioning](https://semver.org/spec/v2.0.0.ht
 ## [Unreleased]
 
 ### NEW
+- Added time dependent coherence analysis #287
 - Added basic statistics (mean, std, var and median) for Syncopy data objects
 - Added additional .info entries for Granger analysis, indicating details about the computation.
 - Added additional .info entries for FoooF results, e.g. Gaussian fit parameters.
diff --git a/syncopy/datatype/methods/statistics.py b/syncopy/datatype/methods/statistics.py
index 7200b19a..126b30ab 100644
--- a/syncopy/datatype/methods/statistics.py
+++ b/syncopy/datatype/methods/statistics.py
@@ -342,7 +342,7 @@ class NumpyStatDim(ComputationalRoutine):
 
         # trial average, needs equal trial lengths.. just copy from 1st
         elif dim != 'time' and not self.keeptrials:
-            trldef = in_data.trialdefinition[0, :][None, :]
+            trldef = in_data.selection.trialdefinition[0, :][None, :]
 
         # each trial has empty time axis, so we attach trivial trialdefinition:
         # 1 sample per trial for stacking
diff --git a/syncopy/nwanalysis/AV_compRoutines.py b/syncopy/nwanalysis/AV_compRoutines.py
index 1f9ece05..a92ec858 100644
--- a/syncopy/nwanalysis/AV_compRoutines.py
+++ b/syncopy/nwanalysis/AV_compRoutines.py
@@ -52,7 +52,7 @@ def normalize_csd_cF(csd_av_dat,
 
     Parameters
     ----------
-    csd_av_dat : (1, nFreq, N, N) :class:`numpy.ndarray`
+    csd_av_dat : (nTime, nFreq, N, N) :class:`numpy.ndarray`
         Cross-spectral densities for `N` x `N` channels
         and `nFreq` frequencies averaged over trials.
     output : {'abs', 'pow', 'fourier', 'real', 'imag'}, default: 'abs'
@@ -69,7 +69,7 @@ def normalize_csd_cF(csd_av_dat,
 
     Returns
     -------
-    CS_ij : (1, nFreq, N, N) :class:`numpy.ndarray`
+    CS_ij : (nTime, nFreq, N, N) :class:`numpy.ndarray`
         Coherence for all channel combinations ``i,j``.
         `N` corresponds to number of input channels.
 
@@ -107,10 +107,9 @@ def normalize_csd_cF(csd_av_dat,
     if noCompute:
         return outShape, fmt
 
-    CS_ij = normalize_csd(csd_av_dat[0], output)
+    CS_ij = normalize_csd(csd_av_dat, output)
 
-    # re-attach dummy time axis
-    return CS_ij[None, ...]
+    return CS_ij
 
 
 class NormalizeCrossSpectra(ComputationalRoutine):
@@ -157,7 +156,9 @@ class NormalizeCrossSpectra(ComputationalRoutine):
 
     def process_metadata(self, data, out):
 
-        propagate_properties(data, out, self.keeptrials)
+        time_axis = np.any(np.diff(data.trialdefinition)[:,0] != 1)
+
+        propagate_properties(data, out, self.keeptrials, time_axis)
         out.freq = data.freq
 
 
diff --git a/syncopy/nwanalysis/ST_compRoutines.py b/syncopy/nwanalysis/ST_compRoutines.py
index 562d1642..4c4a667c 100644
--- a/syncopy/nwanalysis/ST_compRoutines.py
+++ b/syncopy/nwanalysis/ST_compRoutines.py
@@ -23,7 +23,7 @@ from syncopy.shared.kwarg_decorators import process_io
 
 
 @process_io
-def spectral_dyadic_product_cF(trl_dat,
+def spectral_dyadic_product_cF(specs,
                                chunkShape=None,
                                noCompute=False):
     """
@@ -36,8 +36,8 @@ def spectral_dyadic_product_cF(trl_dat,
 
     Parameters
     ----------
-    trl_dat: (K, N) :class:`numpy.ndarray`
-        Complex and frequency aligned multi-channel spectral data.
+    specs: (K, N) :class:`numpy.ndarray`
+        Complex and frequency aligned multi-channel single-trial spectral data.
         The 1st dimension is interpreted as the frequency axis,
         `N` columns represent individual channels.
     noCompute : bool
@@ -68,25 +68,24 @@ def spectral_dyadic_product_cF(trl_dat,
     """
 
     # default dimord for SpectralData is ['time', 'taper', 'freq', 'channel']
-    nFreq = trl_dat.shape[2]
-    nChannels = trl_dat.shape[3]
+    nTime = specs.shape[0]
+    nFreq = specs.shape[2]
+    nChannels = specs.shape[3]
 
     # we always average over tapers here
-    outShape = (1, nFreq, nChannels, nChannels)
+    outShape = (nTime, nFreq, nChannels, nChannels)
 
     # cross spectra are complex, input gets checked in frontend!
     if noCompute:
         return outShape, spectralDTypes["fourier"]
 
-    # split off empty time axis, shape is now (nTapers, nFreq, nChannel)
-    specs = trl_dat[0, ...]
     # dyadic product along channel axes
-    # result has shape (nTapers x nFreq x nChannels x nChannels)
-    CS_ij = specs[:, :, :, np.newaxis] * specs[:, :, np.newaxis, :].conj()
+    # result has shape (nTime, nTapers x nFreq x nChannels x nChannels)
+    CS_ij = specs[..., np.newaxis] * specs[..., np.newaxis, :].conj()
 
-    # now average tapers and attach dummy time axis
-    # result has shape (1 x nFreq x nChannels x nChannels)
-    CS_ij = CS_ij.mean(axis=0)[np.newaxis, ...]
+    # now average tapers
+    # result has shape (nTime x nFreq x nChannels x nChannels)
+    CS_ij = CS_ij.mean(axis=1)
     return CS_ij
 
 
@@ -119,9 +118,16 @@ class SpectralDyadicProduct(ComputationalRoutine):
 
     def process_metadata(self, data, out):
 
-        propagate_properties(data, out, self.keeptrials)
+        # time dependent coherence needs SpectralData input
+        if 'Spectral' in data.__class__.__name__:
+            time_axis = np.any(np.diff(data.trialdefinition)[:,0] != 1)
+        else:
+            time_axis = False
 
+        propagate_properties(data, out, self.keeptrials, time_axis)
+        out.freq = data.freq
 
+        
 @process_io
 def cross_spectra_cF(trl_dat,
                      samplerate=1,
diff --git a/syncopy/nwanalysis/connectivity_analysis.py b/syncopy/nwanalysis/connectivity_analysis.py
index 58fce576..56acec44 100644
--- a/syncopy/nwanalysis/connectivity_analysis.py
+++ b/syncopy/nwanalysis/connectivity_analysis.py
@@ -202,14 +202,6 @@ def connectivityanalysis(data, method="coh", keeptrials=False, output="abs",
         raise SPYValueError(lgl, 'data', act)
     timeAxis = data.dimord.index("time")
 
-    # check that for SpectralData input, we have empty time axes
-    # no time-resolved connectivity supported atm
-    if isinstance(data, SpectralData):
-        if data.data.shape[data.dimord.index('time')] != len(data.trials):
-            lgl = "line spectra without time axis (mtmfft or time average)"
-            act = "time-frequency result"
-            raise SPYValueError(lgl, 'data', act)
-
     # Get everything of interest in local namespace
     defaults = get_defaults(connectivityanalysis)
     lcls = locals()
@@ -323,6 +315,13 @@ def connectivityanalysis(data, method="coh", keeptrials=False, output="abs",
                 act = "real valued spectral data"
                 raise SPYValueError(lgl, 'data', act)
 
+            if method == 'granger':
+                # check that for SpectralData input, we have empty time axes
+                # no time-resolved Granger supported atm
+                if isinstance(data, SpectralData):
+                    if data.data.shape[data.dimord.index('time')] != len(data.trials):
+                        raise NotImplementedError("Time resolved Granger causality from tf-spectra not available atm")
+
             # by constraining to output='fourier', detrimental taper averaging
             # gets already catched by freqanalysis!
 
@@ -385,7 +384,8 @@ def connectivityanalysis(data, method="coh", keeptrials=False, output="abs",
     # now take the trial average from the single trial CR as input
     av_compRoutine.initialize(st_out, out._stackingDim, chan_per_worker=None)
     av_compRoutine.pre_check()   # make sure we got a trial_average
-    av_compRoutine.compute(st_out, out, parallel=False, log_dict=log_dict)
+    av_compRoutine.compute(st_out, out, parallel=kwargs.get("parallel"),
+                           log_dict=log_dict)
 
     # attach potential older cfg's from the input
     # to support chained frontend calls..
diff --git a/syncopy/nwanalysis/csd.py b/syncopy/nwanalysis/csd.py
index d7adedde..905c87c0 100644
--- a/syncopy/nwanalysis/csd.py
+++ b/syncopy/nwanalysis/csd.py
@@ -133,7 +133,7 @@ def normalize_csd(csd_av_dat,
 
     Parameters
     ----------
-    csd_av_dat : (nFreq, N, N) :class:`numpy.ndarray`
+    csd_av_dat : (nTime, nFreq, N, N) :class:`numpy.ndarray`
         Averaged cross-spectral densities for `N` x `N` channels
         and `nFreq` frequencies averaged over trials.
     output : {'abs', 'pow', 'fourier', 'complex', 'angle', 'imag' or 'real'}, default: 'abs'
@@ -146,7 +146,7 @@ def normalize_csd(csd_av_dat,
 
     Returns
     -------
-    CS_ij : (nFreq, N, N) :class:`numpy.ndarray`
+    CS_ij : (nTime, nFreq, N, N) :class:`numpy.ndarray`
         Coherence for all channel combinations ``i,j``.
         `N` corresponds to number of input channels.
 
@@ -156,17 +156,16 @@ def normalize_csd(csd_av_dat,
           data using the imaginary part of coherency."
           Clinical neurophysiology 115.10 (2004): 2292-2307.
     """
-    # re-shape to (nChannels x nChannels x nFreq)
-    CS_ij = csd_av_dat.transpose(1, 2, 0)
 
-    # main diagonal has shape (nFreq x nChannels): the auto spectra
-    diag = CS_ij.diagonal()
+    CS_ij = csd_av_dat.view()
+
+    # channel diagonal has shape (nTime x nFreq x nChannels): the auto spectra
+    diag = CS_ij.diagonal(axis1=-2, axis2=-1)
 
     # get the needed product pairs of the autospectra
-    Ciijj = np.sqrt(diag[:, :, None] * diag[:, None, :]).T
+    Ciijj = np.sqrt(diag[..., None] * diag[..., None, :])
     CS_ij = CS_ij / Ciijj
 
     CS_ij = spectralConversions[output](CS_ij)
 
-    # re-shape to original form
-    return CS_ij.transpose(2, 0, 1)
+    return CS_ij
diff --git a/syncopy/plotting/_helpers.py b/syncopy/plotting/_helpers.py
index 87f3a957..dd1d82f6 100644
--- a/syncopy/plotting/_helpers.py
+++ b/syncopy/plotting/_helpers.py
@@ -150,11 +150,11 @@ def get_method(dataobject):
 def get_output(dataobject):
 
     """
-    Returns the method string from
+    Returns the output string from
     the log of a Syncopy data object
     """
 
-    # get the method string in a capture group
+    # get the output string in a capture group
     pattern = re.compile(r'[\s\w\D]+output = (\w+)')
     match = pattern.match(dataobject._log)
     if match:
@@ -197,3 +197,14 @@ def calc_multi_layout(nAx):
         nrows, ncols = 1, 2
 
     return nrows, ncols
+
+
+def check_if_time_freq(data):
+    """
+    Looks into the first column of the trialdefinition
+    to determine if there is a real time axis, or it is
+    just trial stacking.
+    """
+    is_tf = np.any(np.diff(data.trialdefinition)[:, 0] != 1)
+
+    return is_tf
diff --git a/syncopy/plotting/_plotting.py b/syncopy/plotting/_plotting.py
index 18d1c16d..2b7d1777 100644
--- a/syncopy/plotting/_plotting.py
+++ b/syncopy/plotting/_plotting.py
@@ -183,5 +183,6 @@ def plot_tfreq(ax, data_yx, times, freqs, **pkwargs):
     extent = [times[0], times[-1],
               freqs[0] - df / 2, freqs[-1] - df / 2]
 
-    ax.imshow(data_yx[::-1], aspect='auto', cmap=pltConfig['cmap'],
+    cmap = pkwargs.pop('cmap', pltConfig['cmap'])
+    ax.imshow(data_yx[::-1], aspect='auto', cmap=cmap,
               extent=extent, **pkwargs)
diff --git a/syncopy/plotting/mp_plotting.py b/syncopy/plotting/mp_plotting.py
index 155b30fb..1106cd14 100644
--- a/syncopy/plotting/mp_plotting.py
+++ b/syncopy/plotting/mp_plotting.py
@@ -115,8 +115,8 @@ def plot_SpectralData(data, **show_kwargs):
     elif len(data.trials) == 1:
         trl = 0
 
-    labels = plot_helpers.parse_channel(data, show_kwargs)
-    nAx = 1 if isinstance(labels, str) else len(labels)
+    channels = plot_helpers.parse_channel(data, show_kwargs)
+    nAx = 1 if isinstance(channels, str) else len(channels)
 
     if nAx < 2:
         SPYWarning("Please select at least two channels for a multipanelplot!")
@@ -128,9 +128,10 @@ def plot_SpectralData(data, **show_kwargs):
         # determine axes layout, prefer columns over rows due to display aspect ratio
         nrows, ncols = plot_helpers.calc_multi_layout(nAx)
 
-    # how got the spectrum computed
-    method = plot_helpers.get_method(data)
-    if method in ('wavelet', 'superlet', 'mtmconvol'):
+    # -- check if it is a time-frequency spectrum ----------
+    is_tf = np.any(np.diff(data.trialdefinition)[:, 0] != 1)
+    # ------------------------------------------------------
+    if is_tf:
         fig, axs = _plotting.mk_multi_img_figax(nrows, ncols)
 
         # this could be more elegantly solve by
@@ -147,7 +148,7 @@ def plot_SpectralData(data, **show_kwargs):
             raise SPYValueError(lgl, varname="show_kwargs", actual=act)
 
         maxP = data_cyx.max()
-        for data_yx, ax, label in zip(data_cyx, axs.flatten(), labels):
+        for data_yx, ax, label in zip(data_cyx, axs.flatten(), channels):
             _plotting.plot_tfreq(ax, data_yx, time, freqs, vmax=maxP)
             ax.set_title(label, fontsize=pltConfig['mTitleSize'])
         fig.tight_layout()
@@ -182,11 +183,22 @@ def plot_SpectralData(data, **show_kwargs):
             data_y = data.show(**show_kwargs)
             ylabel = f'{output}'
 
+        taper_labels = None
+        if len(data.taper) != 1:   
+            taper = show_kwargs.get('taper')
+            # multiple tapers are to be plotted
+            if not isinstance(taper, (Number, str)):
+                taper_labels = data.taper
+
         fig, axs = _plotting.mk_multi_line_figax(nrows, ncols, xlabel='frequency (Hz)',
                                                  ylabel=ylabel)
 
-        for chan_dat, ax, label in zip(data_y.T, axs.flatten(), labels):
-            _plotting.plot_lines(ax, data_x, chan_dat, label=label, leg_fontsize=pltConfig['mLegendSize'])
+        for chan_dat, ax, label in zip(data_y.T, axs.flatten(), channels):
+            if taper_labels is not None:
+                _plotting.plot_lines(ax, data_x, chan_dat, label=taper_labels, leg_fontsize=pltConfig['mLegendSize'])
+            else:
+                _plotting.plot_lines(ax, data_x, chan_dat)
+            ax.set_title(label, fontsize=pltConfig['mTitleSize'])
 
         # delete empty plot due to grid extension
         # because of prime nAx -> can be maximally 1 plot
diff --git a/syncopy/plotting/sp_plotting.py b/syncopy/plotting/sp_plotting.py
index 078def8f..2d168d48 100644
--- a/syncopy/plotting/sp_plotting.py
+++ b/syncopy/plotting/sp_plotting.py
@@ -100,11 +100,9 @@ def plot_SpectralData(data, **show_kwargs):
     elif len(data.trials) == 1:
         trl = 0
 
-    # -- how got the spectrum computed --
-    method = plot_helpers.get_method(data)
-    # -----------------------------------
+    is_tf = plot_helpers.check_if_time_freq(data)
 
-    if method in ('wavelet', 'superlet', 'mtmconvol'):
+    if is_tf:
         # multiple channels?
         label = plot_helpers.parse_channel(data, show_kwargs)
         # only relevant for mtmconvol
@@ -119,8 +117,6 @@ def plot_SpectralData(data, **show_kwargs):
         # here we always need a new axes
         fig, ax = _plotting.mk_img_figax()
 
-        # this could be more elegantly solve by
-        # an in-place selection?!
         time = plot_helpers.parse_toi(data, trl, show_kwargs)
         freqs = plot_helpers.parse_foi(data, show_kwargs)
 
@@ -147,6 +143,24 @@ def plot_SpectralData(data, **show_kwargs):
                    "ignoring `toi/toilim` selection!")
             SPYWarning(msg)
 
+        # multiple channels?
+        channels = plot_helpers.parse_channel(data, show_kwargs)
+
+        # just multiple tapers or multiple channels in one plot
+        if len(data.taper) != 1:
+            taper = show_kwargs.get('taper')
+            if not isinstance(taper, (Number, str)) and not isinstance(channels, str):
+                msg = "Please select a single taper or a single channel \nfor plotting multi-taper spectra.. aborting plotting\n"
+                SPYWarning(msg)
+                return None, None
+            # single channel, multiple tapers
+            elif isinstance(channels, str):
+                labels = data.taper
+            # single taper, multiple channels
+            elif isinstance(taper, (Number, str)):
+                labels = channels
+        else:
+            labels = channels
         # get the data to plot
         data_x = plot_helpers.parse_foi(data, show_kwargs)
         output = plot_helpers.get_output(data)
@@ -156,14 +170,15 @@ def plot_SpectralData(data, **show_kwargs):
             data_y = np.log10(data.show(**show_kwargs))
             ylabel = 'power (dB)'
         elif output in ['fourier', 'complex']:
-            SPYWarning("Can't plot complex valued spectra, choose 'real' or 'imag' as output. Aborting plotting.")
+            SPYWarning("Can't plot complex valued spectra, choose 'real' or 'imag' as freqanalysis output.. aborting plotting")
             return None, None
         else:
             data_y = data.show(**show_kwargs)
             ylabel = f'{output}'
 
-        # multiple channels?
-        labels = plot_helpers.parse_channel(data, show_kwargs)
+        # flip if required
+        if data_y.shape[1] == len(data_x):
+            data_y = data_y.T
 
         fig, ax = _plotting.mk_line_figax(xlabel='frequency (Hz)',
                                           ylabel=ylabel)
@@ -245,26 +260,47 @@ def plot_CrossSpectralData(data, **show_kwargs):
     else:
         raise NotImplementedError
 
-    # get the data to plot
-    data_y = data.show(**show_kwargs)
-    if data_y.size == 0:
-        lgl = "Selection with non-zero size"
-        act = "got zero samples"
-        raise SPYValueError(lgl, varname="show_kwargs", actual=act)
+    is_tf = plot_helpers.check_if_time_freq(data)
 
-    # create the axes and figure if needed
-    # persistent axes allows for plotting different
-    # channel combinations into the same figure
-    if not hasattr(data, 'fig') or not _plotting.ppl.fignum_exists(data.fig.number):
-        data.fig, data.ax = _plotting.mk_line_figax(xlabel, ylabel)
-    _plotting.plot_lines(data.ax, data_x, data_y, label=label)
-    # format axes
-    if method in ['granger', 'coh'] and output in ['pow', 'abs']:
-        data.ax.set_ylim((-.02, 1.02))
-    elif method == 'corr':
-        data.ax.set_ylim((-1.02, 1.02))
-    data.ax.legend(ncol=1)
+    # time dependent coherence
+    if method == 'coh' and is_tf:
+        # here we always need a new axes
+        fig, ax = _plotting.mk_img_figax()
 
-    data.fig.tight_layout()
+        time = plot_helpers.parse_toi(data, trl, show_kwargs)
+        freqs = plot_helpers.parse_foi(data, show_kwargs)
 
-    return data.fig, data.ax
+        # custom dimords for SpectralData not supported atm
+        # dimord is time x freq x channel_i x channel_j
+        # need freq x time for plotting
+        data_yx = data.show(**show_kwargs).T
+        _plotting.plot_tfreq(ax, data_yx, time, freqs, cmap='cividis')
+        ax.set_title(f"{method}: " + label, fontsize=pltConfig['sTitleSize'])
+        fig.tight_layout()
+
+        return fig, ax
+        
+    else:
+        # get the data to plot
+        data_y = data.show(**show_kwargs)
+        if data_y.size == 0:
+            lgl = "Selection with non-zero size"
+            act = "got zero samples"
+            raise SPYValueError(lgl, varname="show_kwargs", actual=act)
+
+        # create the axes and figure if needed
+        # persistent axes allows for plotting different
+        # channel combinations into the same figure
+        if not hasattr(data, 'fig') or not _plotting.ppl.fignum_exists(data.fig.number):
+            data.fig, data.ax = _plotting.mk_line_figax(xlabel, ylabel)
+        _plotting.plot_lines(data.ax, data_x, data_y, label=label)
+        # format axes
+        if method in ['granger', 'coh'] and output in ['pow', 'abs']:
+            data.ax.set_ylim((-.02, 1.02))
+        elif method == 'corr':
+            data.ax.set_ylim((-1.02, 1.02))
+        data.ax.legend(ncol=1)
+
+        data.fig.tight_layout()
+
+        return data.fig, data.ax
diff --git a/syncopy/shared/computational_routine.py b/syncopy/shared/computational_routine.py
index 180c113e..6627ff5e 100644
--- a/syncopy/shared/computational_routine.py
+++ b/syncopy/shared/computational_routine.py
@@ -24,7 +24,8 @@ import dask_jobqueue as dj
 
 # Local imports
 import syncopy as spy
-from .tools import get_defaults, check_slurm_available
+from .tools import get_defaults
+from .dask_helpers import check_slurm_available
 from syncopy import __storage__, __acme__
 from syncopy.shared.errors import SPYValueError, SPYTypeError, SPYParallelError, SPYWarning
 if __acme__:
@@ -32,6 +33,7 @@ if __acme__:
     # # In case of problems w/worker-stealing, uncomment the following lines
     # import dask
     # dask.config.set(distributed__scheduler__work_stealing=False)
+
 from syncopy.shared.metadata import parse_cF_returns, h5_add_metadata
 
 __all__ = []
@@ -677,9 +679,10 @@ class ComputationalRoutine(ABC):
             if self.chunkMem >= mem_thresh * memSize:
                 self.chunkMem /= 1024**3
                 memSize /= 1024**3
-                msg = "Single-trial result sizes ({0:2.2f} GB) larger than available " +\
-                      "memory ({1:2.2f} GB) currently not supported"
-                raise NotImplementedError(msg.format(self.chunkMem, memSize))
+                msg = ("Single-trial processing requires {0:2.2f} GB of memory "
+                       "which is larger than the available "
+                       "memory ({1:2.2f} GB)")
+                raise SPYParallelError(msg.format(2 * self.chunkMem, memSize))
 
         # The `method` keyword can be used to override the `parallel` flag
         if method is None:
@@ -843,9 +846,10 @@ class ComputationalRoutine(ABC):
             if self.chunkMem >= self.mem_thresh * workerMemMax:
                 self.chunkMem /= 1024**3
                 workerMemMax /= 1000**3
-                msg = "Single-trial result sizes ({0:2.2f} GB) larger than available " +\
-                    "worker memory ({1:2.2f} GB) currently not supported"
-                raise SPYParallelError(msg.format(self.chunkMem, workerMemMax))
+                msg = ("Single-trial processing requires {0:2.2f} GB of memory "
+                       "which is larger than the available "
+                       "worker memory ({1:2.2f} GB)")
+                raise SPYParallelError(msg.format(2 * self.chunkMem, workerMemMax))
 
         # --- trigger actual computation ---
 
@@ -970,7 +974,7 @@ class ComputationalRoutine(ABC):
                 if self.keeptrials:
                     target[outgrid] = res
                 else:
-                    target[()] = np.nansum([target, res], axis=0)
+                    target[()] += res
 
                 # Flush every iteration to avoid memory leakage
                 h5fout.flush()
@@ -1057,7 +1061,7 @@ class ComputationalRoutine(ABC):
 
 # --- metadata helper functions ---
 
-def propagate_properties(in_data, out_data, keeptrials=True):
+def propagate_properties(in_data, out_data, keeptrials=True, time_axis=False):
 
     """
     Propagating data class properties (channels, trials, time, ...)
@@ -1080,7 +1084,7 @@ def propagate_properties(in_data, out_data, keeptrials=True):
 
     This function is a general approach to unify the
     propagation / manipulation of all needed properties
-    done in the `process_metadata` implementations of the 
+    done in the `process_metadata` implementations of the
     respective ComputationalRoutines.
 
     Parameters
@@ -1089,6 +1093,9 @@ def propagate_properties(in_data, out_data, keeptrials=True):
         Input object to get attributes from
     out_data : Syncopy data object
         Output object to set attributes
+    keeptrials : bool
+    time_axis : bool
+        If `True` `out_data` has a time axis (not just trial stacking)
     """
 
     # instance checkers
@@ -1096,91 +1103,74 @@ def propagate_properties(in_data, out_data, keeptrials=True):
     is_Spectral = lambda data: isinstance(data, spy.SpectralData)
     is_CrossSpectral = lambda data: isinstance(data, spy.CrossSpectralData)
 
-    # simplest case, direct propagation between two AnalogData objects
-    if is_Analog(in_data) and is_Analog(out_data):
-        print(f"Propagating properties AnalogData -> AnalogData")
 
-        # get channels and trial selections
-        if in_data.selection is not None:
-            chanSec = in_data.selection.channel
-            # captures toi selections
-            trldef = in_data.selection.trialdefinition
+    
+    # attach a dummy selection for easier propagation
+    selection_cleanup = False
+    if in_data.selection is None:
+        in_data.selectdata(inplace=True)
+        selection_cleanup = True
+
+    # if preserving the data type, propagation is straightforward
+    if in_data.__class__ == out_data.__class__:
+
+        # Get/set dimensional attributes changed by selection
+        for prop in in_data.selection._dimProps:
+            selection = getattr(in_data.selection, prop)
+            if selection is not None:
+                if np.issubdtype(type(selection), np.number):
+                    selection = [selection]
+                setattr(out_data, prop, getattr(in_data, prop)[selection])
+
+        if keeptrials:
+            out_data.trialdefinition = in_data.selection.trialdefinition
         else:
-            chanSec = slice(None)
-            trldef = in_data.trialdefinition
+            # trial average requires equal length trials, so just copy the 1st
+            out_data.trialsdefinition = in_data.trialdefinition[0, :][None, :]
 
-        out_data.trialdefinition = trldef
-        out_data.channel = np.array(in_data.channel[chanSec])
         out_data.samplerate = in_data.samplerate
+        if selection_cleanup:
+            in_data.selection = None
+        return
+
+    # --- propagate only channels and deal with the rest below---
+
+    elif is_Spectral(out_data):
+        chanSec = in_data.selection.channel
+        out_data.channel = np.array(in_data.channel[chanSec])
+
+    # from one channel to cross-channel data
+    elif (is_Analog(in_data) or is_Spectral(in_data)) and is_CrossSpectral(out_data):
+        chanSec = in_data.selection.channel
+        out_data.channel_i = np.array(in_data.channel[chanSec])
+        out_data.channel_j = np.array(in_data.channel[chanSec])
 
-    # this is for (cross-) spectral data with empty time axis (MultiTaperFFT / CSD)
-    elif (is_Analog(in_data) or is_Spectral(in_data)) and (is_Spectral(out_data) or is_CrossSpectral(out_data)):
+    # --- time and trialdefinition ---
 
-        # Some index gymnastics to get trial begin/end "samples"
-        if in_data.selection is not None:
-            chanSec = in_data.selection.channel
+    if not time_axis:
+        # Note that here the `time` axis is only(!) used as stacking dimension
+        if keeptrials:
             trldef = in_data.selection.trialdefinition
             for row in range(trldef.shape[0]):
                 trldef[row, :2] = [row, row + 1]
-        else:
-            chanSec = slice(None)
-            time = np.arange(len(in_data.trials))
-            time = time.reshape((time.size, 1))
-            trldef = np.hstack((time, time + 1,
-                                np.zeros((len(in_data.trials), 1)),
-                                np.array(in_data.trialinfo)))
-
-        # Attach constructed trialdef-array (if even necessary)
-        # Note that here the `time` axis is only(!) used as stacking dimension
-        if keeptrials:
             out_data.trialdefinition = trldef
+
         else:
             # only single trial on the time stacking dim.
             out_data.trialdefinition = np.array([[0, 1, 0]])
 
-        # Attach remaining meta-data
         out_data.samplerate = in_data.samplerate
 
-        if is_Spectral(out_data):
-            out_data.channel = np.array(in_data.channel[chanSec])
-
-        elif is_CrossSpectral(out_data):
-            out_data.channel_i = np.array(in_data.channel[chanSec])
-            out_data.channel_j = np.array(in_data.channel[chanSec])
-
-        # no further foi manipulation after original freqanalysis!
-        if is_Spectral(in_data):
-            out_data.freq = in_data.freq
-
-    # for the AV_compRoutines and statistics
-    elif is_CrossSpectral(in_data) and is_CrossSpectral(out_data):
-        # Some index gymnastics to get trial begin/end "samples"
-        if in_data.selection is not None:
-            chanSec_i = in_data.selection.channel_i
-            chanSec_j = in_data.selection.channel_j
-            trl = in_data.selection.trialdefinition
-            for row in range(trl.shape[0]):
-                trl[row, :2] = [row, row + 1]
-        else:
-            chanSec_i = slice(None)
-            chanSec_j = slice(None)
-            time = np.arange(len(in_data.trials))
-            time = time.reshape((time.size, 1))
-            trl = np.hstack((time, time + 1,
-                             np.zeros((len(in_data.trials), 1)),
-                             np.array(in_data.trialinfo)))
-
-        # Attach constructed trialdef-array (if even necessary)
+    # this captures active in-place selections on the time axis
+    else:
         if keeptrials:
-            out_data.trialdefinition = trl
+            out_data.trialdefinition = in_data.selection.trialdefinition
+        # just copy 1st trial, have to be all equal anyways
         else:
-            out_data.trialdefinition = np.array([[0, 1, 0]])
+            out_data.trialdefinition = in_data.selection.trialdefinition[0, :][None, :]
 
-        # Attach remaining meta-data
+        # this might bite us
         out_data.samplerate = in_data.samplerate
-        out_data.channel_i = np.array(in_data.channel_i[chanSec_i])
-        out_data.channel_j = np.array(in_data.channel_j[chanSec_j])
 
-    # nothing else supported atm
-    else:
-        raise NotImplementedError
+    if selection_cleanup:
+        in_data.selection = None
diff --git a/syncopy/shared/dask_helpers.py b/syncopy/shared/dask_helpers.py
new file mode 100644
index 00000000..3576ee85
--- /dev/null
+++ b/syncopy/shared/dask_helpers.py
@@ -0,0 +1,51 @@
+# -*- coding: utf-8 -*-
+#
+# Basic checkers to facilitate direct Dask interface
+#
+
+import subprocess
+from time import sleep
+
+# Syncopy import
+from syncopy.shared.errors import SPYWarning, SPYInfo
+
+
+def check_slurm_available():
+    """
+    Returns `True` if a SLURM instance could be reached via
+    a `sinfo` call, `False` otherwise.
+    """
+
+    # Check if SLURM's `sinfo` can be accessed
+    proc = subprocess.Popen("sinfo",
+                            stdout=subprocess.PIPE, stderr=subprocess.PIPE,
+                            text=True, shell=True)
+    _, err = proc.communicate()
+    # Any non-zero return-code means SLURM is not available
+    # so we disable ACME
+    if proc.returncode != 0:
+        has_slurm = False
+    else:
+        has_slurm = True
+
+    return has_slurm
+
+
+def check_workers_available(cluster):
+    """
+    Tries to see the Dask workers and waits
+    until all requested workers are available
+    """
+    
+    totalWorkers = len(cluster.requested)    
+    sec = 0
+    workers = cluster.scheduler_info['workers']
+
+    while len(workers) != totalWorkers:
+        SPYInfo(f"{len(workers)}/{totalWorkers} workers available, waiting.. {sec}s")
+        sleep(1)
+        sec += 2
+        workers = cluster.scheduler_info['workers']
+    # wait a little more to get consistent client print out
+    sleep(0.3)
+    
diff --git a/syncopy/shared/kwarg_decorators.py b/syncopy/shared/kwarg_decorators.py
index 1841f1ab..89446a31 100644
--- a/syncopy/shared/kwarg_decorators.py
+++ b/syncopy/shared/kwarg_decorators.py
@@ -16,7 +16,7 @@ from syncopy.shared.errors import (SPYTypeError, SPYValueError,
                                    SPYError, SPYWarning, SPYInfo)
 from syncopy.shared.tools import StructDict
 from syncopy.shared.metadata import h5_add_metadata, parse_cF_returns
-from .tools import check_slurm_available
+from .dask_helpers import check_slurm_available, check_workers_available
 import syncopy as spy
 
 __all__ = []
@@ -467,17 +467,8 @@ def detect_parallel_client(func):
         parallel = kwargs.get("parallel")
         kill_spawn = False
         has_slurm = check_slurm_available()
-
-        # we are on a HPC but ACME is missing, LocalCluster still got created
-        # but a warning is issued
-        if parallel is True and has_slurm and not spy.__acme__:
-            msg = ("We are apparently on a slurm cluster but Syncopy does not provide an "
-                   "automatic Dask SLURMCluster on its own!"
-                   "\nPlease consider using ACME (https://github.com/esi-neuroscience/acme)"
-                   "\nor configure your own cluster via `dask_jobqueue.SLURMCluster()`"
-                   "\n\nCreating a LocalCluster as fallback.."
-                   )
-            SPYWarning(msg)
+        # warning only emitted if slurm available but no ACME or Dask client
+        slurm_msg = ""
 
         # This effectively searches for a global dask cluster, and sets
         # parallel=True if one was found. If no cluster was found, parallel is set to False,
@@ -485,7 +476,10 @@ def detect_parallel_client(func):
         # this needs explicit `parallel=True`.
         if parallel is None:
             try:
-                dd.get_client()
+                client = dd.get_client()
+                check_workers_available(client.cluster)
+                msg = f"..attaching to running Dask client:\n{client}"
+                SPYInfo(msg)
                 parallel = True
             except ValueError:
                 parallel = False
@@ -496,9 +490,26 @@ def detect_parallel_client(func):
         elif parallel is True and (not has_slurm or not spy.__acme__):
             # if already one cluster is reachable do nothing
             try:
-                dd.get_client()
+                client = dd.get_client()
+                check_workers_available(client.cluster)
+                msg = f"..attaching to running Dask client:\n{client}"
+                SPYInfo(msg)
             except ValueError:
-                # spawn fallback local cluster
+                # we are on a HPC but ACME and Dask client are missing,
+                # LocalCluster still gets created
+                if has_slurm and not spy.__acme__:
+                    slurm_msg = ("We are apparently on a slurm cluster but\n"
+                                 "Syncopy could not find a Dask client.\n"
+                                 "Syncopy does not provide an "
+                                 "automatic Dask SLURMCluster on its own!"
+                                 "\nPlease consider using ACME (https://github.com/esi-neuroscience/acme)"
+                                 "\nor configure your own cluster via `dask_jobqueue.SLURMCluster()`"
+                                 "\n\nCreating a LocalCluster as fallback.."
+                           )
+                    SPYWarning(slurm_msg)
+
+                # -- spawn fallback local cluster --
+
                 cluster = dd.LocalCluster()
                 # attaches to local cluster residing in global namespace
                 dd.Client(cluster)
@@ -518,6 +529,9 @@ def detect_parallel_client(func):
             dd.get_client().close()
             # and kill
             cluster.close()
+        # print again in case it got drowned
+        if slurm_msg:
+            SPYWarning(slurm_msg)
 
         return results
 
@@ -680,7 +694,7 @@ def process_io(func):
                 if keeptrials:
                     main_dset[outgrid] = res
                 else:
-                    main_dset[()] = np.nansum([main_dset, res], axis=0)
+                    main_dset[()] += res
 
                 h5_add_metadata(h5fout, details, unique_key_suffix=call_id)
                 h5fout.flush()
diff --git a/syncopy/shared/tools.py b/syncopy/shared/tools.py
index b72245f8..fa184152 100644
--- a/syncopy/shared/tools.py
+++ b/syncopy/shared/tools.py
@@ -272,22 +272,4 @@ def get_defaults(obj):
     return StructDict(dct)
 
 
-def check_slurm_available():
-    """
-    Returns `True` if a SLURM instance could be reached via
-    a `sinfo` call, `False` otherwise.
-    """
-
-    # Check if SLURM's `sinfo` can be accessed
-    proc = subprocess.Popen("sinfo",
-                            stdout=subprocess.PIPE, stderr=subprocess.PIPE,
-                            text=True, shell=True)
-    _, err = proc.communicate()
-    # Any non-zero return-code means SLURM is not available
-    # so we disable ACME
-    if proc.returncode != 0:
-        has_slurm = False
-    else:
-        has_slurm = True
-
-    return has_slurm
+    
diff --git a/syncopy/specest/compRoutines.py b/syncopy/specest/compRoutines.py
index 4536c11f..702b45b9 100644
--- a/syncopy/specest/compRoutines.py
+++ b/syncopy/specest/compRoutines.py
@@ -211,17 +211,19 @@ class MultiTaperFFT(ComputationalRoutine):
 
         check_freq_hashes(metadata, out)
 
-        # attached to `out` only here for testing of
-        # the whole machinery (see test_metadata.TestMetadataUsingMtmfft)
-        out.metadata = metadata
-
-        # channels, trialdefinition and so on..
+        # channels and trialdefinition
         propagate_properties(data, out, self.keeptrials)
 
-        if self.cfg["method_kwargs"]["taper"] is None:
+        taper_kw = self.cfg["method_kwargs"]["taper"]
+        if taper_kw is None:
             out.taper = np.array(['None'])
+        # multi-tapering
+        elif taper_kw == 'dpss':
+            nTaper =  self.outputShape[out.dimord.index("taper")]
+            out.taper = np.array([taper_kw + str(i) for i in range(nTaper)])
+        # just a single taper
         else:
-            out.taper = np.array([self.cfg["method_kwargs"]["taper"]] * self.outputShape[out.dimord.index("taper")])
+            out.taper = np.array([taper_kw])
         out.freq = self.cfg["foi"]
 
 
@@ -444,10 +446,18 @@ class MultiTaperFFTConvol(ComputationalRoutine):
         out.trialdefinition = trl
         out.samplerate = srate
         out.channel = np.array(data.channel[chanSec])
-        if self.cfg["method_kwargs"]["taper"] is None:
+
+        taper_kw = self.cfg["method_kwargs"]["taper"]
+        if taper_kw is None:
             out.taper = np.array(['None'])
+        # multi-tapering
+        elif taper_kw == 'dpss':
+            nTaper =  self.outputShape[out.dimord.index("taper")]
+            out.taper = np.array([taper_kw + str(i) for i in range(nTaper)])
+        # just a single taper
         else:
-            out.taper = np.array([self.cfg["method_kwargs"]["taper"]] * self.outputShape[out.dimord.index("taper")])
+            out.taper = np.array([taper_kw])
+
         out.freq = self.cfg["foi"]
 
 
@@ -620,6 +630,7 @@ class WaveletTransform(ComputationalRoutine):
         out.freq = 1 / self.cfg["method_kwargs"]["wavelet"].fourier_period(
             self.cfg["method_kwargs"]["scales"]
         )
+        out.taper = np.array(['None'])
 
 
 # -----------------
@@ -783,6 +794,7 @@ class SuperletTransform(ComputationalRoutine):
         out.channel = np.array(data.channel[chanSec])
         # for the SL Morlets the conversion is straightforward
         out.freq = 1 / (2 * np.pi * self.cfg["method_kwargs"]["scales"])
+        out.taper = np.array(['None'])
 
 
 def _make_trialdef(cfg, trialdefinition, samplerate):


Test patch:
diff --git a/syncopy/tests/test_connectivity.py b/syncopy/tests/test_connectivity.py
index fcfe36a3..09e71b74 100644
--- a/syncopy/tests/test_connectivity.py
+++ b/syncopy/tests/test_connectivity.py
@@ -27,6 +27,7 @@ availMem = psutil.virtual_memory().total
 minRAM = 5
 skip_low_mem = pytest.mark.skipif(availMem < minRAM * 1024**3, reason=f"less than {minRAM}GB RAM available")
 
+
 class TestSpectralInput:
     """
     When inputting SpectralData directly into connectivityanalysis, it has to fulfill
@@ -74,13 +75,14 @@ class TestSpectralInput:
         assert "expected AnalogData" in str(err.value)
 
     def test_tf_input(self):
-        """ No time-resolved spectral connectivity implemented yet """
+        """ No time-resolved Granger implemented yet """
         spec = spy.freqanalysis(self.ad, method='mtmconvol',
                                 t_ftimwin=0.01, output='fourier')
 
-        with pytest.raises(SPYValueError) as err:
-            cafunc(spec, method='coh')
-        assert "expected line spectra" in str(err.value)
+        with pytest.raises(NotImplementedError) as err:
+            cafunc(spec, method='granger')
+        assert "Granger causality from tf-spectra" in str(err.value)
+
 
 class TestGranger:
 
@@ -275,8 +277,8 @@ class TestGranger:
 
 class TestCoherence:
 
-    nSamples = 1500
-    nChannels = 6
+    nSamples = 1400
+    nChannels = 4
     nTrials = 100
     fs = 1000
 
@@ -309,6 +311,80 @@ class TestCoherence:
 
     spec = spy.freqanalysis(data, cfg, output='fourier', keeptapers=True)
 
+    def test_timedep_coh(self):
+        """
+        stack together different phase diffusing signals
+        with increasing diffusion strength to emulate non-stationarity
+        """
+
+        trls = []        
+        for _ in range(self.nTrials):
+            pds = []
+            for eps in [1e-3, 1e-2, 1e-1, 2e-1]:
+                pds.append(synth_data.phase_diffusion(freq=self.f2,
+                                                      eps=eps,
+                                                      nChannels=self.nChannels,
+                                                      nSamples=self.nSamples // 8,
+                                                      seed=None)
+                           + .5 * synth_data.white_noise(nChannels=self.nChannels,
+                                                           nSamples=self.nSamples // 8,
+                                                           seed=None)
+                           )
+
+            trial_dat = np.concatenate(pds)
+            trls.append(trial_dat)
+
+        test_data = spy.AnalogData(trls, samplerate=self.fs)
+        # check number of samples
+        nSamples = self.nSamples // 8 * 4
+        assert test_data.time[0].size == nSamples
+
+        # get time-frequency spec for the non-stationary signal
+        spec_tf = spy.freqanalysis(test_data, method='mtmconvol',
+                                   t_ftimwin=0.3, foilim=[5, 100],
+                                   output='fourier')
+
+        # check that we have still the same time axis
+        assert np.all(spec_tf.time[0] == test_data.time[0])
+
+        # abs spectra for plotting
+        spec_tf_abs = spy.SpectralData(data=np.abs(spec_tf.data),
+                                       samplerate=self.fs,
+                                       trialdefinition=spec_tf.trialdefinition)
+        spec_tf_abs.freq = spec_tf.freq
+        spec_tf_abs.multipanelplot(trials=0)
+
+        # rough check that the power indeed drops over time along the 40Hz band
+        profile_40 = spec_tf_abs.show(foi=self.f2, trials=2)[:,1]
+        assert profile_40.argmax() < profile_40.argmin()
+
+        # compute time dependent coherence
+        coh = cafunc(data=spec_tf, method='coh')
+
+        # check that we have still the same time axis
+        assert np.all(coh.time[0] == test_data.time[0])
+
+        # not exactly beautiful but it makes the point
+        coh.singlepanelplot(channel_i=0, channel_j=1)
+
+        # plot the coherence over time just along two different frequency bands        
+        ppl.figure()
+        cprofile40 = coh.show(foi=40, channel_i=0, channel_j=1)
+        ppl.plot(cprofile40, label='40Hz')
+        # here is nothing
+        cprofile10 = coh.show(foi=10, channel_i=0, channel_j=1)
+        ppl.plot(cprofile10, label='10Hz')
+        ppl.xlabel('samples')
+        ppl.ylabel('coherence')
+        ppl.legend()
+
+        # check that the 40 Hz band has high coherence only in the beginning
+        assert cprofile40.max() > 0.9
+        # later coherence goes down
+        assert cprofile40[int(0.8 * nSamples):].max() < 0.4
+        # side band never has high coherence
+        assert cprofile10.max() < 0.4
+
     def test_coh_solution(self, **kwargs):
 
         # re-run spectral analysis
diff --git a/syncopy/tests/test_metadata.py b/syncopy/tests/test_metadata.py
index 18fc1078..0eff04c5 100644
--- a/syncopy/tests/test_metadata.py
+++ b/syncopy/tests/test_metadata.py
@@ -441,143 +441,6 @@ class TestMetadataUsingFooof():
         plt.ion()
 
 
-class TestMetadataUsingMtmfft():
-    tfData = _get_fooof_signal()
-    expected_metadata_dict_entries = ["freqs_hash"]
-
-    @staticmethod
-    def get_mtmfft_cfg():
-        cfg = TestMetadataUsingFooof.get_fooof_cfg()
-        cfg.output = "pow"
-        return cfg
-
-    def test_sequential_mtmfft(self):
-        """
-        Test metadata propagation with mtmfft in sequential compute mode.
-
-        This also demonstrates howto use the hashed frequencies from the 2nd return
-        value to ensure they are identical across trials.
-        """
-        # These are known from the input data and cfg.
-        data_size = 100  # Number of samples (per trial) seen by mtmfftm. The full signal returned by _get_fooof_signal() is
-                         # larger, but the cfg.foilim setting (in get_fooof_cfg()) limits to 100 samples.
-        num_trials_out = 1 # Because of keeptrials = False in cfg.
-
-        cfg = TestMetadataUsingMtmfft.get_mtmfft_cfg()
-        cfg.parallel = False
-
-        spec_dt = freqanalysis(cfg, self.tfData)
-
-        assert spec_dt.data.shape == (num_trials_out, 1, data_size, 1)
-
-        keys_unique = [kv + "__0_0" for kv in self.expected_metadata_dict_entries]
-
-        assert spec_dt.metadata is not None
-        assert isinstance(spec_dt.metadata, dict)  # Make sure it is a standard dict, not a hdf5 group.
-        for kv in keys_unique:
-            assert kv in spec_dt.metadata.keys()
-            assert isinstance(spec_dt.metadata.get(kv), (np.bytes_))
-
-
-        # Demo: Use the extra return values to make sure the hashes of the frequency
-        #       arrays are identicals across all trials:
-        from hmac import compare_digest
-        reference_hash = None
-        trials_with_mismatches = []
-        num_hashes_checked = 0
-        for unique_md_label_rel, v in spec_dt.metadata.items():
-            label, trial_idx, _ = decode_unique_md_label(unique_md_label_rel)
-            if label == "freqs_hash":
-                trial_freqs_hash = v
-                if reference_hash is None:
-                    reference_hash = trial_freqs_hash
-                else:
-                    if not compare_digest(reference_hash, trial_freqs_hash):
-                        trials_with_mismatches.append(trial_idx)
-                num_hashes_checked += 1
-        if trials_with_mismatches:
-            SPYWarning(f"Frequency axes hashes mismatched for {len(trials_with_mismatches)} trials: {trials_with_mismatches} against reference hash from first trial.")
-        else:
-            print(f"Frequency axes hashes are identical across all {num_hashes_checked} trials.")
-
-    def test_sequential_mtmfft_with_selections(self):
-        """
-        Test metadata propagation with mtmfft in sequential compute mode, with selections.
-
-        This also demonstrates howto use the hashed frequencies from the 2nd return
-        value to ensure they are identical across trials.
-        """
-        # These are known from the input data and cfg.
-        data_size = 100  # Number of samples (per trial) seen by mtmfftm. The full signal returned by _get_fooof_signal() is
-                         # larger, but the cfg.foilim setting (in get_fooof_cfg()) limits to 100 samples.
-
-
-        cfg = TestMetadataUsingMtmfft.get_mtmfft_cfg()
-        cfg.parallel = False
-        cfg.keeptrials = True
-        selected_trials = [3, 5, 7]
-        cfg.select = { 'trials': selected_trials }
-        num_trials_out = len(selected_trials)
-
-        spec_dt = freqanalysis(cfg, self.tfData)
-
-        assert spec_dt.data.shape == (num_trials_out, 1, data_size, 1)
-
-        keys_unique = [kv + "__0_0" for kv in self.expected_metadata_dict_entries]
-
-        assert spec_dt.metadata is not None
-        assert isinstance(spec_dt.metadata, dict)  # Make sure it is a standard dict, not a hdf5 group.
-
-        # Check that the metadata keys are absolute, despite the selection.
-        md_trial_indices = []
-        for k, v in spec_dt.metadata.items():
-            label, trial_idx, chunk_idx = decode_unique_md_label(k)
-            md_trial_indices.append(int(trial_idx))
-            assert isinstance(v, (np.bytes_))
-        for ti in md_trial_indices:
-            assert ti in selected_trials, f"Expected trial index '{ti}' not found in selected_trials"
-
-
-
-    def test_par_mtmfft_with_selections(self):
-        """
-        Test metadata propagation with mtmfft in with parallel compute and parallel storage,
-        and trial selections.
-
-        In the case of trial selections, the computation of absolute trial indices
-        from relative ones uses the complex branch for mtmfft, where absolute trial
-        indices are computed from the relative ones.
-        """
-        # These are known from the input data and cfg.
-        data_size = 100  # Number of samples (per trial) seen by mtmfftm. The full signal returned by _get_fooof_signal() is
-                         # larger, but the cfg.foilim setting (in get_fooof_cfg()) limits to 100 samples.
-
-        cfg = TestMetadataUsingMtmfft.get_mtmfft_cfg()
-        cfg.parallel = True  # Enable parallel computation
-        cfg.keeptrials = True  # Enable parallel storage (is turned off when trial averaging is happening)
-
-        selected_trials = [3, 5, 7]
-        cfg.select = { 'trials': selected_trials }
-        num_trials_out = len(selected_trials)
-
-        spec_dt = freqanalysis(cfg, self.tfData)
-
-        assert spec_dt.data.shape == (num_trials_out, 1, data_size, 1)
-
-        assert spec_dt.metadata is not None
-        assert isinstance(spec_dt.metadata, dict)  # Make sure it is a standard dict, not a hdf5 group.
-
-        # Check that the metadata keys are absolute, despite the selection.
-        md_trial_indices = []
-        for k, v in spec_dt.metadata.items():
-            label, trial_idx, chunk_idx = decode_unique_md_label(k)
-            md_trial_indices.append(int(trial_idx))
-            assert isinstance(v, (np.bytes_))
-        for ti in md_trial_indices:
-            assert ti in selected_trials, f"Expected trial index '{ti}' not found in selected_trials"
-
-
 if __name__ == "__main__":
     T1 = TestMetadataHelpers()
     T2 = TestMetadataUsingFooof()
-    T3 = TestMetadataUsingMtmfft()
diff --git a/syncopy/tests/test_plotting.py b/syncopy/tests/test_plotting.py
index 6e610029..9f11ba83 100644
--- a/syncopy/tests/test_plotting.py
+++ b/syncopy/tests/test_plotting.py
@@ -156,6 +156,7 @@ class TestSpectralPlotting():
 
     spec_fft = spy.freqanalysis(adata, tapsmofrq=1)
     spec_fft_imag = spy.freqanalysis(adata, output='imag')
+    spec_fft_mtm = spy.freqanalysis(adata, tapsmofrq=1, keeptapers=True)
     spec_fft_complex = spy.freqanalysis(adata, output='fourier')
 
     spec_wlet = spy.freqanalysis(adata, method='wavelet',
@@ -179,6 +180,12 @@ class TestSpectralPlotting():
             fig1, ax1 = self.spec_fft.singlepanelplot(**kwargs)
             fig2, axs = self.spec_fft.multipanelplot(**kwargs)
 
+            # multi taper
+            assert self.spec_fft_mtm.taper.size > 1
+            _, _ = self.spec_fft_mtm.singlepanelplot(channel=3, **kwargs)
+            _, _ = self.spec_fft_mtm.singlepanelplot(taper=1, **kwargs)
+            _, _ = self.spec_fft_mtm.multipanelplot(**kwargs)
+
             _, _ = self.spec_fft_imag.singlepanelplot(**kwargs)
             _, _ = self.spec_fft_imag.multipanelplot(**kwargs)