简介本资源是一份面向信号处理与非线性时间序列分析初学者及进阶开发者的MATLAB实用工具包聚焦于排列熵Permutation Entropy这一经典复杂度测度算法的完整实现。适用于工控故障诊断、脑电/心电信号分析、振动监测等需量化序列随机性与规律性的实际场景。压缩包仅含1个核心文件——permutations_entropy.m函数脚本代码结构清晰、注释详尽涵盖相空间重构、符号序列生成、概率分布计算及熵值输出全流程无需额外依赖即可直接运行验证。文件大小仅1KB轻量高效便于嵌入现有MATLAB项目或用于教学演示。目前已有702人下载学习是理解排列熵原理、快速复现算法逻辑、开展特征提取研究的可靠参考源码。1. 项目缘起从“混沌”中寻找秩序的信号处理利器最近在整理一个关于轴承故障诊断的老项目时我又翻出了排列熵Permutation Entropy, PE这个老朋友。当时为了从一堆嘈杂的振动信号里揪出早期故障的蛛丝马迹试过不少时频分析和复杂度指标排列熵是其中让我印象最深的一个。它不像小波变换那样需要复杂的基函数选择也不像近似熵、样本熵那样对参数极其敏感。它的核心思想异常简洁优雅不关心信号的具体幅值只关注数据点之间随时间变化的“序”的模式。这种对序列“形态”而非“数值”的专注让它对噪声有天然的鲁棒性计算效率也高特别适合处理工程现场采集的、信噪比不高的非平稳时间序列。简单来说你可以把一段信号想象成一条起伏的曲线。排列熵不关心这座“山”具体有多高它只关心你沿着这条曲线行走时眼前看到的“上坡、下坡、平台”这种局部地形变化的模式是否丰富、是否可预测。如果模式单一、高度重复比如完美的正弦波熵值就很低意味着序列规律性强、复杂度低如果模式杂乱无章、毫无规律比如白噪声熵值就很高意味着序列随机性强、复杂度高。这种特性使得排列熵在机械故障诊断、生物医学信号分析如EEG、ECG、金融时间序列分析乃至气候数据研究中都成了香饽饽。网上能找到的排列熵代码不少但很多是“黑箱”实现参数含义模糊边界条件处理粗糙直接拿来用容易踩坑。所以我决定结合自己多年的使用和调试经验把核心算法、参数选择的门道、实际应用中的注意事项以及一个经过大量实测、鲁棒性更好的MATLAB程序源码整理出来。这份源码不仅仅是算法的简单翻译它包含了数据预处理、多尺度扩展MSPE、以及结果可视化的完整流程希望能帮你绕过我当年走过的弯路直接把这个强大的工具用起来。2. 排列熵的核心原理为什么“排序”比“数值”更聪明要理解排列熵为什么有效得先抛开复杂的公式看看它到底在计算什么。给定一个时间序列{x(1), x(2), ..., x(N)}排列熵的计算可以分解为以下几步2.1 相空间重构从一维序列到多维视图这是几乎所有非线性时间序列分析的第一步目的是还原系统潜在的动力学结构。对于排列熵我们采用最常用的时间延迟嵌入法。嵌入维度 (m) 这是最重要的参数决定了我们每次观察多少个连续的数据点。你可以把它理解为“观察窗口”的长度。例如m3就意味着我们每次截取3个连续的点作为一个观察单元。时间延迟 (τ) 通常取1即连续采样。对于某些特定周期的信号可以取其他值以捕获不同时间尺度上的动力学但在绝大多数通用场景下τ1是默认且安全的选择。重构后我们得到一系列m维的向量X(i) [x(i), x(iτ), ..., x(i(m-1)τ)], 其中i 1, 2, ..., N-(m-1)τ。 每个X(i)都是原始信号的一个“快照”。2.2 序数模式编码将向量转化为“排名符号”这是排列熵的灵魂所在。对于每一个m维向量X(i)我们不再关心它的具体数值而是关心其内部元素的大小排列顺序。例如假设m3我们有一个向量[2.1, 5.6, 1.8]。我们先将这三个数排序1.8 2.1 5.6。排序后它们的原始索引是[3, 1, 2]因为最小值1.8在原数组第3位次小值2.1在第1位最大值5.6在第2位。这个索引序列(3,1,2)就定义了一种序数模式Ordinal Pattern。对于m3所有可能的序数模式是3! 6种即(1,2,3), (1,3,2), (2,1,3), (2,3,1), (3,1,2), (3,2,1)。每一种模式都代表了数据点之间一种特定的局部形态关系比如(1,2,3)代表单调递增(3,2,1)代表单调递减(2,1,3)代表先降后升的“V”形。注意如果向量中存在相等的值平局标准的排列熵算法通常采用其出现顺序的索引但更严谨的做法是给它们分配相同的排名例如使用tiedrank函数。在实际生理或振动信号中由于是连续值完全相等的概率极低但量化后的数据可能出现。我提供的源码中包含了处理这种情况的稳健逻辑。2.3 统计与熵值计算衡量模式的“混乱度”经过上一步我们把长度为N的时间序列转化为了一个由K N - (m-1)个序数模式符号组成的序列每个符号是1到m的一个排列。接下来我们统计这m!种可能的序数模式在整段序列中出现的频率。假设第j种模式出现了C_j次那么其相对频率为p_j C_j / K。最后套用香农熵的公式来计算排列熵H_p(m)H_p(m) - Σ (p_j * log2(p_j))求和范围是所有p_j 0的模式。熵值H_p(m)的范围是[0, log2(m!)]。为了便于在不同m值下比较通常会进行归一化处理H_p_normalized H_p(m) / log2(m!)此时值域为[0, 1]。为什么这样做是聪明的抗噪性强小幅度的噪声可能会改变信号的幅值但不太可能颠覆多个连续点之间的大小顺序关系。因此排列熵对加性噪声不敏感。计算快速核心操作是排序和计数算法复杂度约为O(N*m*log(m))比需要计算距离矩阵的近似熵快得多。概念直观它直接度量了序列局部形态的多样性和不可预测性物理意义明确。3. 关键参数选择让排列熵在你的数据上发挥最大效力排列熵用起来简单但想把它的威力完全发挥出来参数选择是关键。这里没有放之四海而皆准的“黄金值”必须结合你的数据特性来定。3.1 嵌入维度m决定观察的细致程度m的选择是核心中的核心。它直接决定了你能分辨多少种不同的局部模式。m太小如2或3模式种类太少2!2,3!6分辨率不足可能无法区分一些动力学行为相似的复杂系统。例如一个4周期振荡和一个8周期振荡在m2下可能表现出相似的熵值。m太大如7或8模式种类急剧增加7!5040,8!40320。对于有限长度的数据N会导致绝大多数模式只出现一次或根本不出现使得计算出的概率分布不可靠熵值估计产生较大误差。同时计算量也会增加。经验法则Bandt 和 Pompe 在提出该方法的原始论文中建议m在3到7之间。我的实践经验是对于大多数工程和生物信号数据长度N在1000到10000量级m4或m5是一个非常好的起点能在区分度和统计可靠性之间取得很好的平衡。你可以做一个简单的验证逐步增加m观察归一化排列熵值的变化。当m增加到某个值后熵值趋于稳定或开始不规则波动那么前一个m值通常是合适的。3.2 时间延迟τ通常就选1除非你有很强的先验知识知道你的信号在某个特定时间尺度上有重要特征例如心电信号的R-R间期否则τ1即使用连续采样点是最通用、最不容易出错的选择。选择τ1相当于对信号进行下采样可能会丢失高频细节信息。3.3 数据长度N越长越稳但有下限为了保证统计的可靠性序列长度N需要远大于可能出现的模式种类m!。一个常用的经验法则是N 5 * m!。例如当m5时m! 120那么N最好大于600。对于m6(720)则要求N 3600。这也是为什么通常不推荐使用m7的原因——对数据量的要求太高。我提供的源码中包含了数据长度检查的警告提示。3.4 归一化为了跨尺度比较如前所述使用归一化的排列熵(H_p / log2(m!))几乎总是更好的选择。这使得不同m值下计算出的熵值可以放在同一个尺度[0,1]内进行比较结果更直观。4. 实战进阶多尺度排列熵与MATLAB源码深度解析基础的排列熵只能反映单一时间尺度上的复杂度。然而很多真实世界的信号如机械振动、脑电波其复杂性体现在多个时间尺度上。为此Costa等人提出了多尺度排列熵。4.1 多尺度排列熵的计算流程MSPE的核心思想是先对原始序列进行粗粒化得到不同尺度因子下的新序列再对每个新序列计算排列熵。粗粒化过程对于尺度因子s将原始序列{x1, x2, ..., xN}划分为长度为s的不重叠窗口计算每个窗口内数据的均值形成粗粒化序列{y_j^(s)}。y_j^(s) (1/s) * Σ_{i(j-1)s1}^{js} x_i,j 1, 2, ..., floor(N/s)当s1时就是原始序列。计算排列熵对每一个尺度因子s下的粗粒化序列y^(s)计算其排列熵PE(s)。绘制MSPE曲线以尺度因子s为横坐标PE(s)为纵坐标即可得到一条多尺度复杂度曲线。健康的、适应性强的系统如正常人的心率通常在不同尺度上表现出较高的熵值而病态的、僵化的系统如心衰患者的心率其熵值可能在多个尺度上降低。4.2 MATLAB源码实现与关键技巧以下是我提供的PermutationEntropy.m函数的核心部分解析并附上关键技巧说明。function [pe, patterns] PermutationEntropy(data, m, tau, normalize) % 计算时间序列的排列熵 % 输入 % data - 一维时间序列 (行向量或列向量) % m - 嵌入维度 (默认 4) % tau - 时间延迟 (默认 1) % normalize - 是否归一化到[0,1] (默认 true) % 输出 % pe - 排列熵值 % patterns - (可选) 所有序数模式的分布直方图数据 % % 作者基于Bandt Pompe (2002) 算法结合工程实践优化 if nargin 4 || isempty(normalize), normalize true; end if nargin 3 || isempty(tau), tau 1; end if nargin 2 || isempty(m), m 4; end data data(:); % 确保是列向量 N length(data); % 参数有效性检查 if m 2 error(嵌入维度 m 必须大于等于2.); end if tau 1 error(时间延迟 tau 必须为正整数.); end if N m * tau error(数据长度 N 必须大于 m * tau.); end % 数据长度建议检查 if N 5 * factorial(m) warning(数据长度可能不足以保证统计可靠性。建议 N 5*m!。当前 N%d, m!%d, N, factorial(m)); end % 1. 相空间重构 numVec N - (m-1)*tau; % 重构向量的数量 emb zeros(numVec, m); for i 1:m emb(:, i) data((1:numVec) (i-1)*tau); end % 2. 生成序数模式并统计 % 使用 sort 函数获取排序索引效率高且能处理平局情况 [~, idx] sort(emb, 2); % 沿行排序idx 的每一行就是一个序数模式 % 将索引模式转换为唯一的整数标签便于快速统计 % 技巧利用阶乘进制数将排列映射为唯一整数 factVec factorial(m-1:-1:0); patternLabels sum((idx - 1) .* factVec, 2) 1; % 统计每种标签出现的次数 possiblePatterns factorial(m); count histcounts(patternLabels, 1:possiblePatterns1); % 3. 计算概率和熵值 prob count(count 0) / numVec; % 只考虑出现的模式 pe -sum(prob .* log2(prob)); % 4. 归一化 if normalize pe pe / log2(possiblePatterns); end % 5. (可选) 输出模式分布 if nargout 1 patterns count; end end关键技巧与避坑点高效的模式编码代码中没有使用循环去比较每个模式而是利用sort函数一次性得到所有向量的排序索引然后通过阶乘进制数的技巧将每个索引序列如[2,3,1]映射为一个唯一的整数标签。这比用字符串比较或循环查找快了几个数量级尤其是当N很大时。平局处理MATLAB的sort函数在默认情况下是稳定排序即对于相等的值会保持它们在原始数据中的相对顺序。这提供了一种处理平局理论上极少但量化数据可能出现的合理方式。如果你希望更严格地处理平局如视为同排名可以使用tiedrank函数但会显著增加计算量。数据长度警告函数内置了长度检查警告。这是非常重要的调试辅助。如果你发现熵值不稳定或者不同数据段结果差异巨大首先应该怀疑数据长度是否足够。向量化操作整个重构和排序过程完全向量化没有显式循环充分利用了MATLAB的矩阵运算优势速度极快。4.3 多尺度排列熵的封装实现基于上面的核心函数实现多尺度版本就很简单了。以下是MultiscalePermutationEntropy.m的框架function [mspe, scales] MultiscalePermutationEntropy(data, m, tau, maxScale, normalize) % 计算多尺度排列熵 % 输入 % data - 一维时间序列 % m, tau, normalize - 同 PermutationEntropy 函数 % maxScale - 最大尺度因子 (默认 20) % 输出 % mspe - 各尺度下的排列熵值 (向量) % scales - 尺度因子 (向量) if nargin 5, normalize true; end if nargin 4, maxScale 20; end if nargin 3, tau 1; end if nargin 2, m 4; end data data(:); N length(data); scales 1:maxScale; mspe zeros(size(scales)); for i 1:length(scales) s scales(i); % 粗粒化 if s 1 coarseData data; else numCoarse floor(N / s); % 使用 reshape 和 mean 进行快速粗粒化 (需确保能整除) if mod(N, s) 0 coarseData mean(reshape(data(1:numCoarse*s), s, numCoarse), 1); else % 如果不能整除采用循环较慢或截断尾部 % 这里采用截断保证计算速度并给出提示 coarseData mean(reshape(data(1:numCoarse*s), s, numCoarse), 1); if i 1 % 只提示一次 warning(数据长度不能被所有尺度整除尾部数据将被截断用于粗粒化。); end end end % 计算当前尺度下的排列熵 mspe(i) PermutationEntropy(coarseData, m, tau, normalize); end end避坑点粗粒化时的长度处理粗粒化要求数据长度N最好是尺度因子s的整数倍。如果不能整除常见的处理方法是1) 截断尾部多余数据2) 使用滑动平均重叠窗口。我的代码采用了第一种方法截断因为它计算更快且对于尺度分析轻微的尾部截断通常影响不大。但在严谨的研究中你需要明确说明处理方法。第二种方法滑动窗口会引入数据间的相关性可能影响熵值的解释。5. 应用实例轴承故障振动信号分析让我们用一个实际的例子看看如何用这套代码分析轴承振动信号区分正常状态和早期故障。5.1 数据准备与预处理假设我们有两段振动加速度信号采样频率fs 12000 Hz时长2秒即每段N 24000点。一段来自正常轴承 (data_normal)一段来自外圈有轻微故障的轴承 (data_fault)。数据通常会有直流偏移或趋势项需要先去除。% 1. 数据预处理去趋势可选带通滤波 data_normal detrend(data_normal); % 去除线性趋势 data_fault detrend(data_fault); % 可选根据轴承特征频率设计一个带通滤波器保留感兴趣频段 % [b, a] butter(4, [100, 2000]/(fs/2), bandpass); % data_normal filtfilt(b, a, data_normal); % data_fault filtfilt(b, a, data_fault);5.2 计算单一尺度排列熵首先我们用默认参数快速计算一下两个信号的排列熵看看能否直接区分。m 4; tau 1; pe_normal PermutationEntropy(data_normal, m, tau, true); pe_fault PermutationEntropy(data_fault, m, tau, true); fprintf(正常信号排列熵 (m%d): %.4f\n, m, pe_normal); fprintf(故障信号排列熵 (m%d): %.4f\n, m, pe_fault);结果解读通常早期故障会引入周期性冲击使信号的规律性增强局部模式变得不那么“混乱”因此排列熵值会降低。所以我们很可能看到pe_fault pe_normal。但这并非绝对取决于故障类型和位置有时熵值也可能升高。单一尺度的结果只是一个初步的、粗略的指标。5.3 计算多尺度排列熵并可视化单一尺度可能不够稳健我们计算多尺度排列熵来获取更全面的特征。maxScale 30; [mspe_normal, scales] MultiscalePermutationEntropy(data_normal, m, tau, maxScale, true); [mspe_fault, ~] MultiscalePermutationEntropy(data_fault, m, tau, maxScale, true); % 可视化 figure(Position, [100, 100, 800, 400]); subplot(1,2,1); plot(scales, mspe_normal, b-o, LineWidth, 1.5, MarkerSize, 6, DisplayName, 正常); hold on; plot(scales, mspe_fault, r-s, LineWidth, 1.5, MarkerSize, 6, DisplayName, 故障); xlabel(尺度因子 s); ylabel(归一化排列熵 PE(s)); title(多尺度排列熵曲线); legend(Location, best); grid on; subplot(1,2,2); % 计算熵值曲线下的面积 (AUC) 作为一个综合特征 auc_normal trapz(scales, mspe_normal); auc_fault trapz(scales, mspe_fault); bar([1,2], [auc_normal, auc_fault]); set(gca, XTickLabel, {正常, 故障}); ylabel(MSPE曲线下面积 (AUC)); title(综合特征对比); grid on;分析与解读曲线形态正常轴承的MSPE曲线可能在多个尺度上都维持在相对较高的水平表明其动力学行为在不同时间尺度上都具有一定的复杂性。而故障轴承的曲线可能在低尺度反映细节或高尺度反映整体趋势上出现明显的下降表明故障冲击破坏了系统原有的复杂动力学结构使其行为变得更规则、更可预测。综合特征直接比较整条曲线有时不够直观。我们可以提取曲线的统计特征作为故障诊断的输入特征例如曲线下面积如上图所示AUC是一个很好的综合指标故障状态下AUC通常会减小。均值、方差计算所有尺度上熵值的均值和标准差。斜率计算曲线在低尺度段和高尺度段的平均斜率故障可能改变尺度间的复杂度传递特性。参数敏感性测试在实际应用中我强烈建议你对关键参数m进行敏感性分析。固定其他参数让m从3变化到6分别计算MSPE。观察不同m下正常与故障状态的区分度是否稳定。如果某个m值下区分度最好且稳定那么就选用该值。5.4 模式分布的可视化洞察动力学变化除了熵值这个汇总统计量观察序数模式本身的分布变化也能提供深刻见解。% 计算并可视化模式分布 [~, patterns_normal] PermutationEntropy(data_normal, 4, 1, false); [~, patterns_fault] PermutationEntropy(data_fault, 4, 1, false); figure; subplot(2,1,1); bar(patterns_normal / sum(patterns_normal)); xlabel(序数模式索引 (m4, 共24种)); ylabel(相对频率); title(正常轴承信号序数模式分布); xlim([0.5, 24.5]); grid on; subplot(2,1,2); bar(patterns_fault / sum(patterns_fault)); xlabel(序数模式索引 (m4, 共24种)); ylabel(相对频率); title(故障轴承信号序数模式分布); xlim([0.5, 24.5]); grid on;解读在m4时有24种可能的模式。健康系统的模式分布通常相对“平坦”即没有哪种模式占据绝对主导。而故障系统由于周期性冲击的引入可能会使某几种特定的模式例如对应冲击上升沿或下降沿的形态出现频率显著增高导致分布出现尖峰。这种分布的变化是熵值降低的根本原因也为理解故障的物理机制提供了线索。6. 常见问题、局限性与扩展方向即使有了好用的工具明白它的边界同样重要。6.1 常见踩坑点与解决方案熵值结果不稳定每次跑略有差异可能原因数据长度N严重不足。请检查N是否满足N 5*m!。对于短数据排列熵估计的方差会很大。解决方案增加数据长度或使用更小的m如从5降到4。也可以考虑对长数据分多个不重叠的段计算熵值然后取平均和标准差作为最终结果和置信区间。故障信号的熵值反而比正常信号高可能原因这并非不可能。如果故障引入了强烈的随机噪声而非周期性冲击或者使系统进入混沌状态熵值可能会升高。此外参数选择不当如m太大导致正常信号熵值被低估也可能造成反直觉结果。解决方案首先检查参数尤其是m是否合适。其次结合其他诊断方法如频谱分析、包络谱综合判断。排列熵只是一个特征不是万能判据。计算速度慢处理长数据时卡顿可能原因使用了低效的循环实现模式匹配。解决方案使用我提供的基于向量化和整数编码的算法速度会有数量级提升。对于超长序列如 10^6 点可以考虑先对数据分段计算或进行适当的降采样。对脉冲类故障敏感但对磨损类缓慢退化不敏感问题本质排列熵对序列的局部形态变化敏感。明显的冲击脉冲会改变局部形态分布。而均匀磨损可能在整个时间尺度上缓慢地改变信号能量但对局部点之间的序关系影响较小。解决方案将排列熵与反映能量变化的特征如RMS、峰峰值、小波能量熵结合使用构建多维特征向量。6.2 方法局限性丢失幅值信息这是其优点抗噪也是缺点。两个幅值分布完全不同但局部排序模式相似的序列会有相同的排列熵。因此它需要与其他时域、频域特征互补。对慢变趋势敏感一个强烈的线性趋势会主导局部排序导致模式分布严重偏向某几种如单调递增从而低估序列的真实动力学复杂度。因此在计算前对数据去趋势detrend是至关重要的预处理步骤。尺度粗粒化的信息损失MSPE中的粗粒化取平均是一个低通滤波过程会损失高频细节。有学者提出了复合粗粒化、层次熵等改进方案来弥补这一点。6.3 扩展与变体加权排列熵在计算模式概率时不仅考虑模式是否出现还考虑产生该模式的向量中数据点之间的实际数值差异。给差异大的模式赋予更大的权重从而在保留排序信息的同时引入部分幅值信息。模糊排列熵引入模糊隶属度的概念缓和了传统排列熵中“非此即彼”的硬判决使模式划分更平滑对噪声和参数选择的鲁棒性更强。多变量排列熵将算法扩展到多个相关时间序列用于分析系统不同变量间的耦合复杂度。在我提供的完整源码包MATLAB实现排列熵 程序源码.zip中除了上述核心函数还包含了用于批量处理数据、生成特征矩阵、以及绘制出版级质量图形的辅助脚本和示例数据。真正上手时建议你从示例开始替换成自己的数据调整参数观察变化这样才能最深切地体会这个工具的妙处和边界。信号分析没有银弹但排列熵无疑是你工具箱里一把锋利而独特的瑞士军刀。本文还有配套的精品资源点击获取 排列熵原理与MATLAB实现：时间序列复杂度分析利器
简介本资源是一份面向信号处理与非线性时间序列分析初学者及进阶开发者的MATLAB实用工具包聚焦于排列熵Permutation Entropy这一经典复杂度测度算法的完整实现。适用于工控故障诊断、脑电/心电信号分析、振动监测等需量化序列随机性与规律性的实际场景。压缩包仅含1个核心文件——permutations_entropy.m函数脚本代码结构清晰、注释详尽涵盖相空间重构、符号序列生成、概率分布计算及熵值输出全流程无需额外依赖即可直接运行验证。文件大小仅1KB轻量高效便于嵌入现有MATLAB项目或用于教学演示。目前已有702人下载学习是理解排列熵原理、快速复现算法逻辑、开展特征提取研究的可靠参考源码。1. 项目缘起从“混沌”中寻找秩序的信号处理利器最近在整理一个关于轴承故障诊断的老项目时我又翻出了排列熵Permutation Entropy, PE这个老朋友。当时为了从一堆嘈杂的振动信号里揪出早期故障的蛛丝马迹试过不少时频分析和复杂度指标排列熵是其中让我印象最深的一个。它不像小波变换那样需要复杂的基函数选择也不像近似熵、样本熵那样对参数极其敏感。它的核心思想异常简洁优雅不关心信号的具体幅值只关注数据点之间随时间变化的“序”的模式。这种对序列“形态”而非“数值”的专注让它对噪声有天然的鲁棒性计算效率也高特别适合处理工程现场采集的、信噪比不高的非平稳时间序列。简单来说你可以把一段信号想象成一条起伏的曲线。排列熵不关心这座“山”具体有多高它只关心你沿着这条曲线行走时眼前看到的“上坡、下坡、平台”这种局部地形变化的模式是否丰富、是否可预测。如果模式单一、高度重复比如完美的正弦波熵值就很低意味着序列规律性强、复杂度低如果模式杂乱无章、毫无规律比如白噪声熵值就很高意味着序列随机性强、复杂度高。这种特性使得排列熵在机械故障诊断、生物医学信号分析如EEG、ECG、金融时间序列分析乃至气候数据研究中都成了香饽饽。网上能找到的排列熵代码不少但很多是“黑箱”实现参数含义模糊边界条件处理粗糙直接拿来用容易踩坑。所以我决定结合自己多年的使用和调试经验把核心算法、参数选择的门道、实际应用中的注意事项以及一个经过大量实测、鲁棒性更好的MATLAB程序源码整理出来。这份源码不仅仅是算法的简单翻译它包含了数据预处理、多尺度扩展MSPE、以及结果可视化的完整流程希望能帮你绕过我当年走过的弯路直接把这个强大的工具用起来。2. 排列熵的核心原理为什么“排序”比“数值”更聪明要理解排列熵为什么有效得先抛开复杂的公式看看它到底在计算什么。给定一个时间序列{x(1), x(2), ..., x(N)}排列熵的计算可以分解为以下几步2.1 相空间重构从一维序列到多维视图这是几乎所有非线性时间序列分析的第一步目的是还原系统潜在的动力学结构。对于排列熵我们采用最常用的时间延迟嵌入法。嵌入维度 (m) 这是最重要的参数决定了我们每次观察多少个连续的数据点。你可以把它理解为“观察窗口”的长度。例如m3就意味着我们每次截取3个连续的点作为一个观察单元。时间延迟 (τ) 通常取1即连续采样。对于某些特定周期的信号可以取其他值以捕获不同时间尺度上的动力学但在绝大多数通用场景下τ1是默认且安全的选择。重构后我们得到一系列m维的向量X(i) [x(i), x(iτ), ..., x(i(m-1)τ)], 其中i 1, 2, ..., N-(m-1)τ。 每个X(i)都是原始信号的一个“快照”。2.2 序数模式编码将向量转化为“排名符号”这是排列熵的灵魂所在。对于每一个m维向量X(i)我们不再关心它的具体数值而是关心其内部元素的大小排列顺序。例如假设m3我们有一个向量[2.1, 5.6, 1.8]。我们先将这三个数排序1.8 2.1 5.6。排序后它们的原始索引是[3, 1, 2]因为最小值1.8在原数组第3位次小值2.1在第1位最大值5.6在第2位。这个索引序列(3,1,2)就定义了一种序数模式Ordinal Pattern。对于m3所有可能的序数模式是3! 6种即(1,2,3), (1,3,2), (2,1,3), (2,3,1), (3,1,2), (3,2,1)。每一种模式都代表了数据点之间一种特定的局部形态关系比如(1,2,3)代表单调递增(3,2,1)代表单调递减(2,1,3)代表先降后升的“V”形。注意如果向量中存在相等的值平局标准的排列熵算法通常采用其出现顺序的索引但更严谨的做法是给它们分配相同的排名例如使用tiedrank函数。在实际生理或振动信号中由于是连续值完全相等的概率极低但量化后的数据可能出现。我提供的源码中包含了处理这种情况的稳健逻辑。2.3 统计与熵值计算衡量模式的“混乱度”经过上一步我们把长度为N的时间序列转化为了一个由K N - (m-1)个序数模式符号组成的序列每个符号是1到m的一个排列。接下来我们统计这m!种可能的序数模式在整段序列中出现的频率。假设第j种模式出现了C_j次那么其相对频率为p_j C_j / K。最后套用香农熵的公式来计算排列熵H_p(m)H_p(m) - Σ (p_j * log2(p_j))求和范围是所有p_j 0的模式。熵值H_p(m)的范围是[0, log2(m!)]。为了便于在不同m值下比较通常会进行归一化处理H_p_normalized H_p(m) / log2(m!)此时值域为[0, 1]。为什么这样做是聪明的抗噪性强小幅度的噪声可能会改变信号的幅值但不太可能颠覆多个连续点之间的大小顺序关系。因此排列熵对加性噪声不敏感。计算快速核心操作是排序和计数算法复杂度约为O(N*m*log(m))比需要计算距离矩阵的近似熵快得多。概念直观它直接度量了序列局部形态的多样性和不可预测性物理意义明确。3. 关键参数选择让排列熵在你的数据上发挥最大效力排列熵用起来简单但想把它的威力完全发挥出来参数选择是关键。这里没有放之四海而皆准的“黄金值”必须结合你的数据特性来定。3.1 嵌入维度m决定观察的细致程度m的选择是核心中的核心。它直接决定了你能分辨多少种不同的局部模式。m太小如2或3模式种类太少2!2,3!6分辨率不足可能无法区分一些动力学行为相似的复杂系统。例如一个4周期振荡和一个8周期振荡在m2下可能表现出相似的熵值。m太大如7或8模式种类急剧增加7!5040,8!40320。对于有限长度的数据N会导致绝大多数模式只出现一次或根本不出现使得计算出的概率分布不可靠熵值估计产生较大误差。同时计算量也会增加。经验法则Bandt 和 Pompe 在提出该方法的原始论文中建议m在3到7之间。我的实践经验是对于大多数工程和生物信号数据长度N在1000到10000量级m4或m5是一个非常好的起点能在区分度和统计可靠性之间取得很好的平衡。你可以做一个简单的验证逐步增加m观察归一化排列熵值的变化。当m增加到某个值后熵值趋于稳定或开始不规则波动那么前一个m值通常是合适的。3.2 时间延迟τ通常就选1除非你有很强的先验知识知道你的信号在某个特定时间尺度上有重要特征例如心电信号的R-R间期否则τ1即使用连续采样点是最通用、最不容易出错的选择。选择τ1相当于对信号进行下采样可能会丢失高频细节信息。3.3 数据长度N越长越稳但有下限为了保证统计的可靠性序列长度N需要远大于可能出现的模式种类m!。一个常用的经验法则是N 5 * m!。例如当m5时m! 120那么N最好大于600。对于m6(720)则要求N 3600。这也是为什么通常不推荐使用m7的原因——对数据量的要求太高。我提供的源码中包含了数据长度检查的警告提示。3.4 归一化为了跨尺度比较如前所述使用归一化的排列熵(H_p / log2(m!))几乎总是更好的选择。这使得不同m值下计算出的熵值可以放在同一个尺度[0,1]内进行比较结果更直观。4. 实战进阶多尺度排列熵与MATLAB源码深度解析基础的排列熵只能反映单一时间尺度上的复杂度。然而很多真实世界的信号如机械振动、脑电波其复杂性体现在多个时间尺度上。为此Costa等人提出了多尺度排列熵。4.1 多尺度排列熵的计算流程MSPE的核心思想是先对原始序列进行粗粒化得到不同尺度因子下的新序列再对每个新序列计算排列熵。粗粒化过程对于尺度因子s将原始序列{x1, x2, ..., xN}划分为长度为s的不重叠窗口计算每个窗口内数据的均值形成粗粒化序列{y_j^(s)}。y_j^(s) (1/s) * Σ_{i(j-1)s1}^{js} x_i,j 1, 2, ..., floor(N/s)当s1时就是原始序列。计算排列熵对每一个尺度因子s下的粗粒化序列y^(s)计算其排列熵PE(s)。绘制MSPE曲线以尺度因子s为横坐标PE(s)为纵坐标即可得到一条多尺度复杂度曲线。健康的、适应性强的系统如正常人的心率通常在不同尺度上表现出较高的熵值而病态的、僵化的系统如心衰患者的心率其熵值可能在多个尺度上降低。4.2 MATLAB源码实现与关键技巧以下是我提供的PermutationEntropy.m函数的核心部分解析并附上关键技巧说明。function [pe, patterns] PermutationEntropy(data, m, tau, normalize) % 计算时间序列的排列熵 % 输入 % data - 一维时间序列 (行向量或列向量) % m - 嵌入维度 (默认 4) % tau - 时间延迟 (默认 1) % normalize - 是否归一化到[0,1] (默认 true) % 输出 % pe - 排列熵值 % patterns - (可选) 所有序数模式的分布直方图数据 % % 作者基于Bandt Pompe (2002) 算法结合工程实践优化 if nargin 4 || isempty(normalize), normalize true; end if nargin 3 || isempty(tau), tau 1; end if nargin 2 || isempty(m), m 4; end data data(:); % 确保是列向量 N length(data); % 参数有效性检查 if m 2 error(嵌入维度 m 必须大于等于2.); end if tau 1 error(时间延迟 tau 必须为正整数.); end if N m * tau error(数据长度 N 必须大于 m * tau.); end % 数据长度建议检查 if N 5 * factorial(m) warning(数据长度可能不足以保证统计可靠性。建议 N 5*m!。当前 N%d, m!%d, N, factorial(m)); end % 1. 相空间重构 numVec N - (m-1)*tau; % 重构向量的数量 emb zeros(numVec, m); for i 1:m emb(:, i) data((1:numVec) (i-1)*tau); end % 2. 生成序数模式并统计 % 使用 sort 函数获取排序索引效率高且能处理平局情况 [~, idx] sort(emb, 2); % 沿行排序idx 的每一行就是一个序数模式 % 将索引模式转换为唯一的整数标签便于快速统计 % 技巧利用阶乘进制数将排列映射为唯一整数 factVec factorial(m-1:-1:0); patternLabels sum((idx - 1) .* factVec, 2) 1; % 统计每种标签出现的次数 possiblePatterns factorial(m); count histcounts(patternLabels, 1:possiblePatterns1); % 3. 计算概率和熵值 prob count(count 0) / numVec; % 只考虑出现的模式 pe -sum(prob .* log2(prob)); % 4. 归一化 if normalize pe pe / log2(possiblePatterns); end % 5. (可选) 输出模式分布 if nargout 1 patterns count; end end关键技巧与避坑点高效的模式编码代码中没有使用循环去比较每个模式而是利用sort函数一次性得到所有向量的排序索引然后通过阶乘进制数的技巧将每个索引序列如[2,3,1]映射为一个唯一的整数标签。这比用字符串比较或循环查找快了几个数量级尤其是当N很大时。平局处理MATLAB的sort函数在默认情况下是稳定排序即对于相等的值会保持它们在原始数据中的相对顺序。这提供了一种处理平局理论上极少但量化数据可能出现的合理方式。如果你希望更严格地处理平局如视为同排名可以使用tiedrank函数但会显著增加计算量。数据长度警告函数内置了长度检查警告。这是非常重要的调试辅助。如果你发现熵值不稳定或者不同数据段结果差异巨大首先应该怀疑数据长度是否足够。向量化操作整个重构和排序过程完全向量化没有显式循环充分利用了MATLAB的矩阵运算优势速度极快。4.3 多尺度排列熵的封装实现基于上面的核心函数实现多尺度版本就很简单了。以下是MultiscalePermutationEntropy.m的框架function [mspe, scales] MultiscalePermutationEntropy(data, m, tau, maxScale, normalize) % 计算多尺度排列熵 % 输入 % data - 一维时间序列 % m, tau, normalize - 同 PermutationEntropy 函数 % maxScale - 最大尺度因子 (默认 20) % 输出 % mspe - 各尺度下的排列熵值 (向量) % scales - 尺度因子 (向量) if nargin 5, normalize true; end if nargin 4, maxScale 20; end if nargin 3, tau 1; end if nargin 2, m 4; end data data(:); N length(data); scales 1:maxScale; mspe zeros(size(scales)); for i 1:length(scales) s scales(i); % 粗粒化 if s 1 coarseData data; else numCoarse floor(N / s); % 使用 reshape 和 mean 进行快速粗粒化 (需确保能整除) if mod(N, s) 0 coarseData mean(reshape(data(1:numCoarse*s), s, numCoarse), 1); else % 如果不能整除采用循环较慢或截断尾部 % 这里采用截断保证计算速度并给出提示 coarseData mean(reshape(data(1:numCoarse*s), s, numCoarse), 1); if i 1 % 只提示一次 warning(数据长度不能被所有尺度整除尾部数据将被截断用于粗粒化。); end end end % 计算当前尺度下的排列熵 mspe(i) PermutationEntropy(coarseData, m, tau, normalize); end end避坑点粗粒化时的长度处理粗粒化要求数据长度N最好是尺度因子s的整数倍。如果不能整除常见的处理方法是1) 截断尾部多余数据2) 使用滑动平均重叠窗口。我的代码采用了第一种方法截断因为它计算更快且对于尺度分析轻微的尾部截断通常影响不大。但在严谨的研究中你需要明确说明处理方法。第二种方法滑动窗口会引入数据间的相关性可能影响熵值的解释。5. 应用实例轴承故障振动信号分析让我们用一个实际的例子看看如何用这套代码分析轴承振动信号区分正常状态和早期故障。5.1 数据准备与预处理假设我们有两段振动加速度信号采样频率fs 12000 Hz时长2秒即每段N 24000点。一段来自正常轴承 (data_normal)一段来自外圈有轻微故障的轴承 (data_fault)。数据通常会有直流偏移或趋势项需要先去除。% 1. 数据预处理去趋势可选带通滤波 data_normal detrend(data_normal); % 去除线性趋势 data_fault detrend(data_fault); % 可选根据轴承特征频率设计一个带通滤波器保留感兴趣频段 % [b, a] butter(4, [100, 2000]/(fs/2), bandpass); % data_normal filtfilt(b, a, data_normal); % data_fault filtfilt(b, a, data_fault);5.2 计算单一尺度排列熵首先我们用默认参数快速计算一下两个信号的排列熵看看能否直接区分。m 4; tau 1; pe_normal PermutationEntropy(data_normal, m, tau, true); pe_fault PermutationEntropy(data_fault, m, tau, true); fprintf(正常信号排列熵 (m%d): %.4f\n, m, pe_normal); fprintf(故障信号排列熵 (m%d): %.4f\n, m, pe_fault);结果解读通常早期故障会引入周期性冲击使信号的规律性增强局部模式变得不那么“混乱”因此排列熵值会降低。所以我们很可能看到pe_fault pe_normal。但这并非绝对取决于故障类型和位置有时熵值也可能升高。单一尺度的结果只是一个初步的、粗略的指标。5.3 计算多尺度排列熵并可视化单一尺度可能不够稳健我们计算多尺度排列熵来获取更全面的特征。maxScale 30; [mspe_normal, scales] MultiscalePermutationEntropy(data_normal, m, tau, maxScale, true); [mspe_fault, ~] MultiscalePermutationEntropy(data_fault, m, tau, maxScale, true); % 可视化 figure(Position, [100, 100, 800, 400]); subplot(1,2,1); plot(scales, mspe_normal, b-o, LineWidth, 1.5, MarkerSize, 6, DisplayName, 正常); hold on; plot(scales, mspe_fault, r-s, LineWidth, 1.5, MarkerSize, 6, DisplayName, 故障); xlabel(尺度因子 s); ylabel(归一化排列熵 PE(s)); title(多尺度排列熵曲线); legend(Location, best); grid on; subplot(1,2,2); % 计算熵值曲线下的面积 (AUC) 作为一个综合特征 auc_normal trapz(scales, mspe_normal); auc_fault trapz(scales, mspe_fault); bar([1,2], [auc_normal, auc_fault]); set(gca, XTickLabel, {正常, 故障}); ylabel(MSPE曲线下面积 (AUC)); title(综合特征对比); grid on;分析与解读曲线形态正常轴承的MSPE曲线可能在多个尺度上都维持在相对较高的水平表明其动力学行为在不同时间尺度上都具有一定的复杂性。而故障轴承的曲线可能在低尺度反映细节或高尺度反映整体趋势上出现明显的下降表明故障冲击破坏了系统原有的复杂动力学结构使其行为变得更规则、更可预测。综合特征直接比较整条曲线有时不够直观。我们可以提取曲线的统计特征作为故障诊断的输入特征例如曲线下面积如上图所示AUC是一个很好的综合指标故障状态下AUC通常会减小。均值、方差计算所有尺度上熵值的均值和标准差。斜率计算曲线在低尺度段和高尺度段的平均斜率故障可能改变尺度间的复杂度传递特性。参数敏感性测试在实际应用中我强烈建议你对关键参数m进行敏感性分析。固定其他参数让m从3变化到6分别计算MSPE。观察不同m下正常与故障状态的区分度是否稳定。如果某个m值下区分度最好且稳定那么就选用该值。5.4 模式分布的可视化洞察动力学变化除了熵值这个汇总统计量观察序数模式本身的分布变化也能提供深刻见解。% 计算并可视化模式分布 [~, patterns_normal] PermutationEntropy(data_normal, 4, 1, false); [~, patterns_fault] PermutationEntropy(data_fault, 4, 1, false); figure; subplot(2,1,1); bar(patterns_normal / sum(patterns_normal)); xlabel(序数模式索引 (m4, 共24种)); ylabel(相对频率); title(正常轴承信号序数模式分布); xlim([0.5, 24.5]); grid on; subplot(2,1,2); bar(patterns_fault / sum(patterns_fault)); xlabel(序数模式索引 (m4, 共24种)); ylabel(相对频率); title(故障轴承信号序数模式分布); xlim([0.5, 24.5]); grid on;解读在m4时有24种可能的模式。健康系统的模式分布通常相对“平坦”即没有哪种模式占据绝对主导。而故障系统由于周期性冲击的引入可能会使某几种特定的模式例如对应冲击上升沿或下降沿的形态出现频率显著增高导致分布出现尖峰。这种分布的变化是熵值降低的根本原因也为理解故障的物理机制提供了线索。6. 常见问题、局限性与扩展方向即使有了好用的工具明白它的边界同样重要。6.1 常见踩坑点与解决方案熵值结果不稳定每次跑略有差异可能原因数据长度N严重不足。请检查N是否满足N 5*m!。对于短数据排列熵估计的方差会很大。解决方案增加数据长度或使用更小的m如从5降到4。也可以考虑对长数据分多个不重叠的段计算熵值然后取平均和标准差作为最终结果和置信区间。故障信号的熵值反而比正常信号高可能原因这并非不可能。如果故障引入了强烈的随机噪声而非周期性冲击或者使系统进入混沌状态熵值可能会升高。此外参数选择不当如m太大导致正常信号熵值被低估也可能造成反直觉结果。解决方案首先检查参数尤其是m是否合适。其次结合其他诊断方法如频谱分析、包络谱综合判断。排列熵只是一个特征不是万能判据。计算速度慢处理长数据时卡顿可能原因使用了低效的循环实现模式匹配。解决方案使用我提供的基于向量化和整数编码的算法速度会有数量级提升。对于超长序列如 10^6 点可以考虑先对数据分段计算或进行适当的降采样。对脉冲类故障敏感但对磨损类缓慢退化不敏感问题本质排列熵对序列的局部形态变化敏感。明显的冲击脉冲会改变局部形态分布。而均匀磨损可能在整个时间尺度上缓慢地改变信号能量但对局部点之间的序关系影响较小。解决方案将排列熵与反映能量变化的特征如RMS、峰峰值、小波能量熵结合使用构建多维特征向量。6.2 方法局限性丢失幅值信息这是其优点抗噪也是缺点。两个幅值分布完全不同但局部排序模式相似的序列会有相同的排列熵。因此它需要与其他时域、频域特征互补。对慢变趋势敏感一个强烈的线性趋势会主导局部排序导致模式分布严重偏向某几种如单调递增从而低估序列的真实动力学复杂度。因此在计算前对数据去趋势detrend是至关重要的预处理步骤。尺度粗粒化的信息损失MSPE中的粗粒化取平均是一个低通滤波过程会损失高频细节。有学者提出了复合粗粒化、层次熵等改进方案来弥补这一点。6.3 扩展与变体加权排列熵在计算模式概率时不仅考虑模式是否出现还考虑产生该模式的向量中数据点之间的实际数值差异。给差异大的模式赋予更大的权重从而在保留排序信息的同时引入部分幅值信息。模糊排列熵引入模糊隶属度的概念缓和了传统排列熵中“非此即彼”的硬判决使模式划分更平滑对噪声和参数选择的鲁棒性更强。多变量排列熵将算法扩展到多个相关时间序列用于分析系统不同变量间的耦合复杂度。在我提供的完整源码包MATLAB实现排列熵 程序源码.zip中除了上述核心函数还包含了用于批量处理数据、生成特征矩阵、以及绘制出版级质量图形的辅助脚本和示例数据。真正上手时建议你从示例开始替换成自己的数据调整参数观察变化这样才能最深切地体会这个工具的妙处和边界。信号分析没有银弹但排列熵无疑是你工具箱里一把锋利而独特的瑞士军刀。本文还有配套的精品资源点击获取