1024 Python - 锐客网

在1024这个时间点，心中虽然感慨万千，却写不出那万水千山。还是老老实实的写个HHT吧——俗称希尔伯特变换，这里不涉及复杂的理论知识以及数学公式，这里只有简简单单的代码，如果想看详细的原理和公式推导的话请出门左拐。 HHT是在EMD的基础上进行的，首先把信号进行EMD分解（俗称的经验模态分解），然后根据EMD分解得到的各个分量进行Hilbert转换得到最终的时频图。
这里主要用的库就是Scipy这个库太强大了，这个库主要是用于科学计算，Scipy值得各位去研究一番。
友情提示：本次代码是在Jupyter上进行编码运行的
直接上代码了

from pyhht.visualization import plot_imfs import numpy as np from pyhht import EMD from scipy import signal from scipy.fftpack import fft,ifft import matplotlib.pyplot as plt import os import math import pandas as pd from scipy.signal import hilbert, chirp# -*- coding: utf-8 -*- #查询一分钟数据中的几秒数据 def getOneMinutesSecondData(arr,startSecond,endSecond):if endSecond == 0: endSecond = 60 beginIndex = startSecond*10240 print('begin: '+str(beginIndex)) endIndex = beginIndex+(endSecond-startSecond)*10240 print('end: '+str(endIndex)) return arr[beginIndex:endIndex]#绘制EMD分解后的图形 def plotIMFS(rawSingle): decomposer = EMD(rawSingle) imfs = decomposer.decompose() sizeEMD = len(imfs) totalImg = sizeEMD +1 t = np.linspace(0,1,len(rawSingle))plt.subplot(totalImg,1,1) plt.plot(t,rawSingle)for i in range(0,sizeEMD): index = i+2 plt.subplot(totalImg,1,index) plt.plot(t,imfs[i]) return imfsdef plotFFT(data): totalLentt = len(data) yf1tt=abs(np.fft.rfft(data))/totalLentt sizett = int(totalLentt/2)yf2tt = yf1tt[range(sizett)]xftt = np.linspace(0,sizett,sizett)plt.subplot(2,1,1) plt.plot(xftt[1:sizett],20*np.log(yf2tt[1:sizett])) plt.title('FFT db',fontsize=12)plt.subplot(2,1,2) plt.plot(xftt[1:sizett],yf2tt[1:sizett]) plt.title('',fontsize=12) return yf2ttdef plotEachHHT(imfsData):result = [] for i in range (2,len(imfsData)):if i==6: break analytic_signal = hilbert(imfsData[i])t = np.linspace(0,1,len(imfsData[i]))amplitude_envelope = np.abs(analytic_signal) instantaneous_phase = np.unwrap(np.angle(analytic_signal)) instantaneous_frequency = (np.diff(instantaneous_phase) / (2.0*np.pi) * 1024) result.append(instantaneous_frequency) fig=plt.gcf() fig.set_size_inches(20.5, 10.5)plt.subplot(2,1,1) plt.plot(t, imfsData[i], label='signal') plt.plot(t, amplitude_envelope, label='envelope') plt.xlabel("time in seconds")plt.subplot(2,1,2) plt.plot(t[1:], abs(instantaneous_frequency),label=str(i)) plt.xlabel("time in seconds") return resultdef plotSingleHHT(imfsData):analytic_signal = hilbert(imfsData) t = np.linspace(0,1,len(imfsData))amplitude_envelope = np.abs(analytic_signal) instantaneous_phase = np.unwrap(np.angle(analytic_signal)) instantaneous_frequency = (np.diff(instantaneous_phase) / (2.0*np.pi) * 10240)fig=plt.gcf() fig.set_size_inches(20.5, 10.5)plt.subplot(2,1,1) plt.plot(t, imfsData, label='signal') plt.plot(t, amplitude_envelope, label='envelope') plt.xlabel("time in seconds")plt.subplot(2,1,2) plt.plot(t[1:], instantaneous_frequency) plt.xlabel("time in seconds")

以上代码为方法准备阶段，接下来进行读取数据绘制图形

#加载数据arrCar = np.loadtxt("/root/Rone/waveData/carData.txt")xCar = np.linspace(0,30,len(arrCar)) plt.subplot(2,1,1) plt.plot(xCar,arrCar) fig=plt.gcf() fig.set_size_inches(20.5, 10.5)

文章图片
这次主要是对上图红色框内的数据进行EMD and HHT，要不然数据太多了会导致EMD分解的分量太多！

dataA = getOneMinutesSecondData(arrCar,5,7)xDataA = np.linspace(0,30,len(dataA)) plt.subplot(2,1,1) plt.plot(xDataA,dataA) fig=plt.gcf() fig.set_size_inches(18.5, 10.5)

截取到红色框的数据并进行绘图