改进经验模态分解方法-通过迭代方式（IMF振幅加权频率，Python）

一种新颖的改进经验模态分解方法-通过迭代方式（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等。