如何使用 matplotlib 频谱图检测频率？

Question

我有以下压力信号的频谱图：

数据可用here。

我正在尝试计算压力波动的不同谐波的频率。我尝试计算图中白色水平虚线所示的其中一些。但所有频率代码拾取的值介于两者之间。我需要准确检测所有谐波的频率。我用来检测的代码：

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 >= dsmt2[i] and  dsmt2[i + 1] >= 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')

Answer 1

我希望这更接近您的需要

实在看不懂你是怎么计算谐波的，所以我就试着找上次信号已经稳定的时候的频谱值的峰值

sample=np.transpose(np.array(spectrum)[:,10])

ind=[]

wind=10

for i in xrange(7):
    maxind=np.argmax(sample) #find the index of the max value
    ind.append(maxind)       #save that index
    sample[max(0,maxind-(wind/2)):maxind+(wind/2)]=0    #clean around one window size chunk of data to suppress that peak

x=freqs[ind]


frequency_first_harmonic = x[1]  #drop x[0] because its going to be 0
frequency_second_harmonic = x[2]
frequency_third_harmonic = x[3]
frequency_fourth_harmonic = x[4]
frequency_fifth_harmonic = x[5]
frequency_sixth_harmonic = x[6]

如何使用 matplotlib 频谱图检测频率？

How to detect frequencies with matplotlib spectrogram?

python

matplotlib

spectrogram