tsa.periodogram() returns frequencies of all 0s when Fs=1.

[Xunius]

Hi,

I noticed that the following code returns freqs which is all 0s:

    freqs, d_psd = tsa.periodogram(ar_seq, Fs=1., normalize=False)

I believe it is this line (in algorithms/spectral.py) that is causing the issue:

freqs = np.linspace(0, Fs // 2, Fn)

Should it be Fs / 2 instead?

Version: nitime 0.9
Installed via conda

[Xunius]

Update:

I think we should use scipy.fftpack.fftfreq(N) instead. The current way of computing freqs is not quite right for odd numbered sequence lengths. See the following snippet:

import scipy.fftpack
import numpy as np

Fs = 1.

N = 10 # even case
Fn = N // 2 + 1  # = 5
Fl = (N + 1) // 2
freqs = np.linspace(0, Fs / 2, Fn)  # current way of computing freq
freqs2 = scipy.fftpack.fftfreq(N)   # my suggestion
print('\nEven case:', 'N=', N, 'Fn=', Fn, 'Fl=', Fl)
print('freqs:')
print(freqs)
print('freqs that get doubled by 1:Fl indexing:', freqs[1:Fl])
print()
print('freqs2:')
print(freqs2)
print('freqs2 that get doubled by 1:Fl indexing:', freqs2[1:Fl])

N = 9 # odd case
Fn = N // 2 + 1  # = 5
Fl = (N + 1) // 2
freqs = np.linspace(0, Fs / 2, Fn)  # current way of computing freq
freqs2 = scipy.fftpack.fftfreq(N)   # my suggestion
print('\nodd case:', 'N=', N, 'Fn=', Fn, 'Fl=', Fl)
print('freqs:')
print(freqs)
print('freqs that get doubled by 1:Fl indexing:', freqs[1:Fl])
print()
print('freqs2:')
print(freqs2)
print('freqs2 that get doubled by 1:Fl indexing:', freqs2[1:Fl])

The output:

Even case: N= 10 Fn= 6 Fl= 5
freqs:
[0.  0.1 0.2 0.3 0.4 0.5]
freqs that get doubled by 1:Fl indexing: [0.1 0.2 0.3 0.4]

freqs2:
[ 0.   0.1  0.2  0.3  0.4 -0.5 -0.4 -0.3 -0.2 -0.1]
freqs2 that get doubled by 1:Fl indexing: [0.1 0.2 0.3 0.4]

odd case: N= 9 Fn= 5 Fl= 5
freqs:
[0.    0.125 0.25  0.375 0.5  ]
freqs that get doubled by 1:Fl indexing: [0.125 0.25  0.375 0.5  ]

freqs2:
[ 0.          0.11111111  0.22222222  0.33333333  0.44444444 -0.44444444
 -0.33333333 -0.22222222 -0.11111111]
freqs2 that get doubled by 1:Fl indexing: [0.11111111 0.22222222 0.33333333 0.44444444]

Note that for the odd case, the power corresponding to the Nyquist frequency 0.5 gets doubled by this line:

        P[..., 1:Fl] = 2 * (Sk[..., 1:Fl] * Sk[..., 1:Fl].conj()).real

While the scipy.fftpack.fftfreq(N) solution doesn't.

Even if we ignore the doubling, this is consistent with the Fourier coefficients computed using scipy.fftpack so the coefficients match their corresponding frequencies.

[arokem]

Yeah - I would be fine with such a change. Do you want to put up a PR?

[arokem]

Closed through #194