使用经验模态分解（EMD）处理振动信号，并结合样本熵进行特征提取

使用经验模态分解（EMD）处理振动信号，并结合样本熵进行特征提取，是机械故障诊断等领域非常有效的方法。从信号分解 → 阈值去噪 → 特征提取的完整流程

技术路线总览

从原始振动信号到最终特征向量的完整处理流程：

核心步骤详解与MATLAB实现

1. EMD分解振动信号
EMD能将复杂的振动信号自适应地分解为一系列本征模函数（IMF），这是处理非线性、非平稳振动信号的关键。

% 假设已有振动信号 data, 采样频率 fs
[imf, residue] = emd(data, 'Interpolation', 'pchip', 'Display', 0); 
% imf: 分解得到的各阶IMF分量（每一行是一个IMF）
% residue: 残余分量

2. 关键难点：阈值的确定
这是处理中的核心挑战，主要用于从IMF中分离有效成分与噪声。以下是几种实用方法：

3. 样本熵的提取
样本熵用于量化时间序列的复杂性和规律性，非常适合描述振动信号的复杂度特征。

function SampEnVal = compute_sampen(data, m, r)
    % 计算样本熵
    % data: 输入时间序列（如某个IMF）
    % m: 嵌入维度，通常取1或2
    % r: 相似容限，通常取0.1~0.25倍的标准差(std)

    N = length(data);
    lastrun = zeros(1, N);
    run = zeros(1, N);
    A = 0; B = 0;
    p = zeros(1, N);
    e = zeros(1, N);

    % 计算向量匹配个数
    for i = 1:N-m
        for j = i+1:N-m
            if max(abs(data(i:i+m-1) - data(j:j+m-1))) < r
                run(j) = lastrun(j) + 1;
                A = A + run(j);
                B = B + 1;
            else
                run(j) = 0;
            end
            lastrun(j) = run(j);
        end
        for j = 1:N
            lastrun(j) = 0;
            run(j) = 0;
        end
    end

    if A==0 || B==0
        SampEnVal = -log(1/((N-m)*(N-m-1)));
    else
        SampEnVal = -log(A/B);
    end
end

% 对选定的IMF计算样本熵
m = 2; % 嵌入维度
r = 0.2 * std(imf_selected); % 容限
sampen_value = compute_sampen(imf_selected, m, r);

完整MATLAB处理流程示例

将上述步骤整合，形成可直接运行或修改的完整流程：

%% 完整流程：振动信号EMD处理与样本熵特征提取
clear; clc; close all;

% 1. 加载或生成振动信号（示例：含噪声的调频调幅信号）
fs = 1000; % 采样频率
t = 0:1/fs:2-1/fs;
data = 2*sin(2*pi*20*t) .* (1+0.5*sin(2*pi*5*t)) + 0.8*sin(2*pi*80*t) + 0.5*randn(size(t));

% 2. EMD分解
[imf, residue] = emd(data, 'Interpolation', 'pchip', 'Display', 0);
num_imf = size(imf, 1);

% 3. 阈值确定（以能量占比法为例）
energy = sum(imf.^2, 2);
total_energy = sum(energy);
cum_energy_ratio = cumsum(energy) / total_energy;
threshold_ratio = 0.85; % 能量阈值设为85%
k = find(cum_energy_ratio >= threshold_ratio, 1); % 选取前k个主要IMF
imf_primary = imf(1:k, :); % 有效成分
imf_noise = imf(k+1:end, :); % 视为噪声为主的分量（可选后续处理）

% 4. 对每个主要IMF计算样本熵
m = 2; % 嵌入维度
feature_vector = zeros(1, k); % 初始化特征向量
for i = 1:k
    r = 0.2 * std(imf_primary(i, :)); % 容限
    feature_vector(i) = compute_sampen(imf_primary(i, :), m, r);
end

% 5. 可视化
figure('Position', [100, 100, 1200, 800]);
% 5.1 原始信号与IMF
subplot(3, 1, 1);
plot(t, data); title('原始振动信号'); xlabel('时间 (s)'); ylabel('幅值');
subplot(3, 1, 2);
plot(t, imf_primary'); title('主要IMF分量'); xlabel('时间 (s)'); ylabel('幅值');
% 5.2 样本熵特征
subplot(3, 1, 3);
bar(1:k, feature_vector);
title('各主要IMF的样本熵特征值'); xlabel('IMF序号'); ylabel('样本熵');
grid on;

参考代码 emd进行振动信号处理，阀值的确定，样本熵的提取 www.youwenfan.com/contentalh/54949.html

关键参数选择与物理意义

• EMD参数：MATLAB内置emd函数通常默认设置即可。若出现过度分解，可调整'MaxNumIMF'。
• 样本熵参数：
  • 嵌入维度 m：通常取1或2，表示比较向量的长度。
  • 容限 r：通常取0.1到0.25倍时间序列的标准差。r值过小会估计偏差大，过大则会丢失有效信息。建议对同类信号固定r值，以保证特征可比性。
• 物理意义：在故障诊断中，样本熵值降低可能意味着信号复杂度下降，规律性增强，有时对应于早期故障引起的周期性冲击成分增多。