title: How to detect frequencies with matplotlib spectrogram?
body: <p>I have the following spectrogram of a pressure signal:</p>
  
  <p><a href="https://i.stack.imgur.com/1JiMd.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/1JiMd.png" alt="enter image description here"></a></p>
  
  <p>The data is available <a href="https://drive.google.com/file/d/0BwwhEMUIYGyTczFWSk1RX1BTYVU/view?usp=sharing" rel="nofollow noreferrer">here</a>.</p>
  
  <p>I am trying to calculate frequency of the different harmonics of the pressure fluctuation. I made an attempt in calculating some of them shown by the white horizontal dotted lines in the figure. But all frequencies
  picked up by the code are in between. I need to detect the frequency of all the harmonics accurately. The code I used for detection:</p>
  
  <pre><code>from __future__ import division
  from matplotlib import ticker as mtick
  from matplotlib.backends.backend_pdf import PdfPages
  import matplotlib.pyplot as plt
  import numpy as np
  from bisect import bisect
  from scipy import interpolate
  from matplotlib.ticker import MaxNLocator
  import heapq
  import statsmodels.api as sm
  
  data = np.genfromtxt('spectrogram.dat', skiprows = 2, delimiter = ',')
  pressure = data[:, 1] * 0.065
  time = data[:, 0]
  time_start_ignition = 40
  
  
  with PdfPages('Spectorgram.pdf') as spectorgram_pressure:
      _spectorgram_pressure_vs_frequency_ = plt.figure(figsize=(5.15, 5.15))
      _spectorgram_pressure_vs_frequency_.clf()
      spectorgram_pressure_vs_frequency = plt.subplot(111)
      cax = plt.specgram(pressure * 100000, NFFT = 256, Fs = 50000, noverlap=4, cmap=plt.cm.gist_heat, zorder = 1)
      spectorgram_pressure_vs_frequency.grid(False, which="major")
      spectorgram_pressure_vs_frequency.set_xlabel('Time (s)', labelpad=6)
      spectorgram_pressure_vs_frequency.set_ylabel('Frequency (Hz)', labelpad=6)
      spectrum, freqs, t, im = cax
      dB = 10*np.log10(spectrum)
      f = interpolate.RectBivariateSpline(t, freqs,  dB.T) # but this uses xy not ij, hence the .T
      xnew = np.linspace(t[0], t[-1], 100*len(t))
      ynew = np.linspace(freqs[0], freqs[-1], 10*len(freqs)) # was it wider spaced than freqs on purpose?
      znew = f(xnew, ynew).T
      k = len(ynew) / 8
      smt = sm.nonparametric.lowess(znew[:, bisect(xnew, (time_start_ignition * max(xnew))/(max(time)))], np.linspace(0, len(ynew), len(ynew)), frac = 0.2)[:,1]
      i_change = []
      dsmt2 = np.diff(np.diff(smt))
      print dsmt2[10]
      for i in range(len(smt) - 3):
          if 0 &gt;= dsmt2[i] and  dsmt2[i + 1] &gt;= 0:
              i_change = np.append(i_change, i)
      i_change = np.int_(i_change)
      frequency_first_harmonic = ynew[np.argwhere(znew == np.nanmax(znew[i_change[0]:i_change[1], bisect(xnew, (time_start_ignition * max(xnew))/(max(time)))]))[0][0]]
      frequency_second_harmonic = ynew[np.argwhere(znew == np.nanmax(znew[i_change[1]:i_change[2], bisect(xnew, (time_start_ignition * max(xnew))/(max(time)))]))[0][0]]
      frequency_third_harmonic = ynew[np.argwhere(znew == np.nanmax(znew[i_change[2]:i_change[3], bisect(xnew, (time_start_ignition * max(xnew))/(max(time)))]))[0][0]]
      frequency_fourth_harmonic = ynew[np.argwhere(znew == np.nanmax(znew[i_change[3]:i_change[4], bisect(xnew, (time_start_ignition * max(xnew))/(max(time)))]))[0][0]]
      frequency_fifth_harmonic = ynew[np.argwhere(znew == np.nanmax(znew[i_change[4]:i_change[5], bisect(xnew, (time_start_ignition * max(xnew))/(max(time)))]))[0][0]]
      frequency_sixth_harmonic = ynew[np.argwhere(znew == np.nanmax(znew[i_change[5]:i_change[6], bisect(xnew, (time_start_ignition * max(xnew))/(max(time)))]))[0][0]]
  
      spectorgram_pressure_vs_frequency.axhline(frequency_first_harmonic, color = 'white', linewidth = 0.75, linestyle = ':')
      spectorgram_pressure_vs_frequency.axhline(frequency_second_harmonic, color = 'white', linewidth = 0.75, linestyle = ':')
      spectorgram_pressure_vs_frequency.axhline(frequency_third_harmonic, color = 'white', linewidth = 0.75, linestyle = ':')
      # spectorgram_pressure_vs_frequency.axhline(frequency_fourth_harmonic, color = 'white', linewidth = 0.75, linestyle = ':')
      # spectorgram_pressure_vs_frequency.axhline(frequency_fifth_harmonic, color = 'white', linewidth = 0.75, linestyle = ':')
      # spectorgram_pressure_vs_frequency.axhline(frequency_sixth_harmonic, color = 'white', linewidth = 0.75, linestyle = ':')   
      spectorgram_pressure_vs_frequency.axvline((time_start_ignition * max(cax[2]))/(max(time)), color = 'white', linewidth = 0.75, linestyle = ':')
      y_min, y_max = spectorgram_pressure_vs_frequency.get_ylim()
      cbar = plt.colorbar(orientation='vertical', ax = spectorgram_pressure_vs_frequency, fraction = 0.046, pad = 0.2)
      cbar.set_label('Power spectral density (dB)', rotation=90)
      primary_ticks = len(spectorgram_pressure_vs_frequency.yaxis.get_major_ticks())
      pressure_vs_time = spectorgram_pressure_vs_frequency.twinx()
      pressure_vs_time.grid(False)
      pressure_vs_time.plot((time * max(cax[2]))/(max(time)), pressure, linewidth = 0.75, linestyle = '-', color = 'k', alpha = 0.7, zorder = 3)
      pressure_vs_time.set_ylabel('Cylinder pressure (bar)', labelpad=6)
      pressure_vs_time.yaxis.set_major_locator(mtick.LinearLocator(primary_ticks))
      spectorgram_pressure_vs_frequency.set_xlim([0, max(cax[2])])
      spectorgram_pressure.savefig(bbox_inches='tight')
  </code></pre>