nipy
tsa.periodogram() returns frequencies of all 0s when Fs=1.
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
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.
