title: compute coherence in python
body: <p>I'm learning cross-spectrum and coherence. From what I understand, coherence is like the analogue of correlation in that you normalize the cross-spectrum by the product of individual power spectrum:</p>
  
  <p><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/9c64700cf2b13cd0dfa20b45c9a591a29c142d52" alt="wikipedia equation"></p>
  
  <p>Here is my current python implementation</p>
  
  <pre><code>import numpy
  
  def crossSpectrum(x,y):
  
      #-------------------Remove mean-------------------
      xp=x-numpy.mean(x)
      yp=y-numpy.mean(y)
      n=len(x)
  
      # Do FFT
      cfx=numpy.fft.fft(xp)/n
      cfy=numpy.fft.fft(yp)/n
      freq=numpy.fft.fftfreq(n)
  
      # Get cross spectrum
      cross=cfx.conj()*cfy
  
      return cross,freq
  
  
  #-------------Main---------------------------------
  if __name__=='__main__':
  
      x=numpy.linspace(-250,250,500)
      noise=numpy.random.random(len(x))
      y=10*numpy.sin(2*numpy.pi*x/10.)+5*numpy.sin(2*numpy.pi*x/5.)+\
              2*numpy.sin(2*numpy.pi*x/20.)+10
      y+=noise*10
      y2=5*numpy.sin(2*numpy.pi*x/10.)+5+noise*50
  
      p11,freq=crossSpectrum(y,y)
      p22,freq=crossSpectrum(y2,y2)
      p12,freq=crossSpectrum(y,y2)
  
      # coherence
      coh=numpy.abs(p12)**2/p11.real/p22.real
      print coh
  </code></pre>
  
  <p>And my computed coherence is an array of 1s. What am I doing wrong?</p>
  
  <p>Also, sometimes the coherence plot has downward pointing spikes (like the output of <code>scipy.signal.coherence</code>, in other places that are pointing upward (e.g. <a href="http://www2.ocean.washington.edu/flowmow/processing/currents/coh/billtest/cohtest.html" rel="nofollow noreferrer">here</a>). I'm bit confused by the interpretation of coherence, shouldn't a larger coherence implies covariability between the 2 timeseries at that frequency?</p>
  
  <p>Thanks in advance.</p>