title: Power spectrum of real data with fftpack on log axis
body: <p>I already read many discussion about this topic (<a href="https://stackoverflow.com/questions/20658318/my-power-spectra-are-believable-a-comparison-between-lomb-scargle-and-fft-scip?rq=1">comparison between lomb-scargle and fft </a>, <a href="https://stackoverflow.com/questions/15382076/plotting-power-spectrum-in-python">Plotting power spectrum in python</a>, <a href="https://stackoverflow.com/questions/9456037/scipy-numpy-fft-frequency-analysis">Scipy/Numpy FFT Frequency Analysis</a>, and many others), but still can't manage it, so I need some tips.
  I have a list of photon events (detections vs time), the data are available <a href="https://clbin.com/jw5Ls" rel="nofollow noreferrer">here</a>. The columns are <code>time</code>, <code>counts</code> , <code>errors</code>, and counts in different energy bands (you can ignore them). I know the source has a periodicity around <code>8.9 days = 1.3*10^-6 Hz</code>.
  I would like to plot the Power spectrum density showing a peak at this frequency (on a log x-axis, possibly). It would also be nice if I can avoid the half part of the plot (symmetric). This is my code till now, not so far but still something:</p>
  
  <pre><code>import numpy as np
  from scipy.fftpack import fft, rfft, fftfreq
  import pylab as plt
  
  x,y = np.loadtxt('datafile.txt', usecols = (0,1), unpack=True)
  y = y - y.mean() # Removes the large value at the 0 frequency that we don't care about
  
  f_range = np.linspace(10**(-7), 10**(-5), 1000)
  W = fftfreq(y.size, d=x[1]-x[0])
  
  plt.subplot(2,1,1)
  plt.plot(x,y)
  plt.xlabel('Time (days)')
  
  f_signal = fft(y)
  plt.subplot(2,1,2)
  plt.plot(W, abs(f_signal))
  plt.xlabel('Frequency (Hz)')
  </code></pre>
  
  <p>Here the (useless) plot produced:
  <img src="https://i.stack.imgur.com/PUoJP.png" alt="DFT"></p>