简介这是一份基于MATLAB实现的二维经验模态分解BEMD算法资源包面向数值算法与人工智能领域的研究者、工程师及高校师生适用于非线性、非平稳信号及多变量耦合数据的分析与处理在环境、金融、生物医学等场景均有应用。压缩包共39个文件以34个m脚本为主体涵盖BEMD核心分解函数、曲面插值、极值点提取、包络构造等关键模块另有4张lena测试图像与1份license说明文档便于直接运行验证算法效果。整体约221KB结构紧凑已有274人学习下载。通过研读运行这套代码读者可掌握从信号预处理、希尔伯特变换、IMF分解到耦合分析、信号重构与可视化的完整流程理解二维分解与传统EMD的差异还可将其中sift、gridfit、extrema等基础工具直接复用灵活迁移至其他数据分析、模式识别任务中。1. BEMD是什么一个把图像拆成“纹理层残差”的MATLAB工具箱入口如果你手头有一张带噪声的纹理图像想把它拆成一层一层的“细节”和一块“基底”BEMD二维经验模态分解可能是你最快的答案。它不是深度神经网络也不依赖训练数据而是靠数值筛选把图像分解成若干个二维固有模态函数BIMF和残差。这个名叫“BEMD.zip”的MATLAB压缩包常见于数值算法课程和人工智能大作业里本质是一套用插值和筛选实现的非线性自适应分解工具。适合谁打算做图像去噪、纹理特征提取、多聚焦图像融合又不想从头写有限差分或小波栈的工程师。它会告诉你什么是“数据驱动的分解”也可能让你踩进边界延拓和计算量的坑。2. BEMD的数值算法原理从EMD到二维筛选为什么用插值而不是神经网络2.1 一维EMD筛选的回顾包络、均值、停止准则BEMD的前身是Huang等人在1998年提出的经验模态分解EMD。给定一个一维信号EMD把它拆成有限个本征模态函数IMF每个IMF代表信号中的一个振荡模式。筛选过程是递归的先找到信号所有局部极大值和局部极小值用三次样条分别插值出上包络和下包络计算两条包络的平均值m然后让原始信号减去m得到候选IMF。如果候选不满足停止准则就以它为输入继续筛直到满足为止。之后从剩余信号中提取下一个IMF直到剩余信号变成单调趋势或很小。这里有两个细节值得记牢。第一极值点不是简单的“找一个点”而是要处理连续相等值的平台区。如果信号在一段时间内完全相同那么平台两端都可以视为极值点但也有实现取平台中心点。选错了包络会出现台阶IMF里会出现伪频率。第二端点效应特别头疼。三次样条插值在端点处没有数据约束包络往往会翘头或下沉导致分解出来的前几个IMF在端点区域失真。我通常会对信号做两端的镜像延拓延拓长度约为信号长度的5%筛完再裁掉。停止准则最常见的是一维的标准差准则SD sum((h_{k-1} - h_k).^2) / sum(h_{k-1}.^2)当SD小于某个阈值时停止典型的阈值是0.2到0.3。但这是Huang经典论文里的建议实际工程中要按噪声水平调整。如果信号本身很干净SD可以设小一点如果噪声大SD太小会把噪声也筛成一个“IMF”结果看着就像把噪声单独抽出来了其实那并不是你想要的。下面是MATLAB示意代码展示一维EMD的单次筛选核心循环但注意它省略了端点延拓和平台区处理真实实现里这些不能省function imf emd_sift_demo(x, maxIter, sdThresh) % 一维EMD单次筛选示例 % 输入x为信号列向量maxIter为最大筛选次数sdThresh为SD停止阈值 % 输出一个IMF分量 h x(:); [pks, locsMax] findpeaks(h); % 局部极大值幅值和位置 [vls, locsMin] findpeaks(-h); % 局部极小值的相反数及位置 for k 1:maxIter envUp spline([1; locsMax], pks, 1:length(h)); % 上包络 envLo spline([1; locsMin], -vls, 1:length(h)); % 下包络 means (envUp envLo) / 2; newH h - means; sd sum((newH - h).^2) / sum(h.^2); if sd sdThresh break; end h newH; end imf h; end逻辑说明每次筛选都减去上下包络的均值让信号逐渐趋向于局部对称。findpeaks的第一个输出是峰值幅值第二个是峰值位置初学者常把顺序搞反导致spline的节点和值对应不上。参数说明maxIter通常取10到20sdThresh取0.2到0.3。若sdThresh太小噪声会被当成独立IMF分离出来。2.2 二维化改造极值点提取、径向基插值与BIMF把二维图像当作一个“平面信号”BEMD的筛选流程就变成寻找图像中的二维局部极大值和极小值点用这些点拟合上包络面和下包络面计算包络平均面从原图中减去重复筛选。得到的每层结果叫二维固有模态函数BIMF最后剩下的是残差。二维极值点提取的难点在于“局部”和平坦区域。对于像素矩阵一个极值点要同时和上下左右以及两个对角线共8个邻域比较。如果它严格大于所有8邻域就是局部极大值严格小于所有8邻域就是局部极小值。但图像里总有成片的平坦区比如天空、墙面此时会有大量像素同时“等于”邻域无法判定极大极小。常见解决方法是先用imregionalmax和imregionalmin找到极大/极小值区域再用bwmorph(BW,thin,Inf)把每个区域收缩成一个像素也可以求区域重心作为极值点坐标。这样得到的极值点数量可控包络也更平滑。包络面插值的选择直接决定BEMD分解质量。MATLAB里常见做法是F scatteredInterpolant(xs, ys, vals, natural, linear);其中xs、ys是极值点的横纵坐标vals是极值点的灰度值插值方法用natural自然邻域在平滑性和速度之间平衡。如果你追求更光滑的包络面可以自己写薄板样条TPS插值但TPS需要求解一个N*N的线性方程组极值点超过两三千就非常慢且容易内存不足。所以我在实际项目里经常先用“抽稀”策略把极值点数量压到2000以内再做natural插值。筛选停止条件和一维类似但二维SD的计算要对整幅图像逐像素做SD sum(sum((H_prev - H_current).^2)) / sum(sum(H_prev.^2))此外还有Cauchy型条件计算包络平均面的能量与原图能量的比值小于某个百分比就停止。我个人偏好用能量比值因为它更直观也不受边界像素突变影响。下面是二维筛选的一次迭代示意采用natural插值function H bemd_sift2d_demo(img, minDist, sdThresh) % 二维筛选一次迭代的近似示意 % 输入img为double灰度图像minDist为极值点最小间距sdThresh停止阈值 % 输出一个BIMF分量H I img; % 用imregionalmax/min找极值区域再提取中心点 maxReg imregionalmax(I); minReg imregionalmin(I); % 简化直接用find取坐标实际要处理区域收缩 [my, mx] find(maxReg); % 如果极值点太多按minDist抽稀 if ~isempty(my) length(my) 1 keep randperm(length(my), min(length(my), 1000)); % 示例抽稀 mx mx(keep); my my(keep); end valMax I(sub2ind(size(I), my, mx)); Fmax scatteredInterpolant(mx, my, valMax, natural, linear); [xx, yy] meshgrid(1:size(I,2), 1:size(I,1)); envUp Fmax(xx, yy); % 同理计算下包络envLo这里省略 M (envUp envLo) / 2; H I - M; % 停止判断这里简化为前后差异小于阈值 sd sum(sum((H - I).^2)) / sum(sum(I.^2)); if sd sdThresh return; end end注意这个演示代码里我用了randperm抽稀这是非常不可取的它会让分解结果不固定。实际应该用空间距离聚类或者网格化取局部最大值。我只是为了展示结构。后面避坑章节会专门说这个。真正的BEMD库不会这么写但很多粗糙的BEMD.zip里还真有随机抽稀这就是为什么同一幅图跑两次结果不一样。这个过程的复杂度主要在于插值。极值点数量是Nnatural插值构建三角剖分的时间接近O(N log N)之后在每个网格点求插值还要遍历邻域。所以一幅512x512的图像极值点如果不过滤很容易突破5000个一次筛选可能就要几十秒整幅分解几十个筛选循环跑下来10分钟很正常。2.3 BEMD和人工智能的关系当预处理特征提取器用还是当网络层用BEMD近几年在人工智能相关项目里出现多数是以“特征工程模块”身份。图像先做BEMD把得到的高频BIMF和残差当作多通道输入再接CNN或SVM。这样做的好处是BEMD能自适应地抓住图像中不同尺度的振荡模式这些模式比原始像素更有物理意义。比如在医学图像分类中BIMF1的纹理细节、残差的亮度分布分别对应不同诊断特征。也有少数人尝试把筛选过程改造成可微算子放进神经网络里。难点在于极值点搜索和插值梯度回传。极值搜索本质是max操作不可微插值倒是可以线性化。所以常见方案是固定极值点集合比如每层采样网格只让包络插值参与梯度计算。这写起来非常复杂而且BEMD本身迭代次数不固定网络训练时容易不稳。因此我一般不推荐把BEMD嵌入到网络内部而是离线算好特征存盘。在手头没有GPU的机器上BEMD传统分类器是一个性价比很高的方案。你可以先把训练图像库里的所有图像分解提取每层BIMF的统计量存成特征矩阵然后训练一个SVM。整个过程不需要深度网络代码量小还能作为“人工智能大作业”的亮点——把数值算法和AI结合比单纯调一个CNN有说服力。下表简单对比BEMD和常见分解方法方法基函数自适应能力空间局部性典型计算成本傅里叶变换固定正弦基无全局低小波变换固定小波基有限选小波多尺度中BEMD数据驱动包络完全自适应局部极值驱动高从这个表能看出来BEMD的最大特色就是完全自适应和局部性这正好是很多AI任务需要的。但成本高所以实际使用时要控制图像尺寸和极值点数量。3. 在MATLAB里跑通BEMD.zip最小复现代码与文件结构3.1 拿到压缩包先做三件事检查函数依赖、测试入口、看数据格式从网上下载的BEMD.zip打开后别急着运行先看看里面有什么。常见做法是解压后用MATLAB的dir命令遍历所有.m文件cd(D:\workspace\BEMD\); fileList dir(**/*.m); % 递归找所有m文件 for i 1:length(fileList) disp(fileList(i).name); end这段代码会列出压缩包里的所有MATLAB脚本和函数。你要关注三件事有没有demo或test脚本有没有主函数比如bemd2d.m或bemd_main.m有没有明显缺失的依赖。很多压缩包是把论文源码的一部分上传了少一个私有函数比如find_extrema.m或interp_surface.m直接跑主函数会报“Undefined function”。这时你先在整个文件列表里搜一下如果确实没有那就要去报错提示的调用链里逐个补。第二步是测试入口。如果能看到demo.m或example.m直接在MATLAB里运行它。但运行之前建议先看开头几行确认它读的数据文件是否存在。有时候demo里写的图片路径是./data/lena.bmp但压缩包里根本没有data文件夹。我会写一个快速检查文件存在的代码% 检查关键文件是否存在 if exist(demo.m, file) 2 disp(找到demo入口); end if exist(lena.bmp, file) 2 img imread(lena.bmp); disp(size(img)); else warning(缺少演示图像请先准备一张灰度图); end第三步是看数据格式。BEMD一般要求输入是double类型、灰度图。很多代码会默认图像范围是[0,1]如果直接传uint8的[0,255]插值结果会是两三百的样子后面计算能量就爆炸。我习惯统一做一次归一化if nargin 0 img imread(cameraman.tif); end if isinteger(img) img double(img) / double(intmax(class(img))); end3.2 最小运行代码读图、分解、显示BIMF假设我们手头的BEMD.zip里有个函数叫bemd2d它的接口大致是function [bimfs, residue] bemd2d(img, layers, sd, minDist)参数说明img是MxN的double灰度图像layers是分解层数sd是停止阈值minDist是极值点最小间距。这个接口不一定是你下载到的那个但比较通用。最小运行代码如下% 读取并归一化图像 img imread(cameraman.tif); img double(img) / 255; % 设置参数 layers 3; sd 0.3; minDist 3; % 调用分解 [bimfs, residue] bemd2d(img, layers, sd, minDist); % 可视化第一层和残差 figure(Name, BEMD分解结果); subplot(1,3,1); imshow(bimfs{1}); title(BIMF1); subplot(1,3,2); imshow(bimfs{2}); title(BIMF2); subplot(1,3,3); imshow(residue); title(残差);如果bimfs是三维数组比如bimfs(:,:,1)那么你必须按实际返回值方式调整。最好在调用前用which bemd2d和open bemd2d看一下函数定义。这段代码的逻辑说明归一化是必须的因为BEMD的插值过程依赖数值尺度如果输入是[0,255]的uint8包络面的拟合矩阵会有较大的条件数可能在筛选过程中产生NaN。而分解层数设为3是一般经验值第一层是高频细节第二层是中频结构第三层基本就是低频轮廓残差就是全局亮度趋势。如果你只做去噪往往只要BIMF1就够了做融合可能需要前两层做分类可能前三层加残差。3.3 把BIMF转成特征向量供后续聚类/分类使用对于人工智能项目BEMD分解后的BIMF不是直接拍平送给分类器的因为像素维度太高。常见做法是对每层计算一组统计量得到固定长度的特征向量。下面这段代码从bimfs元胞数组中提取特征% 输入bimfs为1*K元胞数组每层是M*N矩阵 numLayers length(bimfs); feat zeros(1, numLayers * 5); for i 1:numLayers b bimfs{i}; b mat2gray(b); % 缩放到0-1便于计算熵和互不相关 feat(i*5-4) mean2(b); feat(i*5-3) std2(b); feat(i*5-2) sum(b(:).^2) / numel(b); % 平均能量 feat(i*5-1) entropy(b); % 灰度共生矩阵对比度需要影像处理工具箱 glcm graycomatrix(b, Offset, [0 1; -1 1; -1 0; -1 -1]); stats graycoprops(glcm, {Contrast, Correlation}); feat(i*5) mean(stats.Contrast); end这段代码里有几个点要注意mat2gray是对数据进行线性归一化让熵函数能够正确统计灰度分布graycomatrix需要Image Processing Toolboxgraycoprops返回的是结构体数组所以用mean取平均。如果没装工具箱可以把熵和对比度换成简单的能量梯度gx diff(b, 1, 2); gy diff(b, 1, 1); feat(i*5-1) mean(gx(:).^2) mean(gy(:).^2);这样整个特征提取就不依赖任何工具箱了。然后按行堆叠所有样本的特征生成训练矩阵% 假设X_train是N*D矩阵N是样本数D是特征维度 X_train [X_train; feat];逻辑说明BIMF特征向量的长度等于分解层数乘以每层特征数如果分解层数波动比如不同图像筛出的层数不一样特征矩阵对不齐。解决办法是固定分解层数让BEMD总是输出前K层加残差不足时补零多余时丢弃。很多MATLAB图像处理大作业里大家就是用这种方式把非定长输出变成定长特征。4. BEMD的三个必调参数筛选次数、停止条件、插值点数4.1 筛选次数sifting iterations直接决定BIMF的“纯度”筛选次数是BEMD中最容易被忽视的参数。很多代码里它叫maxIter或numSift默认值从5到20不等。筛选次数越多每个BIMF越接近“局部对称平衡”的理想状态但同时也意味着计算量线性增加。更重要的是筛选次数过多会发生“过筛选”把噪声从残差中吸入BIMF1。实测经验对一般图像10次够用20次是极限如果超过30次BIMF1会全是椒盐状噪声。判断筛选次数是否合适的办法是看相邻两次筛选结果的差异。你可以在BEMD循环里打点或者自己写一个小脚本% 假设bemd2d内部允许回调iter或者你自己改写筛选循环 for iter [5 10 15 20 25] [C, ~] bemd2d(img, 3, 0.3, 3, iter); % 假设第五个参数是maxIter energy sum(C{1}(:).^2) / sum(img(:).^2); fprintf(iter%d, BIMF1能量占比%.4f\n, iter, energy); end如果energy在iter15之后突然升高说明过筛开始。取拐点前的值比较安全。注意这个前提是BEMD函数真的接受maxIter参数有些实现写死了你需要改代码。我通常保留这个接口因为这是超参搜索的第一步。4.2 停止条件SD阈值调不好就是欠筛或过筛SD阈值表示两次筛选之间的归一化差。阈值越小筛选越苛刻所需迭代次数越多。二维下由于像素数量大SD的数值容易偏小如果沿用一维的0.3往往还没筛完就停了导致BIMF1里混着大尺度结构。我一般对图像设阈值0.4~0.6。不过这也和图像对比度有关如果图像已经归一化到[0,1]阈值0.5比较稳。你可以用下面的脚本快速验证不同SD对残差的影响sds [0.2 0.3 0.4 0.5 0.6]; for i 1:length(sds) [~, res] bemd2d(img, 3, sds(i), 3); resEnergy sum(res(:).^2) / sum(img(:).^2); fprintf(SD%.1f, 残差能量占比%.3f\n, sds(i), resEnergy); end一般来说残差能量占比应该占原图的50%~80%左右。如果你看到SD0.5时残差能量高达95%说明分解深度远远不够——你可能需要增加分解层数而不是继续调小SD。如果SD0.3时残差能量只有20%说明过分解了许多低频信息都被抽到BIMF里去了。4.3 插值点数与边界延拓暗藏玄机的计算量黑洞插值点数指的是参与构造包络面的极值点个数。这个参数通常不直接暴露而是通过minDist或“最大极值点数”间接控制。我见过一个BEMD实现内部用find(imregionalmax(img))直接找极值一张512x512的图能找到6000多个极大值点。对这些点做scatteredInterpolant的natural插值一次包络面就要构造并解一个三角剖分内存开销几GB。如果你发现MATLAB内存飙升八成是极值点太多。设置minDist可以强制极值点之间至少相隔几个像素。在MATLAB里一种实现方式是对二值极大值图做形态学腐蚀% 限制极值点最小间距 maxReg imregionalmax(img); maxReg imdilate(maxReg, strel(disk, minDist)) maxReg; % 再用bwmorph提取中心 maxReg bwmorph(maxReg, shrink, inf); [ys, xs] find(maxReg);这样可以把相邻极值点合并成一个。注意imregionalmax会返回平台区所有像素为1bwmorph shrink会把每一块连续区域收缩为一个像素这是稳妥的做法。如果你想让极值点数量进一步减少可以对xs, ys做一个网格统计每个网格只留最大灰度极值点。边界延拓则是另一个决定结果质量的参数。BEMD在图像边缘处没有相邻像素包络面会在边缘处翘起导致分解出的BIMF在边缘有伪振荡。常见解决方法是在分解前对图像做对称延拓扩展宽度为8~16像素分解后裁掉。下面是一个延拓与裁剪的封装function [bimfs, res] bemd2d_with_extend(img, layers, sd, minDist, ext) % 对图像四周做镜像延拓然后调用bemd2d最后裁剪 extImg padarray(img, [ext ext], symmetric); [bimfs, res] bemd2d(extImg, layers, sd, minDist); for k 1:length(bimfs) bimfs{k} bimfs{k}(ext1:end-ext, ext1:end-ext); end res res(ext1:end-ext, ext1:end-ext); end这段代码里padarray的symmetric是镜像延拓比replicate重复边缘更自然能有效抑制边缘振铃。延拓宽度ext建议设为图像最长边的2%~4%比如512图取16。需要特别提醒有些BEMD源码内部使用griddata(v4)但MATLAB新版已经移除v4方法导致直接报错。遇到这种情况把griddata(x,y,v,xi,yi,v4)改成scatteredInterpolant(x,y,v,natural)后再插值输出会有微小不同但至少能跑。5. BEMD实战避坑MATLAB里最常见的5个翻车现场5.1 现象分解出的BIMF全是噪声没有纹理层你可能看到BIMF1是一幅几乎全是细碎斑点的图看不出任何边缘或纹理。原因通常是筛选次数过多或者插值方法不当。特别是当代码用了线性插值linear时包络面不光滑局部均值面到处是折痕筛出来的残余就是高频伪影。解决方法是改用natural插值同时把最大筛选次数降到10左右。如果图像本身已经有噪声先做一次轻度高斯滤波再分解比如img imgaussfilt(img,0.8)否则BEMD会把噪声当作最高频振荡模式独立筛出来。5.2 现象运行时间爆炸一幅512x512图像跑十分钟这是BEMD最劝退的问题。原因除了极值点太多外还可能是代码里用了双层循环遍历每个像素去判断极值而不是用imregionalmax向量化。解决思路是先剖开看看哪个函数耗时。在MATLAB中运行profile on bemd2d(img, 3, 0.3, 3); profile viewer如果scatteredInterpolant占了大头说明极值点数量失控。先把minDist设到4或5让极值点数量减少到原来的四分之一运行时间能大幅缩短。如果瓶颈是循环那只能改写向量化。另外千万不要在筛选循环内部用imshow或disp打印中间过程这会拖慢几十倍。5.3 现象边界处出现强烈的“振铃”或伪分量边界振铃是插值包络在图像边缘没有足够极值点约束造成的。常见出现在BIMF1的上下左右边缘有像是波浪的亮暗条纹。解决方法是先延拓再分解延拓后裁掉。如果延拓后还有振铃可以尝试对图像边缘做加权淡入淡出或者将图像分块分解后做重叠融合。但重叠融合比较复杂我的经验是先加镜像延拓90%的振铃都能被压下去。5.4 现象不同MATLAB版本结果不一致同一个BEMD.zip在MATLAB R2016a和R2024a上跑出来的BIMF不一样这是正常现象。原因在于griddata和scatteredInterpolant底层算法在不同版本中改动过三角剖分和插值结果会有细微差异。另一个原因是旧版randperm行为可能不同如果代码内部有随机抽样结果天然不稳定。解决方法是固定插值方法为natural固定随机种子并且把关键参数写进注释。如果要求可复现最好在同一版本下运行并在结果保存时记录MATLAB版本和参数。5.5 现象内存不足Out of Memory在二维筛选中高发当极值点数量达到几千甚至一万scatteredInterpolant创建时可能直接抛出“Requested array exceeds maximum array size preference”。原因是natural插值需要构建Delaunay三角剖分复杂度O(N log N)内存需求随N快速上涨。解决方法一是减少极值点数量上面说过二是用分块分解三是如果图像太大先imresize到256x256再分解。注意BEMD本身对分辨率很敏感强行缩放也会改变分解层次但在内存受限的机器上这是必要的妥协。6. 让BEMD真正工作起来图像融合场景的完整验证与进阶技巧最后这一章我拿一个具体场景把整篇文章串起来多聚焦图像融合。假设你有两张同一场景的照片一张左边清晰、右边模糊另一张右边清晰、左边模糊。经典做法是用拉普拉斯金字塔但我们这里用BEMD来验证分解结果有没有意义。融合规则很简单两幅图像分别做BEMD得到各自的前两层BIMF和残差。高频部分BIMF1和BIMF2用局部能量或梯度作为清晰度度量取清晰度大的像素作为融合结果残差部分直接平均因为没有“清晰不清晰”之分残差代表背景强度。实现代码如下% imgA、imgB是已配准的双灰度图尺寸相同 [C_A, r_A] bemd2d(imgA, 2, 0.4, 4); [C_B, r_B] bemd2d(imgB, 2, 0.4, 4); % 初始融合结果用残差平均 fused (r_A r_B) / 2; % 逐层融合高频BIMF for k 1:2 % 用水平Sobel梯度幅值代表清晰度 ga abs(imfilter(C_A{k}, fspecial(sobel))); gb abs(imfilter(C_B{k}, fspecial(sobel))); w ga gb; fused fused w .* C_A{k} (1-w) .* C_B{k}; end这段代码里我简化了清晰度衡量只用了水平Sobel。实际可以用imgradient算多方向梯度效果更稳ga imgradient(C_A{k}); gb imgradient(C_B{k});验证方式也很重要光看图片觉得“好像更清晰”不够。可以计算融合图像与原两张清晰部分的梯度能量对比或者用量化指标信息熵、标准差、平均梯度。MATLAB一行就能算metric std2(fused) entropy(mat2gray(fused)) mean2(abs(imfilter(fused, fspecial(sobel))));如果metric明显高于两张原图各自的对应值说明融合确实保留了更多细节。注意BEMD对两幅图分解时用的是同一组参数否则BIMF层数可能不同导致C_A{1}和C_B{1}尺度不匹配。如果出现这个问题强制设置layers2而且必须保证两幅图的尺寸完全一样。最后分享一个我的个人习惯每次跑BEMD之前我都会把输入图像缩放到最大边不超过512像素。不是算法不行而是工程上值得。BEMD的筛选过程是确定性的但计算代价随图像面积近似平方增长缩小四倍分辨率时间能快十倍以上。先在小图上把参数调好再放大到最终尺寸能少熬夜。如果你一开始就在1024图上慢慢等十有八九会等到怀疑人生。希望这个血泪经验能帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助?