MATLAB熵值计算全攻略：从香农熵到样本熵的工程实践

搞数据分析的朋友应该都遇到过这种情况手里拿到一批复杂的观测数据想衡量它的混乱程度或者有序性但用均值、方差这些经典指标又觉得差点意思。我个人的经验是这种时候直接上信息熵往往能找到别的统计量挖不出来的规律。这篇文章就是围绕我自己在MATLAB里写熵值计算程序的全过程展开的从最基础的香农熵公式怎么落地成代码到后面接上统计分析、结果解读整套流程都会拆开讲清楚。如果你手头正好有需要做信号分析、状态识别或者系统复杂度评估的活儿不管是做研究还是搞工程这篇应该都能给你一个可以直接抄作业的参考方案。我用的环境是MATLAB R2021b核心的熵计算逻辑和统计绘图全部基于MATLAB原生函数没有引入额外的第三方工具箱。这么选纯粹是因为熵值计算的原理并不复杂纯手写反而能保证每一行逻辑都在自己掌控里出了问题也方便排查。1. 熵值计算的整体思路拆解1.1 为什么用熵来衡量数据很多人一听熵就联想到热力学觉得是个物理概念离数据分析很远。实际上在信息论里熵衡量的是随机变量的不确定性。简单说如果一个数据序列的取值非常规律、很容易预测那它的熵值就低反过来如果序列杂乱无章、看不出任何套路熵值就高。举个例子一个恒定的信号序列比如全是1它的熵是0因为没有任何不确定性而一个在每个取值上都均匀分布的序列熵就会达到该取值范围内的最大值。对工程场景来说熵值的变化往往能直接反映系统状态的变化——机械设备正常运转时振动信号可能有一定复杂度一旦出现故障信号可能变得更有规律熵值显著下降也可能变得异常混乱熵值飙升。所以熵值计算的核心价值在于它给数据的混乱程度提供了一个可量化、可比较的数值标尺。我最初做这个程序是想分析某组传感器采集到的时间序列数据判断它在不同工况阶段的状态差异。传统的时域特征比如峰值、均方根值虽然能看出幅值变化但要区分不同阶段的运行模式熵值的区分度明显更好。1.2 常见熵值类型与选型思路熵值计算其实是一个家族除了最基础的香农熵还有近似熵、样本熵、排列熵、多尺度熵等等。我最初的需求是先把最标准的香农熵实现好因为它是最根本的定义所有的其他熵基本上都是在它的思想基础上做改进。熵类型核心思想适用场景计算复杂度香农熵基于取值概率分布离散序列、分布分析低近似熵衡量序列中新模式出现的概率生理信号、时间序列复杂度中样本熵改进近似熵的自匹配偏差生理信号、短序列分析中高排列熵基于相邻值相对顺序时间序列快速复杂度评估低这篇文章里的程序我同时实现了香农熵和样本熵两个版本前者用来分析离散分类数据的分布均匀性后者用来分析连续时间序列的复杂度。如果你的实际需求只关心序列的可预测性和复杂度建议优先做样本熵如果只是分析离散取值的分布特性香农熵就足够没必要引入额外的计算开销。1.3 MATLAB作为实现工具的理由选MATLAB倒不是因为别的语言不能做。熵的计算本身用Python写也很方便但我在实际工作里更依赖MATLAB的原因有两点一是MATLAB天生面向矩阵运算数据导入导出、批量处理这类活儿非常顺手二是后期做统计分析、画频谱图、概率密度图MATLAB的绘图交互体验比写Python脚本调库要直观不少数据探索阶段的效率优势很明显。另一个实际考量是通过MATLAB的App Designer或者脚本封装可以快速把计算过程做成一个小工具团队成员不需要懂编程细节也能使用。这对于需要反复验证数据特征的场景帮助很大。2. 核心实现与参数详解2.1 香农熵的基础实现先看最基本的香农熵实现。假设有一组离散数据data它的取值范围是有限的若干个类别那么香农熵的计算很简单function H shannonEntropy(data) % 原始数据的香农熵计算 % data: 向量元素为离散类别标识数字或可直接比较的值 % H: 香农熵单位比特 % 核心思路估计每个取值类别的概率再求 -p*log2(p) 的累加 data data(~isnan(data) ~isinf(data)); if isempty(data) H NaN; return; end n numel(data); uniqueVals unique(data); % 统计各取值频数 counts zeros(size(uniqueVals)); for i 1:numel(uniqueVals) counts(i) sum(data uniqueVals(i)); end % 计算概率 p counts / n; % 去除0概率项因为 0*log(0) 在信息论中约定为0 p p(p 0); % 信息熵公式底数为2时单位是bit H -sum(p .* log2(p)); end这段代码的关键点我已经写在注释里了但有几个细节一定要强调第一必须先剔除NaN和Inf否则频率统计会出现莫名其妙的结果第二unique函数会默认排序这个不影响计数但是会稍微影响计算顺序对最终熵值没有影响第三也是最重要的0 * log(0)在数学上未定义但在信息论里约定其极限为0所以必须把零概率项过滤掉不然MATLAB会直接返回NaN。2.2 样本熵的实现及参数如何定香农熵虽然好用但它只看取值的概率分布完全忽略了先后顺序。对时间序列这种顺序本身携带信息的数据来说这显然是不够充分的。样本熵就是为了弥补这个缺陷它计算的是序列里生成新模式的概率。如果序列越是随机、越是难以预测样本熵就越大如果序列自我重复性强、有规律样本熵就更小。样本熵的实现代码我在项目里是这么写的function SampEn sampleEntropy(data, m, r) % 样本熵计算 % data: 一维连续时间序列 % m: 嵌入维数通常取2 % r: 相似容限通常取(0.1~0.25) * std(data) data data(:); data data(~isnan(data)); N numel(data); % 提前计算容限避免每次重复 tolerance r * std(data); % 符号模板 phi zeros(1, 2); % phi(1)存m维匹配数phi(2)存m1维匹配数 for dim m:(m1) patterns zeros(N - dim 1, dim); for i 1:(N - dim 1) patterns(i, :) data(i : i dim - 1); end % 用矩阵化方法计算任意两个模板之间的切比雪夫距离 % 这种方法比双重循环快很多 count 0; for i 1:(N - dim) % 计算当前模板与后续所有模板的距离 diffMat abs(patterns(i1:end, :) - patterns(i, :)); maxDiff max(diffMat, [], 2); count count sum(maxDiff tolerance); end phi(dim - m 1) count / ((N - dim) * (N - dim 1) / 2); end % 样本熵定义-ln(phi(m1) / phi(m)) SampEn -log(phi(2) / phi(1)); end这段代码里我对模板匹配的部分做了优化通过矩阵化计算距离替代逐元素的两两循环在序列长度为几千到几万时效率提升很明显。实际计算时有两个参数需要根据数据特征来定嵌入维数m常规工程数据取2居多取1时对波动不够敏感取3以上容易丢失细节。相似容限r我一般取0.15 * std(data)作为起点如果数据噪声较大就适当放宽到0.2倍标准差噪声小就收紧到0.1倍。这两个参数直接影响样本熵的绝对值所以做对比研究时所有数据序列必须用同一组参数不然结果之间没有可比性。这个可以说是踩过坑之后总结出的铁律。2.3 滑动窗口熵值的批处理实现单个序列算一个熵值只是第一步现实问题里往往需要看熵值随时间的变化趋势。我的做法是设计一个滑动窗口批处理函数把整个长序列按照指定窗口长度切成多个片段逐段计算熵值最终得到一条熵值曲线。function entropySeries slidingEntropy(data, windowSize, stepSize, entropyFunc, varargin) % 滑动窗口熵值计算 % data: 输入的长序列 % windowSize: 窗口长度 % stepSize: 每次滑动的步长 % entropyFunc: 函数句柄指向具体熵函数 N numel(data); nWindows floor((N - windowSize) / stepSize) 1; entropySeries zeros(nWindows, 1); startIdx 1; for i 1:nWindows windowData data(startIdx : startIdx windowSize - 1); entropySeries(i) entropyFunc(windowData, varargin{:}); startIdx startIdx stepSize; end end窗口长度的选择很重要。窗口太短统计样本不足熵值波动非常大趋势看不清窗口太长熵值变化滞后时间分辨率下降。我自己用下来对采样率1000Hz的传感器信号窗口长度取1000到2000个点即1到2秒效果比较理想步长取窗口长度的十分之一左右既能保证曲线平滑又不会损失太多细节。2.4 计算效率优化的个人建议刚开始写样本熵代码时我用的是三重循环数据量稍微一大就慢得让人怀疑人生。后来逐步优化总结了两条对提升效率最有用的经验第一能用矩阵运算就不写循环。在对模板距离的计算中用abs(patterns(i1:end, :) - patterns(i, :))替代内层循环代码短了一截计算速度快了好几个数量级。这一点在很多MATLAB算法优化场景里都适用。第二数据预处理能省则省。如果输入信号是模拟采集的连续值直接用原始值做概率统计意义不大建议先做等宽离散化或者取整这样既能加速计算还能提升熵值的稳定性。比如把浮点序列按区间切成若干个等级再计算香农熵运行效率和可解释性都会好很多。3. 统计分析与结果解读3.1 熵值结果如何描述与对比熵值算出来只是第一步怎么分析和解读这些数值才是决定项目成败的关键。我在项目里先对各个工况阶段的熵值做了基本描述性统计包括均值、标准差、中位数和四分位距stats struct; stats.meanVal mean(entropySeries); stats.stdVal std(entropySeries); stats.medianVal median(entropySeries); stats.iqrVal iqr(entropySeries);这里有个细节值得留意熵值分布的均值和中位数如果差异较大说明分布有偏态直接用均值做比较可能会掩盖重要的分组差异。我曾遇到过一组数据均值看着差别不大但中位数和箱线图排出来之后不同工况的差异其实非常明显。可视化方面我用箱线图来表达不同条件下的熵值分布最直观一组数据一个箱子上下四分位线、中位数、离群点全都一目了然。另外也会把熵值曲线和原始信号放在同一张图上对比观察熵值突变点对应的原始信号形态这种对照分析经常能带来额外发现。3.2 分组差异的统计检验方法做对比分析时不能只靠看起来有差异就下结论必须跑统计检验来确认差异的显著性。如果只比较两个组别我会用ttest2直接做两样本t检验前提是数据近似正态分布且方差齐性比较明显[h, p] ttest2(entropySeries_A, entropySeries_B, Alpha, 0.05);如果组别多于两个比如有正常、过渡、异常三种工况就得用单因素方差分析[p, tbl, stats] anova1([entropyGroup1; entropyGroup2; entropyGroup3], ... [groupLabel1; groupLabel2; groupLabel3], off);方差分析的p值只说明各组之间至少存在一组差异具体是哪两组之间有差异需要进一步做多重比较c multcompare(stats);多重比较结果里两组的置信区间如果不包含0就意味着这两组间差异显著。做多组比较时建议用multcompare默认的Tukey法它对多重比较带来的第一类错误控制得比较稳妥比直接做多次t检验更可靠。3.3 熵值与其他特征的相关性分析有时候光看熵值本身还不够我会额外算熵值与原始信号其他特征量比如均方根值、峰值因子、功率谱特征之间的相关系数。这里推荐用Spearman相关系数而不是Pearson因为很多情况下熵值与其他特征之间不是严格的线性关系Spearman基于秩相关对非线性单调关系也能捕捉到一定的趋势[rho, pval] corr(entropyVals, rmsVals, Type, Spearman);这个分析的意义在于如果发现熵值跟某个传统特征高度相关后续工程里可以直接用更容易计算的传统特征作为代理指标如果相关性很弱则说明熵值提供了独立信息具有不可替代性。3.4 评估不同状态识别效果的参考思路如果把熵值当作一个分类特征来用比如判断当前系统处于哪种状态我建议做一步简单的交叉验证来评估熵值的区分能力。我自己通常按如下流程走把样本按时间段切分为训练集和测试集保证同一状态的数据不跨集合重叠。在训练集上计算各状态熵值的统计分布划定简单的判定阈值。在测试集上验证阈值判别的准确率。这么做虽然比较粗糙但能在项目初期快速判断熵值是否有实用价值够用了。如果后续发现区分度不够理想再往两个方向优化一是调整熵计算的参数换窗口长度、改嵌入维数等二是把熵值和传统特征拼接起来组成特征向量用线性判别分析或者简单的机器学习分类器去分类。在某个模拟数据项目里我前期只用香农熵的区分准确率大约在78%后来叠加了样本熵和均方根值两个特征通过线性分类器直接提升到了92%。这说明熵值往往不是单独起作用的它的价值更多时候体现在和别的特征搭配时。4. 常见问题与排查实录4.1 数据预处理导致的零概率与NaN问题熵值计算程序在实际数据上跑最容易碰到的就是数据里混着NaN、Inf和重复异常点。前两类不处理概率统计会收到污染重复异常点如果过多会让某个取值的概率异常偏高熵值被严重拉低。我在程序里用了两种策略组合清理数据先执行data data(~isnan(data) ~isinf(data))把明显的异常值剔除再做一步去重中的过采样判断即如果某个类别取值数量超过整体样本的95%直接从序列中剃掉一部分避免某个异常工况主导整个熵值计算结果。另外一个新手容易踩的坑是直接对float序列的原始值做unique由于浮点数精度问题几乎每个值都不一样导致熵值直接逼近最大值完全失去区分度。解决方法是先用round或者直方图分箱把连续值离散化再做频率统计。在信号分析场景一般把幅值归一化到0到1之间再均分为10到20个箱子效果就很好。4.2 样本熵计算结果不稳定的参数原因样本熵对参数m和r的敏感度非常高这是所有熵类算法的通病不是代码bug。我试过同一段数据r从0.1倍标准差改成0.25倍标准差熵值能差出三四倍。这个现象不是算法本身的缺陷而是相似容限改变了什么叫相似的判定标准。如果发现结果波动剧烈排查顺序是这样的第一步确认数据是否做过去趋势处理带有线性漂移的信号会让模板匹配结果失真第二步确认序列长度是否足够样本熵在短序列少于一两百个点上的统计可靠性较差直接不推荐使用第三步确认参数m和r是否在合理范围如果取m3、r0.05倍标准差大多数实际数据都算不出稳定结果。我自己的经验是在正式批量计算之前先随机挑几段不同工况的代表性数据画一个熵值随r变化的曲线选曲线中相对平稳的区间来定参数这样能显著提升后续结果的可重复性。4.3 大批量数据计算速度过慢批量处理成百上千个文件时MATLAB的运行速度会成为瓶颈。我遇到过连续跑两小时都没算完一半的情况后来做了三件事效果立竿见影第一把循环里的日志输出和中间变量赋值全部去掉只保留计算结果。第二对原始信号做降采样从1000Hz降到500Hz熵值的相对趋势基本不变但计算量减半。第三用parfor替代普通的for循环在有多核处理器的机器上并行处理不同文件。parpool(local, 4); % 开启4个并行工作线程 parfor i 1:numel(fileList) data loadData(fileList{i}); ent(i) sampleEntropy(data, 2, 0.15); end这里要注意样本熵函数内部如果有随机数或者写文件操作parfor会引发不可控的并行冲突。好在我这个计算是纯函数每个文件的计算相互独立并行化改造很顺利。另外如果只是算香农熵数据量百万级别以内MATLAB的原生逻辑运算就很快没必要并行。4.4 结果解读中的常见误区在使用熵值做判断时有一个非常关键但容易忽视的问题熵值高不代表状态好熵值低也不代表状态差。它只描述数据的复杂性和可预测性具体的好坏必须结合领域背景来判断。我见过有人把熵值变化直接等同于故障程度的案例这是不够严谨的。比如系统从正常运行切换到一个新的稳定工况熵值可能降低但这并不意味着设备状态劣化反过来如果没有预判到工况切换点熵值突然升高了也可能被误判为故障其实只是外部输入信号本身的复杂度在增加。所以我在做结果评估时一定会先列出所有可能影响熵值的外部条件比如工况切换、环境噪声变化、采样频率调整等等逐项确认后再把熵值变化与内部状态关联起来。就我踩过坑之后的理解熵值是一个现象描述工具它在状态监测里的定位是指标不是诊断结论。5. 实操总结与体会如果把这篇内容的所有代码和思路浓缩成可复用的实操建议我自己的核心体会大概有这几点开始写熵值计算程序之前先想清楚自己的数据是离散序列还是连续时间序列这决定了用香农熵还是样本熵选错方向后面做的全是无用功。参数一旦选定全部数据都要用同一组参数否则算出来的数值之间没有可比性。预处理优先级高于算法本身数据干净程度直接决定了熵值的可信度。算完熵值不要急着下结论先做统计检验再用可视化图核对最后结合领域知识做解释。这套基于MATLAB的熵值计算程序前后迭代了大约三周时间。它帮我实现了一个完整的分析管线从原始数据导入、滑窗熵值计算、统计检验到结果可视化。现在面对新的不熟悉的数据集我也习惯先算个熵值曲线看看整体形态这个习惯帮我在好几个项目里都迅速锁定了数据中的特殊片段之后再针对这些片段做深入分析效率提升非常明显。如果你也在做类似的状态监测或者复杂度分析这个方案可以直接作为起点然后根据你自己的数据特性再逐步调整参数和流程。