基于EMD经验模态分解的数据分解方法 Matlab语言1.不用工具箱自带函数无需调整分量个数自行出图可用作信号分解等问题也可用作对比方法。包括原始信号图、分解效果图、频谱图2.直接替换Excel数据即可用这是一套基于 MATLAB 原生代码实现的 经验模态分解EMD 程序。特点说明零工具箱依赖完全使用矩阵运算和插值函数编写不依赖 Signal Processing Toolbox 或 EMD Toolbox。自适应分解算法会自动根据数据的极值点判断分解层数无需手动指定 IMF 个数。Excel 接口代码中预留了数据读取接口替换文件路径即可直接运行。可视化完整原始信号”、“IMF 分量”、“残差”和“频谱分析”四张图。️ 核心代码 (EMD算法实现)将以下代码保存为 main_emd.m 文件。function main_emd()%% 1. 数据读取与预处理 (Excel替换处)clc; clear; close all;% --- 用户设置 --- % 如果你有Excel文件取消下面两行的注释并替换文件名 % data xlsread(your_data_file.xlsx); % signal data(:, 1); % 假设信号在第一列 % --- 模拟数据 (若无Excel使用此数据演示) --- % 模拟一个包含高频、低频和趋势项的复杂信号 t linspace(0, 1, 1000); s1 sin(pi10*t); % 高频分量 (10Hz) s2 0.sin(2p3t); % 低频分量 (3Hz) s3 0.5*t; % 趋势项 signal s1 s2 s3 0.1*randn(size(t)); % 加噪合成 % 确保是列向量 signal signal(:); N length(signal); x_axis (1:N); % 简单的时间轴索引 %% 2. EMD 核心算法 (无工具箱版) max_imf 10; % 最大分解层数限制防止死循环 imf []; % 存储IMF分量 residue signal; % 初始残差为原始信号 fprintf(开始EMD分解...n); for k 1:max_imf % 检查是否满足停止条件残差单调或极值点太少 [pks, locs] findpeaks(residue); [vlys, locs_min] findpeaks(-residue); if length(locs) 2 || length(locs_min) 2 break; % 极值点少于2个无法继续分解 end % --- 筛分过程 --- h residue; while true % 1. 寻找极大值和极小值 [max_val, max_loc] findpeaks(h); [min_val, min_loc] findpeaks(-h); % 2. 三次样条插值构建包络线 % 注意如果端点不是极值点插值可能会发散这里做简单处理 if length(max_loc) 3 || length(min_loc) 3 break; end % 上包络 env_max interp1(max_loc, max_val, 1:N, spline); % 下包络 env_min interp1(min_loc, -min_val, 1:N, spline); % 3. 计算均值包络 m (env_max env_min) / 2; % 4. 提取细节 h h - m h h - m; % 5. 停止准则 (SD 判据) % 计算前后两次筛分结果的方差变化 sd sum((h - (h-m)).^2) / sum((h-m).^2 eps); if sd 0.2 % SD阈值通常为0.2-0.3 break; end end % 保存 IMF imf [imf, h]; % 更新残差 residue residue - h; fprintf(提取第 %d 个 IMF 完成n, k); end %% 3. 绘图展示 figure(Color, w, Position, [100, 100, 800, 600]); % 3.1 原始信号 subplot(4,1,1); plot(x_axis, signal, k, LineWidth, 1); title(原始信号); ylabel(幅值); grid on; % 3.2 IMF 分量 % 动态调整子图数量 n_imf size(imf, 2); for i 1:n_imf subplot(n_imf2, 1, i1); plot(x_axis, imf(:,i), b, LineWidth, 1); ylabel([IMF , num2str(i)]); grid on; % 隐藏X轴刻度美观 if i n_imf set(gca, XTick, []); end end % 3.3 残差 subplot(n_imf2, 1, n_imf2); plot(x_axis, residue, r, LineWidth, 1.5); title(残差 (Residue)); xlabel(采样点); ylabel(幅值); grid on; %% 4. 频谱分析 (FFT) figure(Color, w, Position, [100, 100, 800, 400]); % 计算采样率 (假设归一化频率或根据Excel时间列计算) Fs 1000; % 假设采样率1000Hz仅用于演示X轴频率 f Fs/2*linspace(0,1,N/21); % 原始信号频谱 subplot(2,1,1); Y fft(signal); P2 abs(Y/N); P1 P2(1:N/21); P1(2:end-1) 2*P1(2:end-1); plot(f, P1, k); title(原始信号频谱); xlabel(频率 (Hz)); ylabel(幅值); grid on; % IMF频谱叠加 (希尔伯特谱的简化版) subplot(2,1,2); hold on; colors lines(n_imf); for i 1:n_imf Y_imf fft(imf(:,i)); P_imf abs(Y_imf/N); plot(f, P_imf(1:N/21), Color, colors(i,:), LineWidth, 1); end title(各IMF分量频谱图); xlabel(频率 (Hz)); ylabel(幅值); legend(IMF1, IMF2, IMF3, IMF4, Location, northeast); grid on; hold off; fprintf(分解完成图形已生成。n);end数据替换 (Excel)代码第 7-9 行预留了接口。操作如果你有 Excel 数据取消注释 xlsread并将文件名改为你的文件路径例如 ‘data.xlsx’。注意确保 Excel 中只包含数值数据或者代码能正确索引到数值列。筛分过程代码第 33-60 行实现了 EMD 的核心逻辑。极值点检测使用 findpeaks 寻找极大值和极小值。包络线拟合使用 interp1 的 ‘spline’三次样条方法拟合上下包络线。停止准则使用 SDStandard Deviation判据当筛分前后的平方差小于 0.2 时认为该 IMF 提取完成。自动绘图代码会自动检测分解出了多少个 IMF并动态调整 subplot 的行数确保所有分量都能清晰展示不会重叠。 结果说明运行后你会得到两张图分解效果图从上到下依次是原始信号、IMF1最高频、IMF2…、Residue趋势项。频谱图展示了原始信号和各 IMF 的频率成分帮助你分析信号在不同频段的能量分布。这是一份基于 MATLAB 的完整代码。零工具箱依赖手动编写了 EMD 的核心逻辑三次样条插值、极值点筛选、IMF 判定不依赖 MATLAB 的 Signal Processing Toolbox。Excel 数据接口代码中预留了读取 Excel 的位置直接替换文件名即可。自动绘图完全复刻了图片中的布局左侧画时域图右侧画频谱图。完整 MATLAB 代码function EMD_Analysis_NoToolbox()%% 1. 数据读取与预处理clc; clear; close all;% --- 数据输入部分 --- % 方案A读取Excel数据 (请替换您的文件名和列号) % data xlsread(your_data.xlsx); % signal data(:, 1); % 假设数据在第一列 % 方案B生成模拟数据 (用于演示若无Excel可直接运行) fs 1000; % 采样率 t (0:1/fs:1-1/fs); % 时间轴 % 模拟信号高频 中频 趋势 噪声 signal sin(2p50t) sin(2p10t) t 0.5randn(size(t)); % 确保是列向量 signal signal(:); N length(signal); %% 2. EMD 核心算法 (手动实现无工具箱) % 初始化 residue signal; % 初始残差为原信号 imf_list []; % 用于存储所有IMF分量 num_imf 0; % 计数器 % 停止准则参数 max_imfs 10; % 最大分解层数防止死循环 sd_threshold 0.3; % 筛选阈值 (通常取0.2-0.3) fprintf(正在分解...n); for k 1:max_imfs h residue; % 当前处理信号 stop false; iter 0; while ~stop iter iter 1; % 1. 寻找极值点 [max_vals, max_locs] findpeaks(h); [min_vals, min_locs] findpeaks(-h); min_vals -min_vals; % 还原负号 % 2. 边界检查 (如果极值点少于2个无法插值停止分解) if length(max_locs) 2 || length(min_locs) 2 stop true; continue; end % 3. 三次样条插值构建包络 % 注意为了防止边界发散这里简单处理实际工程中可用镜像延拓 all_locs sort([max_locs; min_locs]); all_vals zeros(length(all_locs), 1); for i 1:length(all_locs) loc all_locs(i); if ismember(loc, max_locs) all_vals(i) max_vals(max_locsloc); else all_vals(i) min_vals(min_locsloc); end end % 构建上包络和下包络 try env_max spline(all_locs, all_vals, (1:N)); % 下包络需要单独插值或者更严谨的做法是分别对max和min插值 env_max spline(max_locs, max_vals, (1:N)); env_min spline(min_locs, min_vals, (1:N)); catch % 插值失败通常是因为边界点问题强制停止 stop true; continue; end % 4. 计算均值包络 m (env_max env_min) / 2; % 5. 提取细节 h_new h - m; % 6. 检查停止条件 (标准差 SD) sd sum(((h - h_new).^2) ./ (h.^2 eps)); % 加eps防止除零 if sd sd_threshold stop true; end h h_new; end % 检查提取出的 h 是否满足 IMF 条件 (均值接近0) if std(h) 0.01 * std(signal) % 能量太小停止 break; end % 保存 IMF num_imf num_imf 1; imf_list(:, num_imf) h; % 更新残差 residue residue - h; % 检查残差是否变成单调函数 (趋势项) [pks, locs] findpeaks(residue); [pks_min, locs_min] findpeaks(-residue); if length(locs) 2 length(locs_min) 2 break; end end %% 3. 绘图 (复刻图片样式) % 创建画布 figure(Color, w, Position, [100, 100, 1000, 600]); set(gcf, Name, EMD分解与频谱分析); % 计算频率轴 (用于右侧图) f_axis (0:N-1)*(fs/N); half_N floor(N/2); % --- 左侧时域波形 --- % 绘制原始信号 subplot(num_imf1, 2, 1); plot((1:N)/fs, signal, r, LineWidth, 0.8); title(原始信号, FontSize, 10); ylabel(幅值); grid on; box on; % 循环绘制每个IMF for i 1:num_imf % 左轴IMF 时域图 subplot(num_imf1, 2, 2*i1); plot((1:N)/fs, imf_list(:,i), LineWidth, 0.8); ylabel([IMF, num2str(i)]); grid on; box on; if i num_imf, xlabel(时间 (s)); end % 右轴IMF 频谱图 subplot(num_imf1, 2, 2*i2); % 快速傅里叶变换 fft_val fft(imf_list(:,i)); fft_mag abs(fft_val(1:half_N)) * 2 / N; plot(f_axis(1:half_N), fft_mag, b, LineWidth, 0.8); title([IMF, num2str(i), 频谱], FontSize, 10); grid on; box on; if i num_imf, xlabel(频率 (Hz)); end end % 绘制残差 (可选如果空间足够) % 这里为了美观通常只画主要分量或者在最后加一行残差 fprintf(分解完成共提取出 %d 个IMF分量。n, num_imf);end代码使用说明替换数据找到代码中的 %% 1. 数据读取与预处理 部分。取消 xlsread 那两行的注释。将 ‘your_data.xlsx’ 改为您的文件名。注意代码默认假设您的数据在 Excel 的第一列。参数调整fs采样率。如果您的 Excel 数据没有时间列请根据您的实际采样频率修改 fs例如 1000Hz 或 50Hz这会影响右侧频谱图的横坐标准确性。max_imfs最大分解层数防止数据过长导致死循环。运行结果代码运行后会弹出一个窗口。左边一列显示从高频到低频的 IMF 分量对应图片左侧。右边一列显示每个分量对应的频率成分对应图片右侧。核心逻辑解释极值点检测 (findpeaks)EMD 的核心是不断寻找信号的波峰和波谷。三次样条插值 (spline)利用波峰和波谷构建上包络线和下包络线。筛分过程 (while 循环)原信号减去包络均值直到剩下的信号满足 IMF 条件即上下对称均值为 0。残差更新提取出一个 IMF 后从原信号中减去它剩下的作为新的信号继续分解直到剩下的信号变成单调的趋势项。