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)存m+1维匹配数 for dim = m:(m+1) 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(i+1: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(m+1) / 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(i+1: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是否在合理范围，如果取m=3、r=0.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的熵值计算程序，前后迭代了大约三周时间。它帮我实现了一个完整的分析管线：从原始数据导入、滑窗熵值计算、统计检验到结果可视化。现在面对新的不熟悉的数据集，我也习惯先算个熵值曲线看看整体形态，这个习惯帮我在好几个项目里都迅速锁定了数据中的"特殊片段"，之后再针对这些片段做深入分析，效率提升非常明显。如果你也在做类似的状态监测或者复杂度分析，这个方案可以直接作为起点，然后根据你自己的数据特性再逐步调整参数和流程。