一种新颖的改进经验模态分解方法-通过迭代方式(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 path
sample_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+ nSamples
timeAx_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 = 4
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, 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=None
wTot = 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等。