#Создаем тестовый сигнал
import nump as np
Fs=1000   #частота выборки
t = np.arange(0, 1.0, 1.0 / Fs)
f1,f2,f3 = 100,200,300
signal = np.piecewise(t, [t < 1, t < 0.8, t < 0.3],
[lambda t: np.sin(2 * np.pi * f1 * t), lambda t: np.sin(2 * np.pi * f2 * t),
lambda t: np.sin(2 * np.pi * f3 * t)])                       #моделируемый сигнал1

# -*- coding: utf-8 -*-
""" Created on Fri Dec 17 21:18:48 2021 @author: lw """
import matplotlib.pyplot as plt
import numpy as np
from PyEMD import EEMD,EMD,Visualisation
from scipy.signal import hilbert
# from fftlw import fftlw
from vmdpy import VMD
#Метод декомпозиции (emd, eemd, vmd)
def decompose_lw(signal,t,method='eemd',K=5,draw=1):
names=['emd','eemd','vmd']
idx=names.index(method)
#emd разложение
if idx==0:
emd = EMD()
IMFs= emd.emd(signal)
#vmddecomposition
elif idx==2:
alpha = 2000       # moderate bandwidth constraint
tau = 0.            # noise-tolerance (no strict fidelity enforcement)
DC = 0             # no DC part imposed
init = 1           # initialize omegas uniformly
tol = 1e-7
# Run actual VMD code
IMFs, _, _ = VMD(signal, alpha, tau, K, DC, init, tol)
#eeemd разложение
else:
eemd = EEMD()
emd = eemd.EMD
emd.extrema_detection="parabol"
IMFs= eemd.eemd(signal,t)
#визуализация
if draw==1:
plt.figure()
for i in range(len(IMFs)):
plt.subplot(len(IMFs),1,i+1)
plt.plot(t,IMFs[i])
if i==0:
plt.rcParams['font.sans-serif']='Times New Roman'
plt.title('Decomposition Signal',fontsize=14)
elif i==len(IMFs)-1:
plt.rcParams['font.sans-serif']='Times New Roman'
plt.xlabel('Time/s')
# plt.tight_layout()
return IMFs       
#Преобразование Гильберта и рисование временного спектра
def hhtlw(IMFs,t,f_range=[0,500],t_range=[0,1],ft_size=[128,128],draw=1):
fmin,fmax=f_range[0],f_range[1]         #Диапазон частот, отображаемый частотно-временной диаграммой
tmin,tmax=t_range[0],t_range[1]         #диапазон времени
fdim,tdim=ft_size[0],ft_size[1]         #Размер (разрешение) частотно-временного графика
dt=(tmax-tmin)/(tdim-1)
df=(fmax-fmin)/(fdim-1)
vis = Visualisation()
#hilbertchange
c_matrix=np.zeros((fdim,tdim))
for imf in IMFs:
imf=np.array([imf])
#Находим мгновенную частоту
freqs = abs(vis._calc_inst_freq(imf, t, order=False, alpha=None))
#Находим мгновенную амплитуду
amp= abs(hilbert(imf))
#Удалить размерность 1
freqs=np.squeeze(freqs)
amp=np.squeeze(amp)
#Преобразовать в матрицу
temp_matrix=np.zeros((fdim,tdim))
n_matrix=np.zeros((fdim,tdim))
for i,j,k in zip(t,freqs,amp):
if i>=tmin and i<=tmax and j>=fmin and j<=fmax:
temp_matrix[round((j-fmin)/df)][round((i-tmin)/dt)]+=k
n_matrix[round((j-fmin)/df)][round((i-tmin)/dt)]+=1
n_matrix=n_matrix.reshape(-1)
idx=np.where(n_matrix==0)[0]
n_matrix[idx]=1
n_matrix=n_matrix.reshape(fdim,tdim)
temp_matrix=temp_matrix/n_matrix
c_matrix+=temp_matrix
t=np.linspace(tmin,tmax,tdim)
f=np.linspace(fmin,fmax,fdim)
#визуализация
if draw==1:
fig,axes=plt.subplots()
plt.rcParams['font.sans-serif']='Times New Roman'
plt.contourf(t, f, c_matrix,cmap="jet")
plt.xlabel('Time/s',fontsize=16)
plt.ylabel('Frequency/Hz',fontsize=16)
plt.title('Hilbert spectrum',fontsize=20)
x_labels=axes.get_xticklabels()
[label.set_fontname('Times New Roman') for label in x_labels]
y_labels=axes.get_yticklabels()
[label.set_fontname('Times New Roman') for label in y_labels]
# plt.show()
return t,f,c_matrix
#%%Тестовая функция
if __name__=='__main__':
#Создаем тестовый сигнал
Fs=1000   #частота выборки
t = np.arange(0, 1.0, 1.0 / Fs)
f1,f2,f3 = 100,200,300
signal = np.piecewise(t, [t < 1, t < 0.8, t < 0.3],
[lambda t: np.sin(2 * np.pi * f1 * t), lambda t: np.sin(2 * np.pi * f2 * t),
lambda t: np.sin(2 * np.pi * f3 * t)])                       #моделируемый сигнал1
# signal = 3*np.sin(2*np.pi*f1*t)+6*np.sin(2*np.pi*f2*t)+5*np.sin(2*np.pi*f3*t) #моделируемый сигнал2
# signal = 3*t*np.sin(2*np.pi*f1*t) #моделируемый сигнал3
#рисовать диаграмму во временной области
plt.figure()
plt.plot(t,signal)
plt.rcParams['font.sans-serif']='Times New Roman'
plt.xlabel('Time/s',fontsize=16)
plt.title('Original Signal',fontsize=20) plt.show()
#рисованиемоделируемый сигнальная спектрограмма
# _,_=fftlw(Fs,signal,1)
IMFs=decompose_lw(signal,t,method='vmd',K=10)                                    #сигнал разложения
tt,ff,c_matrix=hhtlw(IMFs,t,f_range=[0,500],t_range=[0,1],ft_size=[128,128])     #drawhilbertspectrum

# -*- coding: utf-8 -*-
""" Created on Fri Dec 17 20:17:42 2021 @author: lw """
import matplotlib.pyplot as plt
import numpy as np
import pywt
# from matplotlib.font_manager import FontProperties
# chinese_font = FontProperties(fname='/usr/share/fonts/truetype/wqy/wqy-microhei.ttc')
sampling_rate = 1000
t = np.arange(0, 1.0, 1.0 / sampling_rate)
f1 = 100
f2 = 200
f3 = 300
data = np.piecewise(t, [t < 1, t < 0.8, t < 0.3],
[lambda t: np.sin(2 * np.pi * f1 * t), lambda t: np.sin(2 * np.pi * f2 * t),
lambda t: np.sin(2 * np.pi * f3 * t)])
wavename = 'cgau8'
totalscal = 256
fc = pywt.central_frequency(wavename)
cparam = 2 * fc * totalscal
scales = cparam / np.arange(totalscal, 1, -1)
[cwtmatr, frequencies] = pywt.cwt(data, scales, wavename, 1.0 / sampling_rate)
plt.figure(figsize=(8, 8))
plt.subplot(211)
plt.plot(t, data)
# plt.xlabel(u"Время (секунды)", fontproperties=chinese_font)
# plt.title(u"Сегментированные сигналы и временной спектр 300 Гц, 200 Гц и 100 Гц", fontproperties=chinese_font, fontsize=20)
plt.subplot(212)
plt.contourf(t, frequencies, abs(cwtmatr))
# plt.ylabel(u"частота (Гц)", fontproperties=chinese_font)
# plt.xlabel(u"Время (секунды)", fontproperties=chinese_font)
plt.subplots_adjust(hspace=0.4)
plt.tight_layout()
plt.show()