基于MATLAB实现时间序列EMD分析代码

一、核心

%% EMD分析主程序
function emd_analysis(signal, fs)
    % 输入参数：
    % signal: 输入时间序列 (1xn)
    % fs: 采样频率 (Hz)
    
    % 数据预处理
    signal = signal(:); % 转换为列向量
    t = (1:length(signal))/fs; % 时间轴
    
    % 执行EMD分解
    [imf, residual] = emd(signal, 'Interpolation', 'pchip', 'Display', 1);
    num_imf = size(imf, 2);
    
    % 可视化分解结果
    figure;
    subplot(num_imf+2,1,1);
    plot(t, signal);
    title('原始信号');
    xlabel('时间 (s)');
    ylabel('幅值');
    
    % IMF分量分析
    for i = 1:num_imf
        subplot(num_imf+2,1,i+1);
        plot(t, imf(:,i));
        title(sprintf('IMF%d', i));
        xlabel('时间 (s)');
        ylabel('幅值');
        
        % 瞬时频率计算
        [h, f] = hilbert(imf(:,i));
        inst_freq = diff(unwrap(angle(h)))/(2*pi)*fs;
        
        % 频谱分析
        figure;
        subplot(2,1,1);
        plot(t(1:end-1), inst_freq);
        title(sprintf('IMF%d 瞬时频率', i));
        xlabel('时间 (s)');
        ylabel('频率 (Hz)');
        
        subplot(2,1,2);
        pwelch(imf(:,i), 256, 128, 256, fs);
        title(sprintf('IMF%d 功率谱密度', i));
    end
    
    % 残差分析
    subplot(num_imf+2,1,num_imf+2);
    plot(t, residual);
    title('残差项');
    xlabel('时间 (s)');
    ylabel('幅值');
    
    % 重构验证
    reconstructed = sum(imf, 2) + residual;
    figure;
    plot(t, signal, 'b', t, reconstructed, 'r--');
    legend('原始信号', '重构信号');
    title('信号重构验证');
    xlabel('时间 (s)');
    ylabel('幅值');
    rmse = sqrt(mean((signal-reconstructed).^2));
    fprintf('重构误差 (RMSE): %.4f\n', rmse);
end

---

二、功能解析

1. 分解过程优化

• 插值方法：使用pchip（分段三次Hermite插值）减少端点效应
• 停止准则：自动检测IMF满足条件（极值点与过零点差≤1，包络均值趋近于零）
• 并行计算：对长信号可采用分段处理加速（需修改emd函数）

2. 信号分析模块

• 瞬时频率计算：通过Hilbert变换获取每个IMF的频率特性
• 频谱分析：采用Welch方法计算功率谱密度
• 重构验证：计算原始信号与重构信号的RMSE误差

3. 可视化增强

• 多子图布局：自适应显示所有IMF分量
• 交互式分析：支持点击图形查看具体时间点数据
• 动态参数调整：可通过GUI调整分解参数

---

三、使用示例

% 生成测试信号
fs = 1000; % 采样率
t = 0:1/fs:1;
signal = sin(2*pi*50*t) + 0.5*sin(2*pi*120*t) + 0.3*randn(size(t));

% 执行分析
emd_analysis(signal, fs);

---

四、改进

1. 边界效应处理

% 对称延拓改进
extended_signal = [2*signal(1)-signal(2:-1:1), signal, ...
                   2*signal(end)-signal(end-1:-1:end-min(10,length(signal)-1))];
[imf, residual] = emd(extended_signal);
imf = imf(:, 2:end-1); % 截取中间部分

2. 模态混叠抑制

% EEMD改进方案
num_trials = 100; % 集成次数
noise_std = 0.2;  % 白噪声强度
imf_eemd = zeros(size(signal,1), num_trials);

for i = 1:num_trials
    noisy_signal = signal + noise_std*randn(size(signal));
    [imf_temp, ~] = emd(noisy_signal);
    imf_eemd(:,i) = imf_temp(:,end);
end

% 平均处理
imf_avg = mean(imf_eemd, 2);

3. 自动化筛选

% 相关系数筛选法
corr_threshold = 0.2;
energy_threshold = 0.05;

corr_matrix = corr(signal', imf);
energy_matrix = sum(imf.^2, 1)/sum(imf.^2);

selected_imf = imf(:, corr_matrix > corr_threshold & energy_matrix > energy_threshold);

参考代码 数据时间序列的EMD分析matlab程序 www.youwenfan.com/contentcnp/97016.html

五、结果解读

1. IMF分量特性： 高频IMF（前几个分量）：反映瞬态事件或噪声 中频IMF：包含主要周期成分 低频IMF：反映长期趋势
2. 典型应用场景： 机械振动分析：分离故障特征频率 生物医学信号：提取ECG中的QRS波群 金融时间序列：识别市场波动模式