Re: [matplotlib-devel] Comparison of mlab.csd and Matlab's cpsd

[Ludwig S. <lud...@gm...>]

Hi,

> From: Ariel Rokem <ar...@be...>
> However - two elements are off by a factor of approximately 2 - the
> very first element and the very last. ... Does anyone have any idea
> why this would be the case?

>From a quick look at the mlab code, it looks like a bug in
mlab._spectral_helper.

The default spectrum is 'onesided' (same as for Matlab's cpsd). A
single-sided spectrum of a real signal has double the magnitude of a
double-sided spectrum, *except* at the origin (frequency index n = 0)
and Nyquist frequency (n = NFFT / 2), where it is the *same* as the
double-sided one [1]_.

In the mlab code, all the spectral values are simply scaled by a
factor of 2 (among other factors) in this line:

    # Scale the spectrum by the norm of the window to compensate for
    # windowing loss; see Bendat & Piersol Sec 11.5.2.  Also include
    # scaling factors for one-sided densities and dividing by the sampling
    # frequency, if desired.
    Pxy *= scaling_factor / (np.abs(windowVals)**2).sum()

This should be easy to fix (although the function probably needs a
little rework).

Regards,
Ludwig

Quick reference from my bookshelf:
---------------------------------------------------
.. [1] W. L. Briggs, V. E. Henson, "The DFT: An Owner's Manual for the
Discrete Fourier Transform," Section 1.3, Problem 6 (a), p. 13.