Inconsistency in function PSD when the NFFT parameter is an odd number

[cvr]

Function mlab.psd appears to have one inconsistency if the NFFT parameter is an odd number. In such case, as the highest frequency is not the Nyquist frequency, instead is int(.5*NFFT)/float(dt*NFFT), there are two distinct values showing both at the positive and negative frequencies in the FFT result. Consider the following example:

import numpy as np
from matplotlib import mlab
u = np.array([0, 1, 2, 3, 1, 2, 1])
dt = 1.0
Su = np.abs(np.fft.fft(u) * dt)**2 / float(dt * u.size)
P, f = mlab.psd(u, NFFT=u.size, Fs=1/dt, window=mlab.window_none,\
    detrend=mlab.detrend_none, noverlap=0, pad_to=None, scale_by_freq=None,\
    sides='onesided')

The resulting arrays Su and P are:

Su = [ 14.29, 1.61, 0.69, 0.55, 0.55, 0.69, 1.61 ]
P  = [ 14.29, 3.23, 1.39, 0.55 ]

To convert Su to one-sided output, it is common to sum the Fourier coefficients corresponding to negative and positive frequencies. Considering that u.size is odd, the frequencies which correspond to the coefficients in Su are:

df = 1 / float(u.size * dt)
f = np.array([ 0, 1, 2, 3, -3, -2, -1 ]) * df

thus, as the Su[1:4] are the positive side coefficients and Su[4:] the negative side, all coefficients are summed except for Su[0] (which in this case, as u is real is the same as Su[1:4] * 2):

Su_1side = np.append([Su[0]], Su[1:4] + Su[4:][::-1])

This results in Su_1side = [ 14.29, 3.23, 1.39, 1.10 ], whose last coefficient differs from the result of P, i.e. the mlab.psd result.

Note that the lack of scaling of the last coefficient is not a problem when the NFFT is an even number. In such case, as the as the frequency array is:

f = [ 0, df, 2*df, ..., (fnyquist - df), -fnyquist, -(fnyquist - df), ..., -2*df, -df ]

there is only one coefficient for the highest frequency (fnyquist), thus it is not summed up with another coefficient.

[tacaswell]

Thank you for the detailed bug report and explanation of what we were doing wrong!

Should be fixed by #4326

In [6]: P
Out[6]: array([ 14.28571429,   3.22739913,   1.38941047,   1.09747611])