Welch Periodogram in Wolfram Language

This is a continuation of my post on the wave spectrum in Wolfram Language. I discovered that my Wolfram Language code could be significantly improved given some limitations in the PeriodogramArray function and an error on my part:

First, PeriodogramArray does not allow for detrending each partition before Fourier transformation. I have to detrend the whole record before passing it as an argument to the function. Function welch in SciPy applies the specified detrending technique to each partition. I consider this to be a good practice, even if this is not the case for MATLAB's pwelch function. I liked Octave's approach even more: you can specify if you want “short” or “long” detrending, i.e., remove the trend from each partition or from the full record, respectively.

Second, PeriodogramArray computes a symmetric window, not periodic. According to both MATLAB and SciPy documentations, a symmetric window is used for filter design and a periodic window is used for spectral analysis. The record partitions are multiplied with the data window before Fourier transformation, and the resultant spectrum is divided by a normalization constant computed from the window values. Consequently, a difference in the window values, resulting from it being symmetric rather than periodic, can have a non-negligible impact on the overall result.

Finally, my implementation doubles all spectral densities to account for both sides of the spectrum, but, actually, some frequencies have unpaired power spectral density (PSD), i.e., they are symmetry points, and should not be multiplied by two. For instance, at frequency zero, the PSD is always unpaired because there is symmetry about zero. If the record length is even, the last discrete frequency equals the Nyquist frequency, a symmetry point, so it is not doubled. However, if the record length is odd the Nyquist frequency lies between two discrete frequency points, where one is the mirrored counterpart of the other. Thus, if half the spectrum is considered, the last frequency point would be the one to the left of the Nyquist frequency, so its PSD must be doubled.

computePeriodogram[list_, w_] := 
  With[{n = Length@list}, {m = n/2 + 1}, 
    list // Fourier[# w] & // #[[1 ;; m]] & // Abs[#]^2 &]

doublePeriodogram[list_] := 
  With[{m = Length@list}, 
    list //
      MapIndexed[If[First[#2] == 1 || First[#2] == m, #1, 2 #1] &]]

computeWindow[wfun_, n_] := 
  wfun /@ N[Subdivide[-1/2, 1/2 - 1/n, n - 1]]

welchPeriodogram[list_, p_, wfun_, dfun_, fs_] := 
  With[{n = 2^p}, {w = computeWindow[wfun, n]}, 
    {f = Range[0, n/2] fs/n}, 
    Partition[list, n, n/2] // Map[dfun] // 
          Map[computePeriodogram[#, w] &] // Mean // 
        doublePeriodogram // (# n/(fs Total[w^2])) & // 
      Transpose[{f, #}] &]

Function welchPeriodogram returns the periodogram of the record list sampled at frequency fs. This implementation uses partitions of size 2^p and data window defined by function wfun. Each partition is detrended with function dfun before transformation. The code above gives the same result as SciPy's welch function for the same input and parameters.

Function wfun can be any window function already implemented in Wolfram Language. For example, HannWindow. The function calling, however, is modified to get a periodic window instead of symmetric.

The application of dfun corresponds to a “short” or partition detrending. If the user wants to apply a “long” or full-record detrending, they can pass function Identity and perform detrending to the full record, list, before evaluating welchPeriodogram. The same idea applies to wfun.

An example of detrending function is the following. It works for both full record and partitions.

linearDetrend[list_] := 
 With[{t = Range@Length@list}, 
  list - (LinearModelFit[list, xx, xx] /@ t)]

Additionally, the code can be adapted to use partitions with length different to powers of 2, but, at least in my computer, function Fourier runs twice as fast for powers of 2. This is expected and documented in the language reference. Note that the last PSD is always doubled in my function because powers of 2 are even numbers; fix accordingly.

Interestingly, there is already a resource function called WelchSpectralEstimate in the Wolfam Function Repository. This function aims to replicate SciPy's welch function but it has a minor bug. The last PSD of the periodogram is doubled when the partition length is even instead of when it is odd. So, if you run WelchSpectralEstimate and welch with the same input and parameters, the periodograms are identical for all the points but the last. I already reported this minor bug.