基于EMD分解的IMF分量优选及多尺度熵特征提取方法研究

经验模态分解EMDMatlab 对信号EMD分解 IMF分量图频谱图功率谱能量谱等 选择最优Imf分量重构信号也可以在EMD的基础上提取相关时域、频域、熵值特征 多尺度熵咱们今天拿个很接地气的例子开搞——就是你运动手表腕上晃出来的原始三轴加速度。这种信号哪是单纯的摆臂/步频啊混了手腕甩、手表松垮、甚至可能有静电干扰或者地面坑洼的碎抖直接用来算配速或者能量消耗肯定飘。这时候经验模态分解EMD就该上场救场了它不像傅里叶那种硬分“正弦波块”而是顺着信号本身的“振动习惯”拆成一个个固有模态函数IMF拆完再挑好用的分量拼回去干净得一批。先上咱这次模拟的原始信号吧别真去戴手表测了写个Matlab代码凑个带真实感的杂信号clear; clc; close all; fs 100; % 采样率100Hz运动手表差不多这水平 t 0:1/fs:10-1/fs; % 取10秒数据 % 真实步频中年慢跑约1.8Hz对应摆臂频率0.9Hz arm_swing 0.5*sin(2*pi*0.9*t); % 手表没带紧的高频小抖大概20Hz watch_wobble 0.2*sin(2*pi*20*t 0.3*pi); % 偶尔踩个小石子的脉冲尖峰随机放3个 spike_loc randsample(length(t),3); spike_signal zeros(size(t)); spike_signal(spike_loc) 2; % 再加一点点白噪声意思意思 white_noise 0.05*randn(size(t)); % 拼起来就是咱们的原始加速度x轴假设主要动在x轴 raw_signal arm_swing watch_wobble spike_signal white_noise; % 先画原始信号瞅瞅多乱 figure(Color,w,Position,[100 100 1000 400]); plot(t, raw_signal, k-, LineWidth,0.8); xlabel(时间 (s),FontSize,12); ylabel(归一化加速度,FontSize,12); title(模拟中年慢跑手表松小石子白噪的x轴加速度,FontSize,14); grid on;你看这图步频的小波浪都快被淹没在杂波里了对吧接下来就直接调Matlab自带的emd函数就行不用自己写那套复杂的筛选循环了当然要写的话可以练手但自带的效率高还稳。记得2021b之后的版本好像默认是带改进的EEMD修正的没关系先直接跑基础EMD看看结果% 跑EMD分解成IMF和残差r [imf, r] emd(raw_signal, Interpolation,pchip); % 用pchip插值更适合带尖峰的信号 n_imf size(imf,2); % 画所有IMF和残差 figure(Color,w,Position,[100 100 1200 200*(n_imf1)]); subplot(n_imf1,1,1); plot(t, raw_signal, k-, LineWidth,0.8); ylabel(原始信号,FontSize,11); title(EMD分解结果,FontSize,14); grid on; for i 1:n_imf subplot(n_imf1,1,i1); plot(t, imf(:,i), b-, LineWidth,0.7); ylabel([IMF,num2str(i)],FontSize,11); grid on; end subplot(n_imf1,1,n_imf2); % 不对刚才循环加1了直接i循环到n_imf1哦刚才写错了循环的subplot索引 % 重新修正一下画残差的部分 end subplot(n_imf1,1,n_imf1); plot(t, r, r-, LineWidth,0.8); ylabel(残差 r,FontSize,11); xlabel(时间 (s),FontSize,11); grid on;修正后出来的IMF大概长啥样IMF1肯定是最高频的手表抖和小石子尖峰、白噪的一部分应该在这儿然后往下频率越来越低大概IMF3或者IMF4就是咱们要的0.9Hz左右的摆臂/步频了残差就是信号的整体趋势比如模拟的时候没加趋势残差应该是个很小的接近0的直线。接下来到关键步骤怎么挑最优的IMF不能光靠眼睛瞅也得有点客观指标对吧咱们可以搞几个常用的相关系数Correlation Coefficient看这个IMF和原始信号的线性相关程度越高说明包含的原始有效信息越多。但太高频的比如IMF1虽然和噪声尖峰相关但和步频没用所以得结合下一个。平均频率Mean Frequency算IMF的功率谱重心大概对应这个分量的主要振动频率咱们要0.9±0.2Hz的范围。排列熵Permutation Entropy, PE熵值越大数据越乱步频的IMF排列熵应该比高频杂波低很多因为它是周期性的。先写这三个指标的计算代码顺便连每个IMF的频谱图、功率谱、能量谱都画一下能量谱其实就是功率谱对时间积分或者说傅里叶变换幅值的平方有时候大家会混淆但运动场景看功率谱重心就行% 先定义计算排列熵的小函数Matlab自带的话可能在Statistics and Machine Learning Toolbox没有的话自己写这个简单的就行 function pe permutation_entropy(x, m, tau) % x: 输入信号 % m: 嵌入维度一般取3-7 % tau: 延迟时间一般取1 N length(x); x_embed zeros(N - (m-1)*tau, m); for i 1:m x_embed(:,i) x((i-1)*tau1 : N - (m-1)*tau (i-1)*tau); end % 对每个嵌入向量排序得到符号序列 symbol_seq zeros(size(x_embed,1),1); for i 1:size(x_embed,1) [~, idx] sort(x_embed(i,:)); symbol_seq(i) sum((idx - 1) .* (m.^(0:m-1))); % 把符号序列转成唯一的数字索引 end % 计算每个符号的概率 prob histcounts(symbol_seq, 0:m^m) / length(symbol_seq); prob prob(prob 0); % 去掉概率为0的 pe -sum(prob .* log2(prob)) / log2(m^m); % 归一化到0-1 end % 开始批量算指标、画图 m_pe 5; tau_pe 1; corr_vals zeros(n_imf,1); mean_freqs zeros(n_imf,1); pe_vals zeros(n_imf,1); figure(Color,w,Position,[150 150 1400 600*(ceil(n_imf/2))]); tiledlayout(ceil(n_imf/2),4,TileSpacing,compact,Padding,compact); for i 1:n_imf % 取当前IMF current_imf imf(:,i); % 1. 时域图不用单独再画大的了嵌在频谱图旁边 nexttile; plot(t, current_imf, b-, LineWidth,0.6); xlabel(时间 (s),FontSize,9); ylabel([IMF,num2str(i)],FontSize,10); title([时域], FontSize,10); grid on; % 2. 频谱图用fft直接算加汉宁窗防泄漏 nexttile; n_fft 2^nextpow2(length(current_imf)); f fs/2 * linspace(0,1,n_fft/21); win hann(length(current_imf)); current_imf_win current_imf .* win; fft_vals fft(current_imf_win, n_fft); amp_spectrum abs(fft_vals(1:n_fft/21)) * 2 / sum(win); plot(f, amp_spectrum, g-, LineWidth,0.6); xlabel(频率 (Hz),FontSize,9); ylabel(幅值,FontSize,10); title([频谱], FontSize,10); xlim([0, fs/2]); grid on; % 3. 功率谱幅值平方单位是能量/Hz这里归一化时间/信号就行 nexttile; power_spectrum amp_spectrum.^2; plot(f, power_spectrum, m-, LineWidth,0.6); xlabel(频率 (Hz),FontSize,9); ylabel(功率,FontSize,10); title([功率谱], FontSize,10); xlim([0, fs/2]); grid on; % 4. 能量谱功率谱对频率积分的累积还是幅值平方的总和在时域分布哦这里说的能量谱通常是指短时傅里叶的能量或者就简单说总能量累积到残差前越来越少算了换个大家更常用的每个IMF的**归一化能量占比**直接标在这个tile的标题或者备注里吧同时算另外三个指标 nexttile; % 算指标 corr_vals(i) corr(current_imf, raw_signal); mean_freqs(i) sum(f .* power_spectrum) / sum(power_spectrum); % 功率谱重心就是平均频率 pe_vals(i) permutation_entropy(current_imf, m_pe, tau_pe); % 归一化能量占比 total_energy sum(raw_signal.^2); imf_energy sum(current_imf.^2); energy_ratio imf_energy / total_energy * 100; % 画个简单的指标可视化柱形图小合集或者直接文字标出来方便看 axis off; text(0.1,0.8, [相关系数: , num2str(corr_vals(i), %.3f)], FontSize,10, Color,k); text(0.1,0.6, [平均频率: , num2str(mean_freqs(i), %.2f), Hz], FontSize,10, Color,k); text(0.1,0.4, [排列熵: , num2str(pe_vals(i), %.3f)], FontSize,10, Color,k); text(0.1,0.2, [能量占比: , num2str(energy_ratio, %.2f), %], FontSize,10, Color,k); title([指标汇总], FontSize,10); end好现在看指标和小图的话应该很清楚了比如刚才模拟的情况IMF1平均频率18-22Hz排列熵0.7-0.9很高乱相关系数可能0.3左右IMF2可能10Hz左右还是杂IMF3或者IMF4平均频率刚好0.8-1.1Hz排列熵0.2-0.4很低周期强相关系数0.6-0.9很高有效信息多能量占比可能20-30%。那最优的就是这个IMF经验模态分解EMDMatlab 对信号EMD分解 IMF分量图频谱图功率谱能量谱等 选择最优Imf分量重构信号也可以在EMD的基础上提取相关时域、频域、熵值特征 多尺度熵接下来就重构信号只加这个最优IMF就行如果有趋势残差有用的话也可以加但咱们模拟的没有趋势% 假设刚才算出来最优的是IMF3实际运行看你的指标选索引哦 best_imf_idx 3; reconstructed_signal imf(:, best_imf_idx); % 对比原始、最优IMF、重构后的信号 figure(Color,w,Position,[100 100 1000 600]); subplot(3,1,1); plot(t, raw_signal, k-, LineWidth,0.8); ylabel(原始信号,FontSize,12); title(信号对比原始、最优IMF、重构,FontSize,14); grid on; subplot(3,1,2); plot(t, imf(:,best_imf_idx), b-, LineWidth,0.8); ylabel([最优IMF,num2str(best_imf_idx)],FontSize,12); grid on; subplot(3,1,3); plot(t, reconstructed_signal, r-, LineWidth,0.8); ylabel(重构信号,FontSize,12); xlabel(时间 (s),FontSize,12); grid on;完美重构后的信号就是干净的摆臂小波浪配速或者能量消耗的算法往里面一套肯定准多了。既然提到了EMD基础上提取特征最后再补个多尺度熵MSE的小尾巴——单尺度排列熵只能看某一个“观察粒度”下的信号复杂度多尺度就是先把信号做不同尺度的粗粒化比如尺度2就是每2个点取平均尺度3每3个点然后算每个粗粒化信号的排列熵这样能看出信号在不同时间尺度下的变化比如步频信号的多尺度熵在小尺度可能有点抖因为采样的小波动但大尺度就稳定低而纯噪声的多尺度熵所有尺度都高有趋势的信号多尺度熵会随尺度下降。多尺度熵的Matlab代码同样基于刚才的排列熵小函数function mse multiscale_entropy(x, m, tau, max_scale) % x: 输入信号 % m: 嵌入维度 % tau: 延迟时间 % max_scale: 最大粗粒化尺度 mse zeros(max_scale,1); for s 1:max_scale % 粗粒化每s个点取平均 N length(x); x_coarse zeros(1, floor(N/s)); for i 1:floor(N/s) x_coarse(i) mean(x((i-1)*s1 : i*s)); end % 算这个粗粒化信号的排列熵 mse(s) permutation_entropy(x_coarse, m, tau); end end % 测试一下算原始信号、最优IMF、重构信号的多尺度熵 max_scale 10; m_pe 5; tau_pe 1; mse_raw multiscale_entropy(raw_signal, m_pe, tau_pe, max_scale); mse_best multiscale_entropy(imf(:,best_imf_idx), m_pe, tau_pe, max_scale); mse_rec multiscale_entropy(reconstructed_signal, m_pe, tau_pe, max_scale); % 画多尺度熵曲线 figure(Color,w,Position,[100 100 800 500]); plot(1:max_scale, mse_raw, k-o, LineWidth,1.2, MarkerSize,6, DisplayName,原始信号); hold on; plot(1:max_scale, mse_best, b-s, LineWidth,1.2, MarkerSize,6, DisplayName,最优IMF); plot(1:max_scale, mse_rec, r-^, LineWidth,1.2, MarkerSize,6, DisplayName,重构信号); hold off; xlabel(粗粒化尺度 s,FontSize,12); ylabel(归一化多尺度排列熵,FontSize,12); title(不同信号的多尺度熵对比,FontSize,14); legend(Location,best); grid on;这个曲线应该能很明显看出三者的区别原始信号的熵所有尺度都比另外两个高最优IMF和重构信号的熵几乎重合因为重构就是它自己小尺度可能有个小波动因为原始采样的一点点残留但粗粒化后很快稳定在0.2左右的低熵。好了今天这一套从模拟杂信号、EMD分解、指标选最优、重构、到时域/频域/熵值特别是多尺度特征提取应该够用了吧有兴趣的话可以自己换个真实的生理信号比如心电、脑电或者机械振动信号试试EMD对非平稳非线性信号真的友好