CHG: Normalize spectra independent of padding

CHANGELOG.md

@@ -13,6 +13,7 @@ All notable changes to this project will be documented in this file.
 - concatenation of syncopy data objects along trials
 
 ### CHANGED
+- spectral power for `mtmfft` now independent of padding as originally intended
 - support unequal trial sizes for `load_ft_raw`
 - major performance improvements for DiscreteData #403 #418, #424
 

syncopy/specest/freqanalysis.py

@@ -56,7 +56,7 @@ def freqanalysis(data, method='mtmfft', output='pow',
                  taper_opt=None, tapsmofrq=None, nTaper=None, keeptapers=False,
                  toi="all", t_ftimwin=None, wavelet="Morlet", width=6, order=None,
                  order_max=None, order_min=1, c_1=3, adaptive=False,
-                 out=None, fooof_opt=None, **kwargs):
+                 out=None, fooof_opt=None, ft_compat=False, **kwargs):
     """
     Perform (time-)frequency analysis of Syncopy :class:`~syncopy.AnalogData` objects
 
@@ -179,7 +179,7 @@ def freqanalysis(data, method='mtmfft', output='pow',
         For the default `maxperlen`, no padding is performed in case of equal
         length trials, while trials of varying lengths are padded to match the
         longest trial. If `pad` is a number all trials are padded so that `pad` indicates
-        the absolute length of all trials after padding (in seconds). For instance
+        the absolute length of all trials after padding in seconds. For instance
         ``pad = 2`` pads all trials to an absolute length of 2000 samples, if and
         only if the longest trial contains at maximum 2000 samples and the
         samplerate is 1kHz. If `pad` is `'nextpow2'` all trials are padded to the
@@ -284,6 +284,10 @@ def freqanalysis(data, method='mtmfft', output='pow',
         `FOOOF docs <https://fooof-tools.github.io/fooof/generated/fooof.FOOOF.html#fooof.FOOOF>`_
         for the meanings and the defaults.
         See the FOOOF reference [Donoghue2020]_ for details.
+    ft_compat : bool, optional
+        Set to `True` to use Field Trip's spectral normalization for FFT based methods
+        (``method='mtmfft'`` and ``method='mtmconvol'``). So spectral power is NOT
+        independent of the padding size!
 
     Returns
     -------
@@ -387,7 +391,7 @@ def freqanalysis(data, method='mtmfft', output='pow',
         raise SPYValueError(legal=lgl, varname="output", actual=output)
 
     # Parse all Boolean keyword arguments
-    for vname in ["keeptrials", "keeptapers", "demean_taper"]:
+    for vname in ["keeptrials", "keeptapers", "demean_taper", "ft_compat"]:
         if not isinstance(lcls[vname], bool):
             raise SPYTypeError(lcls[vname], varname=vname, expected="Bool")
 
@@ -605,7 +609,8 @@ def freqanalysis(data, method='mtmfft', output='pow',
             'taper': taper,
             'taper_opt': taper_opt,
             'nSamples': minSampleNum,
-            'demean_taper': demean_taper
+            'demean_taper': demean_taper,
+            'ft_compat': ft_compat
         }
 
         # Set up compute-class

syncopy/specest/mtmconvol.py

@@ -15,7 +15,7 @@
 
 
 def mtmconvol(data_arr, samplerate, nperseg, noverlap=None, taper="hann",
-              taper_opt={}, boundary='zeros', padded=True, detrend=False):
+              taper_opt=None, boundary='zeros', padded=True, detrend=False):
 
     """
     (Multi-)tapered short time fast Fourier transform. Returns
@@ -37,7 +37,7 @@ def mtmconvol(data_arr, samplerate, nperseg, noverlap=None, taper="hann",
     taper : str or None
         Taper function to use, one of `scipy.signal.windows`
         Set to `None` for no tapering.
-    taper_opt : dict
+    taper_opt : dict or None
         Additional keyword arguments passed to the `taper` function.
         For multi-tapering with ``taper='dpss'`` set the keys
         `'Kmax'` and `'NW'`.
@@ -91,6 +91,9 @@ def mtmconvol(data_arr, samplerate, nperseg, noverlap=None, taper="hann",
 
     taper_func = getattr(signal.windows, taper)
 
+    if taper_opt is None:
+        taper_opt = {}
+
     # this parameter mitigates the sum-to-zero problem for the odd slepians
     # as signal.stft has hardcoded scaling='spectrum'
     # -> normalizes with win.sum() :/

syncopy/specest/mtmfft.py

@@ -18,7 +18,8 @@ def mtmfft(data_arr,
            nSamples=None,
            taper="hann",
            taper_opt=None,
-           demean_taper=False):
+           demean_taper=False,
+           ft_compat=False):
     """
     (Multi-)tapered fast Fourier transform. Returns
     full complex Fourier transform for each taper.
@@ -45,6 +46,9 @@ def mtmfft(data_arr,
         `SciPy docs <https://docs.scipy.org/doc/scipy/reference/signal.windows.html>`_
     demean_taper : bool
         Set to `True` to perform de-meaning after tapering
+    ft_compat : bool
+        Set to `True` to use Field Trip's normalization,
+        which is NOT independent of the padding size
 
     Returns
     -------
@@ -107,8 +111,13 @@ def mtmfft(data_arr,
         if demean_taper:
             win -= win.mean(axis=0)
         ftr[taperIdx] = np.fft.rfft(win, n=nSamples, axis=0)
-        # FT uses potentially padded length `nSamples`
-        ftr[taperIdx] = _norm_spec(ftr[taperIdx], nSamples, samplerate)
+        # FT uses potentially padded length `nSamples`, which dilutes the power
+        if ft_compat:            
+            ftr[taperIdx] = _norm_spec(ftr[taperIdx], nSamples, samplerate)
+        # here the normalization adapts such that padding is NOT changing power
+        else:
+            ftr[taperIdx] = _norm_spec(ftr[taperIdx], signal_length * np.sqrt(nSamples / signal_length), samplerate)
+
     return ftr, freqs
 
 

syncopy/tests/test_specest.py

@@ -103,7 +103,7 @@ class TestMTMFFT():
     fband = np.linspace(1, fs/2, int(np.floor(fs/2)))
     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.]
+             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
@@ -190,6 +190,87 @@ def test_solution(self):
             # ensure amplitude is consistent across all channels/trials
             assert np.all(np.diff(amps) < self.tols[sk])
 
+    def test_normalization(self):
+
+        nSamples = 1000
+        fsample = 500  # 2s long signal
+        Ampl = 4  # amplitude
+        # 50Hz harmonic, spectral power is given by: Ampl^2 / 2 = 8
+        signal = Ampl * np.cos(2 * np.pi * 50 * np.arange(nSamples) * 1 / fsample)
+
+        # single signal/channel is enough
+        ad = AnalogData([signal[:, None]], samplerate=fsample)
+
+        cfg = StructDict()
+        cfg.foilim = [40, 60]
+        cfg.output = 'pow'
+        cfg.taper = None
+
+        # -- syncopy's default, padding does NOT change power --
+
+        cfg.ft_compat = False
+        cfg.pad = 'maxperlen'  # that's the default -> no padding
+        spec = freqanalysis(ad, cfg)
+        peak_power = spec.show().max()
+        df_no_pad = np.diff(spec.freq)  # freq. resolution
+        assert np.allclose(peak_power, Ampl**2 / 2, atol=1e-5)
+
+        cfg.pad = 4  # in seconds, double the size
+        spec = freqanalysis(ad, cfg)
+        df_with_pad = np.diff(spec.freq)
+        # we double the spectral resolution
+        assert np.allclose(df_no_pad[0], 2 * df_with_pad[0])
+        # yet power stays the same
+        peak_power = spec.show().max()
+        assert np.allclose(peak_power, Ampl**2 / 2, atol=1e-5)
+
+        # -- FT compat mode, padding does dilute the power --
+
+        cfg.ft_compat = True
+        cfg.pad = 'maxperlen'  # that's the default
+        spec = freqanalysis(ad, cfg)
+        peak_power = spec.show().max()
+        df_no_pad = np.diff(spec.freq)
+        # default padding is no padding if all trials are equally sized,
+        # so here the results are the same
+        assert np.allclose(peak_power, Ampl**2 / 2, atol=1e-5)
+
+        cfg.pad = 4  # in seconds, double the size
+        spec = freqanalysis(ad, cfg)
+        df_with_pad = np.diff(spec.freq)
+        # we double the spectral resolution
+        assert np.allclose(df_no_pad[0], df_with_pad[0] * 2)
+        # here half the power is now lost!
+        peak_power = spec.show().max()
+        assert np.allclose(peak_power, Ampl**2 / 4, atol=1e-5)
+
+        # -- works the same with tapering --
+
+        cfg.ft_compat = False
+        cfg.pad = 'maxperlen'  # that's the default
+        cfg.taper = 'kaiser'
+        cfg.taper_opt = {'beta': 10}
+        spec = freqanalysis(ad, cfg)
+        peak_power_no_pad = spec.show().max()
+
+        cfg.pad = 4
+        spec = freqanalysis(ad, cfg)
+        peak_power_with_pad = spec.show().max()
+        assert np.allclose(peak_power_no_pad, peak_power_with_pad, atol=1e-5)
+
+        cfg.ft_compat = True
+        cfg.pad = 'maxperlen'  # that's the default
+        cfg.taper = 'kaiser'
+        cfg.taper_opt = {'beta': 10}
+        spec = freqanalysis(ad, cfg)
+        peak_power_no_pad = spec.show().max()
+
+        cfg.pad = 4
+        spec = freqanalysis(ad, cfg)
+        peak_power_with_pad = spec.show().max()
+        # again half the power is lost with FT compat
+        assert np.allclose(peak_power_no_pad, 2 * peak_power_with_pad, atol=1e-5)
+
     def test_foi(self):
         for select in self.sigdataSelections: