Merge branch '233-job-failed-965-testmtmfft' into 'dev'

FIX: Yet another round of specest-related bugfixes

See merge request it/syncopy!88

syncopy/examples/ex_specest.py

@@ -43,7 +43,7 @@
     cfg.tapsmofrq = 20
     cfg.output = 'pow'
     cfg.keeptrials = False
-    cfg.keeptapers = False
+    cfg.keeptapers = True
     overallSpectrum = spy.freqanalysis(cfg, data)
     sys.exit()
 

syncopy/specest/freqanalysis.py

@@ -4,7 +4,7 @@
 # 
 # Created: 2019-01-22 09:07:47
 # Last modified by: Stefan Fuertinger [[email protected]]
-# Last modification time: <2019-10-08 16:49:32>
+# Last modification time: <2019-10-09 12:09:31>
 
 # Builtin/3rd party package imports
 import numpy as np
@@ -82,7 +82,9 @@ def freqanalysis(data, method='mtmfft', output='fourier',
         with respect to the longest trial found in `data`. For instance, 
         `pad = 'nextpow2'` pads all trials in `data` to the next power of 2 higher 
         than the sample-count of the longest trial in `data`. See :func:`syncopy.padding` 
-        for more information. If `pad` is `None`, no padding is performed. 
+        for more information. If `pad` is `None`, no padding is performed and
+        all trials have to have approximately the same length (up to next even 
+        sample-count). 
     padtype : str
         Values to be used for padding. Can be 'zero', 'nan', 'mean', 
         'localmean', 'edge' or 'mirror'. See :func:`syncopy.padding` for 
@@ -184,7 +186,6 @@ def freqanalysis(data, method='mtmfft', output='fourier',
         trialList = list(range(len(data.trials)))
         sinfo = data.sampleinfo
     lenTrials = np.diff(sinfo)
-    minSampleNum = lenTrials.min()
 
     # Ensure padding selection makes sense: do not pad on a by-trial basis but 
     # use the longest trial as reference acn compute `padlength` from there
@@ -205,6 +206,13 @@ def freqanalysis(data, method='mtmfft', output='fourier',
         minSamplePos = lenTrials.argmin()
         minSampleNum = padding(data._preview_trial(trialList[minSamplePos]), padtype, pad=pad,
                                padlength=padlength, prepadlength=True).shape[timeAxis]
+
+    else:
+        if np.unique((np.floor(lenTrials / 2))).size > 1:
+            lgl = "trials of approximately equal length"
+            act = "trials of unequal length"
+            raise SPYValueError(legal=lgl, varname="data", actual=act)
+        minSampleNum = lenTrials.min()
 
     # Construct array of maximally attainable frequencies
     minTrialLength = minSampleNum/data.samplerate
@@ -304,8 +312,8 @@ def freqanalysis(data, method='mtmfft', output='fourier',
         log_dct["nTaper"] = nTaper
 
         # Set up compute-kernel
-        specestMethod = MultiTaperFFT(nTaper, 
-                                      timeAxis, 
+        specestMethod = MultiTaperFFT(nTaper=nTaper, 
+                                      timeAxis=timeAxis, 
                                       taper=taper, 
                                       taperopt=taperopt,
                                       tapsmofrq=tapsmofrq,

syncopy/specest/mtmfft.py

@@ -4,7 +4,7 @@
 # 
 # Created: 2019-09-02 14:25:34
 # Last modified by: Stefan Fuertinger [[email protected]]
-# Last modification time: <2019-10-08 16:50:35>
+# Last modification time: <2019-10-09 12:12:19>
 
 # Builtin/3rd party package imports
 import numpy as np
@@ -19,7 +19,7 @@
 
 # Local workhorse that performs the computational heavy lifting
 @unwrap_io
-def mtmfft(trl_dat, nTaper, timeAxis,
+def mtmfft(trl_dat, nTaper=1, timeAxis=0,
            taper=spwin.hann, taperopt={}, tapsmofrq=None,
            pad="nextpow2", padtype="zero", padlength=None, foi=None,
            keeptapers=True, polyorder=None, output_fmt="pow",
@@ -83,7 +83,7 @@ def mtmfft(trl_dat, nTaper, timeAxis,
         dat = padding(dat, padtype, pad=pad, padlength=padlength, prepadlength=True)
     nSamples = dat.shape[0]
     nChannels = dat.shape[1]
-
+    
     # Determine frequency band and shape of output (time=1 x taper x freq x channel)
     nFreq = int(np.floor(nSamples / 2) + 1)
     fidx = slice(None)
@@ -102,13 +102,8 @@ def mtmfft(trl_dat, nTaper, timeAxis,
     if noCompute:
         return outShape, freq.spectralDTypes[output_fmt]
 
-    # Get final output shape from `chunkShape` keyword modulo per-worker channel-count
     # In case tapers aren't kept allocate `spec` "too big" and average afterwards
-    shp = list(chunkShape)
-    shp[-1] = nChannels
-    shp[1] = nTaper
-    chunkShape = tuple(shp)
-    spec = np.full(chunkShape, np.nan, dtype=freq.spectralDTypes[output_fmt])
+    spec = np.full((1, nTaper, nFreq, nChannels), np.nan, dtype=freq.spectralDTypes[output_fmt])
     fill_idx = tuple([slice(None, dim) for dim in outShape[2:]])
 
     # Actual computation

syncopy/tests/test_specest.py

@@ -36,7 +36,11 @@ class TestMTMFFT():
     nTrials = 8
     fs = 1024
     fband = np.linspace(1, fs/2, int(np.floor(fs/2)))
-    freqs = np.random.choice(fband[:-2], size=nChannels, replace=False)
+    freqs = [88.,  35., 278., 104., 405., 314., 271., 441., 343., 374., 428., 
+             367., 75., 118., 289., 310., 510., 102., 123., 417., 273., 449., 
+             416.,  32., 438., 111., 140., 304., 327., 494.,  23., 493.]
+    freqs = freqs[:nChannels]
+    # freqs = np.random.choice(fband[:-2], size=nChannels, replace=False)
     amp = np.pi
     phases = np.random.permutation(np.linspace(0, 2 * np.pi, nChannels))
     t = np.linspace(0, nTrials, nTrials * fs)
@@ -157,6 +161,11 @@ def test_dpss(self):
         cfg.method = "mtmfft"
         cfg.taper = "dpss"
         cfg.tapsmofrq = 9.3
+
+        # trigger error for non-equidistant trials w/o padding       
+        cfg.pad = None 
+        with pytest.raises(SPYValueError):
+            spec = freqanalysis(cfg, artdata)
 
         for sk, select in enumerate(self.artdataSelections):