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

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

技术路线总览

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

flowchart LR A[“原始振动信号”] --> B[“EMD分解<br>获得若干IMF分量”] B --> C{“选择需要处理的IMF”} C --> D[“阈值去噪处理<br>（关键步骤）”] C --> E[“直接计算样本熵”] D --> F[“计算去噪后IMF的样本熵”] E --> G[“形成特征向量”] F --> G

核心步骤详解与MATLAB实现

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

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

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

方法	核心思想	适用场景与MATLAB关键步骤
基于能量占比法	认为前几个能量高的IMF主要包含有效信号，后续IMF主要为噪声。	最常用，适用于大多数旋转机械振动信号。 计算各IMF能量：energy = sum(imf.^2, 2); 设定累积能量阈值（如85%），选择能量占比达到阈值的前k个IMF。
基于统计特性的固定阈值	对每个IMF，计算其标准差，设定一个乘数（如2~3倍）作为阈值。	适合噪声特性相对稳定的场景。 threshold = k * std(imf(i,:)); 对小波阈值函数（如wthresh）也可调整使用。
基于排列熵的自适应阈值	排列熵能衡量序列的随机性，值越大越可能是噪声。	适合噪声复杂或信号非平稳性强的场景。 计算各IMF排列熵，设定熵阈值（如0.6），高于阈值的IMF被视为需处理的噪声主导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/contentcnq/54949.html

关键参数选择与物理意义

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