1

I am trying to understand the frequency bins returned by the matplotlib.mlab.psd() function.

Using the following code I can inspect the frequencies which are returned and I'm not convinced they are correct. 

import matplotlib.mlab as ml
import numpy as np
sampf=500.
nfft=2**4
testdat=np.random.randn(10000,)
p2,f2=ml.psd(testdat, nfft,sampf,sides='twosided')
p1,f1=ml.psd(testdat, nfft,sampf,sides='onesided')

print testdat.shape
print "Twosided"
print "\tbin1     : {:f} ".format(f2[0])
print "\tbin2     : {:f} ".format(f2[1])
print "\tbinlast  : {:f} ".format(f2[-1])

print "onesided"
print "\tbin1     : {:f} ".format(f1[0])
print "\tbin2     : {:f} ".format(f1[1])
print "\tbinlast  : {:f} ".format(f1[-1])

print "recreate"
f3=np.arange(nfft)*(sampf/2.)/nfft
print "\tbin1     : {:f} ".format(f3[0])
print "\tbin2     : {:f} ".format(f3[1])
print "\tbinlast  : {:f} ".format(f3[-1])

which gives this output:

Twosided
    bin1     : -250.000000 
    bin2     : -218.750000 
    binlast  : 218.750000 
onesided
    bin1     : 0.000000 
    bin2     : 31.250000 
    binlast  : 250.000000 
recreate
    bin1     : 0.000000 
    bin2     : 15.625000 
    binlast  : 234.375000

Am I right in thinking that the maximum frequency (binlast) for the 2 sided case should be half the sampling frequency?

Following this SO post I think it should range to sampf/2.

Community
  • 1
  • 1
mor22
  • 1,372
  • 1
  • 15
  • 32

1 Answers1

2

All the one sided does is not return the negative side.

Because you are handing in a real signal f_hat(w) = conj(f_hat(-w)) (that is the Fourier component at negative omega is the complex conjugate of the component at omega) thus they will have the same magnitude and thus in terms of the power spectrum are redundant.

If you are missing the exactly sampf/2 it is because of off-by-one issues related to having an even number of steps but needing and odd number of points if you are going to include 0 and be perfectly symmetric. Note that in your two sided case, the most negative frequency is -sampf/2 and your maximum misses sampf/2 by one bin step. Your reconstruction bin last is (nfft-1)/nfft * (sampf/2) and misses the value due to I suspect float rounding error.

tacaswell
  • 84,579
  • 22
  • 210
  • 199
  • Yes, I understand that the data which I am passing is Real and therefore I only need one half of the FFT. My question is what values should the frequency bins be labeled with? If I do a double sided fft then should the maximum bin frequency be (sampf/2) or should it be ((sampf/2) - (2/nfft)) ? – mor22 Jul 12 '13 at 12:47
  • Read the last paragraph of my answer again and then read https://en.wikipedia.org/wiki/Fast_Fourier_transform . – tacaswell Jul 12 '13 at 14:11
  • Ok so I have recalculated the values with a lower nfft to show that there is no rounding error. The wiki page (and others) indicate that the last bin of the single sided fft should be ((sampf/2) - (1/nfft)). Your point regarding the symmetric nature for a double sided fft of real data therefore indicates that the first bin should be labelled -((sampf/2) - (1/nfft)). This is not what is returned by the matplotlib function, which is want I'm interested in. I want to know whether I can trust the stock function. – mor22 Jul 12 '13 at 15:44
  • You have an _even_ number of bins (if you are using a power of 2 for your fft length). You need to include 0 -> you *must* have one more frequency bin on either the positive or negative side. This is in the structure of the DFFT. I assure you that the numpy implementation of the FFT is 1) correct 2) bug free 3) more efficient than anything you will write – tacaswell Jul 12 '13 at 16:10