Remove overpadding from FFT computation. Extra zeros were causing erroneous estimation of first and last PSD in timeseries

Author: jmcvey3

mhkit/dolfyn/binned.py

@@ -370,7 +370,6 @@ def _psd_base(
         noise=0,
         n_bin=None,
         n_fft=None,
-        n_pad=None,
         step=None,
     ):
         """
@@ -393,8 +392,6 @@ def _psd_base(
         n_fft : int
           n_fft of veldat2, number of elements per bin if 'None' is taken
           from VelBinner
-        n_pad : int (optional)
-          The number of values to pad with zero. Default = 0
         step : int (optional)
           Controls amount of overlap in fft. Default: the step size is
           chosen to maximize data use, minimize nens, and have a
@@ -413,11 +410,9 @@ def _psd_base(
         fs = self._parse_fs(fs)
         n_bin = self._parse_nbin(n_bin)
         n_fft = self._parse_nfft(n_fft)
-        if n_pad is None:
-            n_pad = min(n_bin - n_fft, n_fft)
         out = np.empty(self._outshape_fft(dat.shape, n_fft=n_fft, n_bin=n_bin))
         # The data is detrended in psd, so we don't need to do it here.
-        dat = self.reshape(dat, n_pad=n_pad)
+        dat = self.reshape(dat)
 
         for slc in slice1d_along_axis(dat.shape, -1):
             out[slc] = psd_1D(dat[slc], n_fft, fs, window=window, step=step)
@@ -474,8 +469,8 @@ def _csd_base(self, dat1, dat2, fs=None, window="hann", n_fft=None, n_bin=None):
         oshp[-2] = np.min([oshp[-2], int(dat2.shape[-1] // n_bin2)])
 
         # The data is detrended in psd, so we don't need to do it here:
-        dat1 = self.reshape(dat1, n_pad=n_fft)
-        dat2 = self.reshape(dat2, n_pad=n_fft)
+        dat1 = self.reshape(dat1)
+        dat2 = self.reshape(dat2)
         out = np.empty(oshp, dtype="c{}".format(dat1.dtype.itemsize * 2))
         if dat1.shape == dat2.shape:
             cross = cpsd_1D

mhkit/dolfyn/tools/fft.py

@@ -246,11 +246,13 @@ def cpsd_1D(a, b, nfft, fs, window="hann", step=None):
     window = _getwindow(window, nfft)
     fft_inds = slice(1, int(nfft / 2.0 + 1))
     wght = 2.0 / (window**2).sum()
+    # Take FFT of first segment
     s1 = fft(detrend_array(a[0:nfft]) * window)[fft_inds]
     if auto_psd:
         pwr = np.abs(s1) ** 2
     else:
         pwr = s1 * np.conj(fft(detrend_array(b[0:nfft]) * window)[fft_inds])
+    # take FFT of remaining segments
     if nens - 1:
         for i in range(step, l - nfft + 1, step):
             s1 = fft(detrend_array(a[i : (i + nfft)]) * window)[fft_inds]

mhkit/dolfyn/velocity.py

@@ -306,7 +306,7 @@ def u(
         """
         try:
             return self.ds["vel"][0].drop_vars("dir")
-        except:
+        except KeyError:
             return self.ds["vel_avg"][0].drop_vars("dir")
 
     @property
@@ -327,7 +327,7 @@ def v(
         """
         try:
             return self.ds["vel"][1].drop_vars("dir")
-        except:
+        except KeyError:
             return self.ds["vel_avg"][1].drop_vars("dir")
 
     @property
@@ -348,7 +348,7 @@ def w(
         """
         try:
             return self.ds["vel"][2].drop_vars("dir")
-        except:
+        except KeyError:
             return self.ds["vel_avg"][2].drop_vars("dir")
 
     @property
@@ -967,7 +967,6 @@ def power_spectral_density(
         noise=0,
         n_bin=None,
         n_fft=None,
-        n_pad=None,
         step=None,
     ):
         """
@@ -991,8 +990,6 @@ def power_spectral_density(
           The bin-size. Default = `self.n_bin`
         n_fft : int (optional)
           The fft size. Default = `self.n_fft`
-        n_pad : int (optional)
-          The number of values to pad with zero. Default = 0
         step : int (optional)
           Controls amount of overlap in fft. Default: the step size is
           chosen to maximize data use, minimize nens, and have a
@@ -1056,7 +1053,6 @@ def power_spectral_density(
                     noise=noise[idx],
                     window=window,
                     n_bin=n_bin,
-                    n_pad=n_pad,
                     n_fft=n_fft,
                     step=step,
                 )
@@ -1076,7 +1072,6 @@ def power_spectral_density(
                 noise=noise,
                 window=window,
                 n_bin=n_bin,
-                n_pad=n_pad,
                 n_fft=n_fft,
                 step=step,
             )