[Bug]: Possible problem with parameter NFFT and pad_to of mlab._spectral_helper #24822

Activity

1. jklymak commented on Dec 27, 2022

   

   jklymak

   on Dec 27, 2022

   Member

   More actions

   This implementation is pretty clunky, and we strongly recommend just using scipy. Indeed, I was hoping to deprecate the whole thing, but there wasn't sufficient support (#22920)

   You are not quite using the parameters correctly, but I agree the results are still incorrect:

   import numpy as np
   from matplotlib import mlab
   import matplotlib.pyplot as plt
   
   fs = 1000
   t = np.arange(0, 1, 1/fs)
   rng = np.random.default_rng(42)
   s = np.sin(2 * np.pi * 100 * t) + 0.5 * rng.normal(size=t.shape)
   
   nfft = 2048
   nsignal = len(s) #1000
   psd, freqs = mlab.psd(s, pad_to=nfft, NFFT=nfft, Fs=fs, window=mlab.window_none,
                         scale_by_freq=True)
   psd2, freqs2 = mlab.psd(s, pad_to=nsignal, NFFT=nsignal, Fs=fs, window=mlab.window_none,
   scale_by_freq=True)
   
   fig, ax = plt.subplots()
   df = np.median(np.diff(freqs))
   ax.loglog(freqs, psd, label=f'NFFT = 2048: {np.sum(psd[1:] * df)}')
   
   df2 = np.median(np.diff(freqs2))
   ax.loglog(freqs2, psd2, label=f'NFFT=len(data): {np.sum(psd2[1:] * df2)}')
   ax.legend()
   ax.set_title(f'Signal variance: {np.var(s)}')
   plt.show()

   yields:

   

   So, the padded normalization is incorrect by a factor of 1000/2048.
2. gapplef commented on Dec 28, 2022

   

   gapplef

   on Dec 28, 2022

   Last edited by gapplef

   ContributorAuthor

   More actions

   If NFFT was meant to be the actual length of signal (segment) for DFT, I think the following code is the source of the bug.
   

   matplotlib/lib/matplotlib/mlab.py

   Lines 342 to 351 in 3418bad

   	# zero pad x and y up to NFFT if they are shorter than NFFT
   	if len(x) < NFFT:
   	n = len(x)
   	x = np.resize(x, NFFT)
   	x[n:] = 0
   	if not same_data and len(y) < NFFT:
   	n = len(y)
   	y = np.resize(y, NFFT)
   	y[n:] = 0

   And the problem should be fixed with:

   # Check if x and y are the same length, zero-pad if necessary
   if not same_data:
       if len(x) != len(y):
           if len(x) < len(y):
               n = len(x) 
               x = np.resize(x, len(y)) 
               x[n:] = 0
           else:
               n = len(y) 
               y = np.resize(y, len(x)) 
               y[n:] = 0
   
   if len(x) < NFFT: 
      NFFT = len(x) 

   ---

   Also the following code need be fixed for correct one-sided psd:
   

   matplotlib/lib/matplotlib/mlab.py

   Lines 411 to 416 in 3418bad

   	# if we have a even number of frequencies, don't scale NFFT/2
   	if not NFFT % 2:
   	slc = slice(1, -1, None)
   	# if we have an odd number, just don't scale DC
   	else:
   	slc = slice(1, None, None)

    # if we have a even number of frequencies, don't scale pad_to/2 
    if not pad_to % 2: 
        slc = slice(1, -1, None) 
    # if we have an odd number, just don't scale DC 
    else: 
        slc = slice(1, None, None) 
3. jklymak commented on Dec 28, 2022

   

   jklymak

   on Dec 28, 2022

   Member

   More actions

   @gapplef Maybe? _spectral_helper is written so that pad_to ends up as n to FFT. That is not what I would have done, and I don't understand why it is that way. However, as above, I think we should not be dragging around a different API from scipy.signal, and should just remove all this, the practical problem is that trying to fix the inconsistencies will introduce fixes that are not necessarily back-compatible.
4. QuLogic commented on Dec 28, 2022

   

   QuLogic

   on Dec 28, 2022

   Member

   More actions

   cc @oscargus
5. jklymak

   added

   topic: spectral

   on Jan 28, 2023
6. misjob commented on Jan 9, 2024

   

   misjob

   on Jan 9, 2024

   More actions

   Also the zero padding if NFFT is less then the signal length should be performed after the detrend. Otherwise the zero padding would screw up the "mean" and "linear" calculations.