一种新颖的改进经验模态分解方法-通过迭代方式(IMF振幅加权频率)有效缓解了模态混叠缺陷,以后慢慢讲,先占坑。
import numpy as np
import matplotlib.pyplot as plt
import os
import seaborn as sns
from scipy import stats
import pickle as pk
import IterEMD as temd
plt.style.use('Solarize_Light2')
# Load dummy signal
X, signals, signal_colors, sample_rate = temd.get_figure_1()
# Plot the dummy signal, which is contructed
plt.figure(figsize=(30,8))
temd.plot_emd(signals, sample_rate, X=X, spaceFactor=0.8, imfs_shift=12, lw_imfs=2, lw_X=2, imf_cols=signal_colors)
# Run for a single mask frequency - see the agrgument f_set
Xs = [X]
mixScore_func = temd.get_modeMixScore_corr
it_mask_freqs, it_mix_scores, it_adj_mix_scores, it_consistency_scores, it_is, optimised_mask_freqs, converged = \
temd.run_tmEMD(Xs, sample_rate, mixScore_func=mixScore_func, max_imfs=2, n_per_it=100, top_n=5, nprocesses=4, f_set=[None, 0])
# Visualise and compare mode mixing to EMD variants
variants = ['EMD', 'eEMD', 'mEMD_zc', 'itEMD']
temd.figplot_tmEMD(Xs, 0, it_mask_freqs, it_mix_scores, sample_rate, mixScore_func, spaceFactor=0.8, show_variants=True, variants=variants,cmap='Set2', eg_percs=[100, 80, 50, 0],ms=8, ms_=10, large=True, show_egs=True, pad_egs=True)
# Plot mixing scores as a function
_, label = mixScore_func(None, None, None, compute=False, return_label=True)it_mix_scores_M = it_mix_scores.mean(axis=1)plt.figure(figsize=(20, 4))
plt.plot(np.arange(it_mix_scores.shape[0])+1, it_mix_scores_M, color='k')
plt.xlabel('Sub-iteration')
plt.ylabel(label)
plt.yscale('log')
rootFolder = '' # change this pathsample_rate = 1250.
dataFolder = rootFolder+'data/'
outputFolder = rootFolder+'output/'
nX = 8 # number of animals to load and optimise
region = 'NAc' # brain region recorded
nSecs = 10 # up to 300 secs, for a faster runtime, make lower# Load the local field potential recordings
length = int(nSecs*sample_rate)
X_paths = [dataFolder+'ani_'+str(i+1)+'.'+region+'.lfp'+'.npy' for i in range(nX)]
Xs = [np.load(path)[:length] for path in X_paths]
nSecs = 3
color='k'
w, h = 30, 4
#
nSamples = int(nSecs*sample_rate)plt.figure(figsize=(w, h*nX))
for i, X in enumerate(Xs):plt.subplot(nX, 1, i+1)plt.title('ani_'+str(i+1), loc='left', fontweight='bold')st = np.random.choice(np.arange(len(X)-nSamples))en = st+ nSamplestimeAx_secs = np.linspace(st/sample_rate, en/sample_rate, nSamples)plt.plot(timeAx_secs, X[st:en], color=color)
alpha, lw = 0.4, 2
for X in Xs:X = stats.zscore(X)f, p = temd.get_psd(X, sample_rate)plt.plot(f, p, color=color, alpha=alpha, lw=lw)plt.xscale('log')
nprocesses = 4it_mask_freqs, it_mix_scores, it_adj_mix_scores, it_consistency_scores, it_is, optimised_mask_freqs, converged = \
temd.run_tmEMD(Xs, sample_rate, n_per_it=20, nprocesses=nprocesses, max_iterations=4, max_imfs=8)
plt.plot(it_mix_scores.mean(axis=1))
xi=0
temd.figplot_tmEMD(Xs, xi, it_mask_freqs, it_mix_scores, sample_rate, mixScore_func, show_variants=False, eg_percs=[0, 30, 80], opt2xi=True)
xi = 0
variants = ['EMD', 'eEMD', 'mEMD_zc', 'itEMD']
temd.figplot_tmEMD(Xs, xi, it_mask_freqs, it_mix_scores, sample_rate, mixScore_func, log_mixScore=True, show_egs=True, opt2xi=True, eg_percs=[0, 40, 95], variants=variants, cmap='husl',fontsize=20, show_variants=True, window=[22220, 25220])
w_emd = 18
w_psd = 6
h = 6
facecolor=NonewTot = w_emd + w_psd
hTot = h * len(variants)plt.figure(figsize=(wTot, hTot))
grid = plt.GridSpec(hTot, wTot, hspace=3, wspace=3)currH = 0
for variant in variants:imfs = temd.run_emd(Xs[xi], sample_rate, variant)imf_cols = sns.color_palette('husl', imfs.shape[1])[::-1]plt.subplot(grid[currH:(currH+h), :w_emd], facecolor=facecolor)plt.title(variant, loc='left', fontweight='bold', fontsize=18)temd.plot_emd(imfs, sample_rate, window=[22220, 25220], imf_cols=imf_cols)plt.subplot(grid[currH:(currH+h), w_emd:], facecolor=facecolor)freqAx_psd, imfPSDs = temd.get_imfPSDs(imfs, sample_rate)temd.plot_imfPSDs(freqAx_psd, imfPSDs, imf_cols=imf_cols)currH += h
知乎学术咨询:
https://www.zhihu.com/consult/people/792359672131756032?isMe=1
擅长领域:现代信号处理,机器学习,深度学习,数字孪生,时间序列分析,设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。