简介这份代码包提供 Lomb-Scargle 周期图在 Matlab 中的完整实现面向需要分析非均匀时间序列的天文、地球科学及工程领域用户可用于探测周期性信号并评估其统计显著性。包内共 2 个文件均为 m 脚本lombscargle.m 实现核心算法覆盖数据预处理、线性谐波模型构建、最小二乘拟合、功率谱计算与结果返回支持自定义频率范围可处理不等间隔观测数据inputtolomb.m 则为输入样例演示如何准备时间戳、观测值与感兴趣频率区间帮助使用者快速理解调用方式。压缩包仅 7KB非常轻量。目前已有 493 人学习下载。读者可据此掌握 LSP 分析的完整流程包括频率扫描、周期识别与显著性判断代码结构清晰、便于替换和扩展可直接迁移到天体测光、地球物理监测或工程振动分析等场景中在多周期识别与显著性阈值对比时获得可靠参考结果。 第一次看到lombscargle.m这个文件名是在处理一批不规则的振动监测数据时。传感器每隔几秒采样一次中途还因为网络抖动断了好几次数据摆在那儿明明感觉里面有个周期可当我把时间轴对齐到均匀网格再跑fft结果完全是另一回事。后来接触了Lomb-Scargle周期图才明白对这类非均匀采样时间序列直接做傅里叶变换是行不通的。这个MATLAB脚本做的事情就是绕过“均匀采样”这个前提用最小二乘正弦拟合的方式直接估算频率谱。今天我从原理、函数使用到自编实现把这套东西完整梳理一遍给同样在处理非均匀时序的朋友做个参考。1. lombscargle.m到底解决什么问题1.1 非均匀采样在真实场景中其实很常见很多人在学校用MATLAB做信号处理时接触到的数据都是规规矩矩的等间隔采样t 0:0.001:1然后x sin(2*pi*50*t)最后fft一把梭。这种操作在实验室里确实没问题但到了现场工程中数据往往没有这么听话。就拿我接触过的旋转机械振动监测来说传感器读数虽然名义上是等间隔的但采集系统负载一高报文就会延迟时间戳记录下来的实际采样间隔可能在100ms到300ms之间随机抖动。更极端的情况是天文观测望远镜只有在夜间才能观测遇到阴天还要停工得到的时间序列天然就是“这里一段、那里一段”的非均匀采样。这种数据如果直接做FFT第一步就卡住了——FFT要求输入数据在时间轴上是等距的。为了能用FFT很多人会选择先插值再重采样把数据强行映射到均匀时间网格上。这个操作本身就会引入误差插值相当于在数据里面添加了人为的低通滤波会抑制高频成分更麻烦的是它可能把真实信号的能量“抹”到相邻频点上让本来干净的频谱多出一堆假峰。1.2 Lomb-Scargle周期图的定位不等间隔数据的频谱分析工具Lomb-Scargle周期图最早由Lomb在1976年提出后来Scargle在1982年做了系统化整理主要就是为了解决天文学中不等间隔观测数据的周期搜寻问题。它的核心思路很直接对每一个候选频率用最小二乘法去拟合一条正弦曲线看这个频率能在多大程度上解释数据中的波动。所有频率下的拟合优度连起来就是周期图。这个做法不需要数据等间隔也不需要对缺失部分做任何插值补全。相比于插值后FFTLomb-Scargle有一个非常实在的优点它把采样时刻本身当成已知信息参与运算所以不会“凭空捏造”数据点的位置。也就是说只要原始时间戳是准的周期图的结果就能真实反映数据中的周期性成分。这也是为什么后来它被大量用在生物节律研究、气候分析、心电图异常检测以及我所在领域的设备故障诊断中——只要你手上的时间序列不是规规矩矩等间隔的lombscargle这一类的工具就大概率比FFT靠谱。2. 原理拆解为什么“拟合正弦波”能替代傅里叶变换2.1 从最小二乘拟合到Lomb-Scargle公式要理解Lomb-Scargle的工作原理可以把它看成是一系列单频正弦拟合的循环。假设我们有一组观测值(x_i)对应的时间戳是(t_i)其中(i 1,2,\dots,N)。对于给定的频率(f)我们希望找到一组系数(a,b)使得模型[ x_i a \cos(2\pi f t_i) b \sin(2\pi f t_i) ]能最大程度地接近实际观测值。这里的“接近”用最小二乘来度量也就是让残差平方和最小。如果某个频率对应的拟合效果特别好说明数据里确实存在这个频率的周期分量如果拟合效果很差说明该频率在数据中不明显。Scargle后来把这个问题整理成了一个便于计算的统计量。先对数据做中心化处理令(X_i x_i - \bar{x})其中(\bar{x})是均值。然后对于每个频率(f)定义一个与时间偏移有关的相位量(\tau)[ \tan(2\omega \tau) \frac{\sum_i \sin(2\omega t_i)}{\sum_i \cos(2\omega t_i)} ]其中(\omega 2\pi f)。这个偏移的目的是让正弦和余弦基函数在采样时刻上保持正交性从而让功率计算变得稳定。最终的周期图功率定义为[ P(f) \frac{1}{2\sigma^2}\left{ \frac{\left[\sum_i X_i \cos\omega(t_i-\tau)\right]^2}{\sum_i \cos^2\omega(t_i-\tau)} \frac{\left[\sum_i X_i \sin\omega(t_i-\tau)\right]^2}{\sum_i \sin^2\omega(t_i-\tau)} \right} ]其中(\sigma^2)是数据的方差。这个式子的计算复杂度看起来是(O(N))每个频率如果我们扫几千个频率点计算量就是(O(N \times M))。对现代计算机来说几万条数据配几千个频率点完全不是问题。2.2 和FFT的关系到底谁更“准”有人可能会问既然傅里叶变换在等间隔采样下运行极快(O(N\log N))为什么非均匀采样不直接套FFT偏要用Lomb-Scargle这种逐频率拟合的笨办法关键在于FFT的数学推导依赖一个前提采样点均匀分布在时间轴上。一旦这个前提被打破FFT的基函数之间不再正交计算出来的频谱中会混入各频率之间的“串扰”。用插值硬凑均匀只是牺牲分辨率换来了一个“能跑的FFT”而不是真正准确的频谱。Lomb-Scargle的拟合过程虽然慢一些但它在每个频率上做的是真正的投影不会因为其他频率成分的存在而干扰当前频率的功率估计。尤其是当数据存在大量缺口、采样间隔差异明显时Lomb-Scargle往往能恢复出真实的谱峰而插值FFT会在谱峰旁边拖出一长片泄漏。当然如果数据本身已经是均匀采样的FFT无论在速度还是精度上都是更好的选择Lomb-Scargle只是为“非均匀”这个场景兜底的方案。3. MATLAB实操用plomb函数实现周期估计3.1 数据准备先模拟一组非均匀采样信号在MATLAB里最简单的落地方案是直接用信号处理工具箱自带的plomb函数。它的名字虽然不叫lombscargle但底层就是Lomb-Scargle周期图的标准实现。先用一个仿真例子跑通流程方便说明参数设置。rng(42); % 构造均匀时间轴目标频率5 Hz t_uniform 0:0.01:60; f_true 5; x_clean sin(2*pi*f_true*t_uniform) 0.3*randn(size(t_uniform)); % 随机抽取20%的点得到非均匀采样序列 keep rand(size(t_uniform)) 0.2; t sort(t_uniform(keep)); x x_clean(keep);这里故意只保留20%的点相当于采样率很低而且时间间隔完全不均匀。如果用插值后FFT插值点的密度不够结果会很糟糕但Lomb-Scargle可以直接利用这些不规则的采样时刻。画个散点图先看一眼数据形态figure; stem(t, x, filled, MarkerSize, 3); xlabel(Time (s)); ylabel(Amplitude); title(Non-uniformly sampled data);你会看到时间轴上点分布集中在某些区域、又缺掉某些区域这就是典型的非均匀采样。3.2 调用plomb并解读输出调用方式非常简单[p, f] plomb(x, t); plot(f, p); xlabel(Frequency (Hz)); ylabel(Power); xlim([0, 10]);plomb默认会做去均值处理所以即使信号有直流分量也不会在0频处出现大尖峰。输出p是功率谱密度f是对应的频率轴单位取决于输入时间戳的单位。运行后你应该能在5 Hz附近看到一个明显的尖峰这正是我们埋进去的真实频率。这里有个重要经验如果数据里混有线性趋势光靠plomb默认的去均值是不够的。比如传感器的缓慢漂移会在低频段形成一个“陡坡”把真正的周期峰淹没。所以实际操作中我习惯于先对数据做一次线性去趋势x_detrend detrend(x, linear); [p, f] plomb(x_detrend, t);去趋势之后低频段的能量会明显减小中高频的周期峰才能浮现出来。这个细节不处理很多时候你会看到周期图整体呈“右下倾斜”的形状第一个尖峰永远出现在最低频处其实那个不代表真实周期只是趋势的贡献。3.3 频率范围与结果解读plomb默认会基于数据的时间跨度自动选择频率范围。不过在某些场景下自动范围可能包含过高的频率导致计算量偏大且产生大量伪峰。我们可以手动指定频率范围freq linspace(0, 10, 2000); [p, f] plomb(x_detrend, t, freq);这里的freq实际上是一个频率网格plomb会在这些频率点上计算功率。网格越密频率分辨率越细但计算量线性增加。我的经验是先根据物理背景确定“可能存在的周期范围”再在这个范围内均匀取1000~5000个频率点就足够用了。拿到周期图后找峰通常用findpeaks[pks, locs] findpeaks(p, f, SortStr, descend, NPeaks, 3);但要注意一点findpeaks找出来的峰值不一定都显著。这时需要评估“误报概率”False Alarm Probability, FAP。plomb内部提供了相关计算具体参数可以查MATLAB文档中关于FAP/Pd的选项。核心思想是如果某峰值的高度远高于噪声基底那它偶然出现的概率就极低对应FAP很小说明周期可信。不要只看峰的高度差要结合统计显著性判断尤其是低信噪比数据。4. 自己编写lombscargle.m时的关键细节4.1 去均值、去趋势与归一化虽然直接用plomb很方便但有时候我们得在旧版MATLAB或者没有信号处理工具箱的环境下工作这时自己写一个lombscargle.m就很有必要。我在自编时踩过几个坑第一个就是数据中心化。如果数据有一个较大的直流偏置而公式中直接用原始(x_i)去计算正弦拟合的分子计算结果会在低频段出现巨大峰值因为直流分量对任意低频正弦都有很强的投影。正确的做法是先把均值减掉变成(X_i x_i - \bar{x})。更进一步如果数据存在缓慢漂移最好先做一次线性去趋势否则周期图低频部分会被趋势主导。我习惯把函数入口设计成function [p, f] lombscargle(t, x, fgrid) % LOMBSCARGLE 计算非均匀采样时间序列的Lomb-Scargle周期图 % 输入 % t - 时间戳列向量 % x - 观测值列向量长度与t相同 % fgrid - 频率网格用于计算的频率点 % 输出 % p - 周期图功率 % f - 与p对应的频率4.2 频率网格过采样因子和最高频率的选择自编时最影响结果质量的是频率网格的选择。很多文献里提到两个参数ofacoversampling factor和hifachighest frequency factor。ofac用来控制频率网格的过采样程度通常取值4~20hifac用来控制最大扫描频率通常取平均采样间隔倒数的一半乘以某个系数。没有现成工具箱的话可以按下面逻辑设定频率上限% 平均采样间隔 dt_mean mean(diff(t)); % 一个保守的频率上限 fmax 1 / (2 * dt_mean); % 频率网格步长取单个频点对应的扫描宽度 df 1 / (range(t) * 4); fgrid (0:df:fmax).;这里range(t)是数据的总时间跨度df相当于在傅里叶分析中“频率分辨率”的四分之一目的是让谱峰更平滑便于找峰。如果数据里有很短的采样间隔用1/(2*dt_min)算出来的上限会非常高导致计算量爆炸。实际操作时我通常会结合信号物理特征把fmax限制在一个合理范围内比如机械故障诊断只看0~500 Hz天文周期搜索只扫0~1/天没必要扫描所有理论可用的频率。4.3 自编实现的计算效率窍门最朴素的实现是双层循环外层扫频率内层对每个频率做正弦拟合。这在数据量只有几百、频率点只有几百时没问题但数据量一大就非常慢。推荐的做法是向量化把时间向量写成行向量频率网格写成列向量然后用矩阵运算一次性算出所有频率下的正弦和余弦值。以核心的分子计算为例可以这样向量化omega 2 * pi * fgrid; % 频率网格角频率列向量 T t(:).; % 时间行向量 cos_wt cos(omega .* T); % 矩阵每个频率×每个时间点 sin_wt sin(omega .* T); % 重心化后的数据 xw x(:) - mean(x); % 计算每个频率下的余弦基投影 cos_sum cos_wt * xw; % 相当于 sum(xw * cos(omega_i * t_k))不过要注意Lomb-Scargle公式中还有一个相位偏移(\tau)也就是说直接用cos(omega .* T)和sin(omega .* T)做投影在非均匀采样时存在基函数不正交的问题。简单实现可以忽略(\tau)直接算得到的结果和标准版有细微差别严格实现需要先对每个频率求出(\tau)再做一次旋转。我建议初学者先用plomb对比验证自编结果如果差异在可接受范围内可以接受简化版毕竟很多实际判断只看峰位不看绝对功率。5. 常见问题与避坑经验5.1 伪峰从哪里来非均匀采样下最大的坑是“混叠”。由于采样间隔不均匀某些频率的真实信号可能会在另一个频率处产生虚假的响应。这有点像街头摄影中飞驰而过的车轮拍到的轮辐看起来转得很慢甚至反向转因为你采样时刻刚好“欺骗”了视觉。要降低伪峰风险有几个实操手段。一是尽量保证数据中缺失的块不要太大不要出现长时间完全无观测的“空洞”二是多个频率网格参数交叉验证在ofac4和ofac16下都算一遍真实峰的位置应该稳定不动伪峰则经常随网格变化三是利用FAP筛选峰值。三者结合基本能避免把随机波动当成周期。5.2 数据长度与信噪比阈值Lomb-Scargle统计量在数据量较小时有较高的不确定性。如果总观测时长只覆盖了目标周期的两三个周期即使数据确实存在周期性谱峰也会很宽、很难精确定位。我的经验是至少要让时间跨度覆盖目标周期的5倍以上峰位才能可靠。假如你的数据只有60秒就别指望能从里面可靠提取周期大于30秒的信号。另外噪声不能太大。一个粗略的参考是如果数据信噪比低于3 dB周期图上的峰就可能淹没在噪声中肉眼几乎看不出。这时候不要盲目相信plomb的输出而是尝试对数据做带通滤波预处理或者分段后重新计算。5.3 工具箱依赖与版本问题plomb函数需要Signal Processing Toolbox如果用的是轻量级MATLAB授权或者手头只有MATLAB Online基础版有可能调用时报错“未定义函数或变量”。这时候自写的lombscargle.m就派上用场了。还有一个坑是频率单位如果时间戳t的单位是小时那么plomb返回的频率单位是“周/小时”不是默认的Hz。我一度因为时间单位用错把算出来的周期读成了完全不对的秒数。建议在代码开头就把时间统一换算成秒后面所有频率解读都基于单位Hz。最后再分享一个我自己调试的小习惯拿到一批非均匀采样数据后我会先用一个已知周期的人工合成信号做一遍全流程测试确认代码、参数、解读方式没问题再处理真实数据。这样做的好处是你对结果长什么样心里有数真数据跑出奇怪结果时能更快判断是算法问题、参数问题还是数据本身的问题。非均匀采样的周期分析并不神秘关键的思路就是尊重每个采样点的真实时刻用拟合代替插值这也是lombscargle.m这类代码能站住脚的根本原因。本文还有配套的精品资源点击获取 
长河君 · 2026/9/3 3:02:16 · 阅读 1.2k MATLAB实现Lomb-Scargle周期图：非均匀采样时间序列频谱分析
简介这份代码包提供 Lomb-Scargle 周期图在 Matlab 中的完整实现面向需要分析非均匀时间序列的天文、地球科学及工程领域用户可用于探测周期性信号并评估其统计显著性。包内共 2 个文件均为 m 脚本lombscargle.m 实现核心算法覆盖数据预处理、线性谐波模型构建、最小二乘拟合、功率谱计算与结果返回支持自定义频率范围可处理不等间隔观测数据inputtolomb.m 则为输入样例演示如何准备时间戳、观测值与感兴趣频率区间帮助使用者快速理解调用方式。压缩包仅 7KB非常轻量。目前已有 493 人学习下载。读者可据此掌握 LSP 分析的完整流程包括频率扫描、周期识别与显著性判断代码结构清晰、便于替换和扩展可直接迁移到天体测光、地球物理监测或工程振动分析等场景中在多周期识别与显著性阈值对比时获得可靠参考结果。 第一次看到lombscargle.m这个文件名是在处理一批不规则的振动监测数据时。传感器每隔几秒采样一次中途还因为网络抖动断了好几次数据摆在那儿明明感觉里面有个周期可当我把时间轴对齐到均匀网格再跑fft结果完全是另一回事。后来接触了Lomb-Scargle周期图才明白对这类非均匀采样时间序列直接做傅里叶变换是行不通的。这个MATLAB脚本做的事情就是绕过“均匀采样”这个前提用最小二乘正弦拟合的方式直接估算频率谱。今天我从原理、函数使用到自编实现把这套东西完整梳理一遍给同样在处理非均匀时序的朋友做个参考。1. lombscargle.m到底解决什么问题1.1 非均匀采样在真实场景中其实很常见很多人在学校用MATLAB做信号处理时接触到的数据都是规规矩矩的等间隔采样t 0:0.001:1然后x sin(2*pi*50*t)最后fft一把梭。这种操作在实验室里确实没问题但到了现场工程中数据往往没有这么听话。就拿我接触过的旋转机械振动监测来说传感器读数虽然名义上是等间隔的但采集系统负载一高报文就会延迟时间戳记录下来的实际采样间隔可能在100ms到300ms之间随机抖动。更极端的情况是天文观测望远镜只有在夜间才能观测遇到阴天还要停工得到的时间序列天然就是“这里一段、那里一段”的非均匀采样。这种数据如果直接做FFT第一步就卡住了——FFT要求输入数据在时间轴上是等距的。为了能用FFT很多人会选择先插值再重采样把数据强行映射到均匀时间网格上。这个操作本身就会引入误差插值相当于在数据里面添加了人为的低通滤波会抑制高频成分更麻烦的是它可能把真实信号的能量“抹”到相邻频点上让本来干净的频谱多出一堆假峰。1.2 Lomb-Scargle周期图的定位不等间隔数据的频谱分析工具Lomb-Scargle周期图最早由Lomb在1976年提出后来Scargle在1982年做了系统化整理主要就是为了解决天文学中不等间隔观测数据的周期搜寻问题。它的核心思路很直接对每一个候选频率用最小二乘法去拟合一条正弦曲线看这个频率能在多大程度上解释数据中的波动。所有频率下的拟合优度连起来就是周期图。这个做法不需要数据等间隔也不需要对缺失部分做任何插值补全。相比于插值后FFTLomb-Scargle有一个非常实在的优点它把采样时刻本身当成已知信息参与运算所以不会“凭空捏造”数据点的位置。也就是说只要原始时间戳是准的周期图的结果就能真实反映数据中的周期性成分。这也是为什么后来它被大量用在生物节律研究、气候分析、心电图异常检测以及我所在领域的设备故障诊断中——只要你手上的时间序列不是规规矩矩等间隔的lombscargle这一类的工具就大概率比FFT靠谱。2. 原理拆解为什么“拟合正弦波”能替代傅里叶变换2.1 从最小二乘拟合到Lomb-Scargle公式要理解Lomb-Scargle的工作原理可以把它看成是一系列单频正弦拟合的循环。假设我们有一组观测值(x_i)对应的时间戳是(t_i)其中(i 1,2,\dots,N)。对于给定的频率(f)我们希望找到一组系数(a,b)使得模型[ x_i a \cos(2\pi f t_i) b \sin(2\pi f t_i) ]能最大程度地接近实际观测值。这里的“接近”用最小二乘来度量也就是让残差平方和最小。如果某个频率对应的拟合效果特别好说明数据里确实存在这个频率的周期分量如果拟合效果很差说明该频率在数据中不明显。Scargle后来把这个问题整理成了一个便于计算的统计量。先对数据做中心化处理令(X_i x_i - \bar{x})其中(\bar{x})是均值。然后对于每个频率(f)定义一个与时间偏移有关的相位量(\tau)[ \tan(2\omega \tau) \frac{\sum_i \sin(2\omega t_i)}{\sum_i \cos(2\omega t_i)} ]其中(\omega 2\pi f)。这个偏移的目的是让正弦和余弦基函数在采样时刻上保持正交性从而让功率计算变得稳定。最终的周期图功率定义为[ P(f) \frac{1}{2\sigma^2}\left{ \frac{\left[\sum_i X_i \cos\omega(t_i-\tau)\right]^2}{\sum_i \cos^2\omega(t_i-\tau)} \frac{\left[\sum_i X_i \sin\omega(t_i-\tau)\right]^2}{\sum_i \sin^2\omega(t_i-\tau)} \right} ]其中(\sigma^2)是数据的方差。这个式子的计算复杂度看起来是(O(N))每个频率如果我们扫几千个频率点计算量就是(O(N \times M))。对现代计算机来说几万条数据配几千个频率点完全不是问题。2.2 和FFT的关系到底谁更“准”有人可能会问既然傅里叶变换在等间隔采样下运行极快(O(N\log N))为什么非均匀采样不直接套FFT偏要用Lomb-Scargle这种逐频率拟合的笨办法关键在于FFT的数学推导依赖一个前提采样点均匀分布在时间轴上。一旦这个前提被打破FFT的基函数之间不再正交计算出来的频谱中会混入各频率之间的“串扰”。用插值硬凑均匀只是牺牲分辨率换来了一个“能跑的FFT”而不是真正准确的频谱。Lomb-Scargle的拟合过程虽然慢一些但它在每个频率上做的是真正的投影不会因为其他频率成分的存在而干扰当前频率的功率估计。尤其是当数据存在大量缺口、采样间隔差异明显时Lomb-Scargle往往能恢复出真实的谱峰而插值FFT会在谱峰旁边拖出一长片泄漏。当然如果数据本身已经是均匀采样的FFT无论在速度还是精度上都是更好的选择Lomb-Scargle只是为“非均匀”这个场景兜底的方案。3. MATLAB实操用plomb函数实现周期估计3.1 数据准备先模拟一组非均匀采样信号在MATLAB里最简单的落地方案是直接用信号处理工具箱自带的plomb函数。它的名字虽然不叫lombscargle但底层就是Lomb-Scargle周期图的标准实现。先用一个仿真例子跑通流程方便说明参数设置。rng(42); % 构造均匀时间轴目标频率5 Hz t_uniform 0:0.01:60; f_true 5; x_clean sin(2*pi*f_true*t_uniform) 0.3*randn(size(t_uniform)); % 随机抽取20%的点得到非均匀采样序列 keep rand(size(t_uniform)) 0.2; t sort(t_uniform(keep)); x x_clean(keep);这里故意只保留20%的点相当于采样率很低而且时间间隔完全不均匀。如果用插值后FFT插值点的密度不够结果会很糟糕但Lomb-Scargle可以直接利用这些不规则的采样时刻。画个散点图先看一眼数据形态figure; stem(t, x, filled, MarkerSize, 3); xlabel(Time (s)); ylabel(Amplitude); title(Non-uniformly sampled data);你会看到时间轴上点分布集中在某些区域、又缺掉某些区域这就是典型的非均匀采样。3.2 调用plomb并解读输出调用方式非常简单[p, f] plomb(x, t); plot(f, p); xlabel(Frequency (Hz)); ylabel(Power); xlim([0, 10]);plomb默认会做去均值处理所以即使信号有直流分量也不会在0频处出现大尖峰。输出p是功率谱密度f是对应的频率轴单位取决于输入时间戳的单位。运行后你应该能在5 Hz附近看到一个明显的尖峰这正是我们埋进去的真实频率。这里有个重要经验如果数据里混有线性趋势光靠plomb默认的去均值是不够的。比如传感器的缓慢漂移会在低频段形成一个“陡坡”把真正的周期峰淹没。所以实际操作中我习惯于先对数据做一次线性去趋势x_detrend detrend(x, linear); [p, f] plomb(x_detrend, t);去趋势之后低频段的能量会明显减小中高频的周期峰才能浮现出来。这个细节不处理很多时候你会看到周期图整体呈“右下倾斜”的形状第一个尖峰永远出现在最低频处其实那个不代表真实周期只是趋势的贡献。3.3 频率范围与结果解读plomb默认会基于数据的时间跨度自动选择频率范围。不过在某些场景下自动范围可能包含过高的频率导致计算量偏大且产生大量伪峰。我们可以手动指定频率范围freq linspace(0, 10, 2000); [p, f] plomb(x_detrend, t, freq);这里的freq实际上是一个频率网格plomb会在这些频率点上计算功率。网格越密频率分辨率越细但计算量线性增加。我的经验是先根据物理背景确定“可能存在的周期范围”再在这个范围内均匀取1000~5000个频率点就足够用了。拿到周期图后找峰通常用findpeaks[pks, locs] findpeaks(p, f, SortStr, descend, NPeaks, 3);但要注意一点findpeaks找出来的峰值不一定都显著。这时需要评估“误报概率”False Alarm Probability, FAP。plomb内部提供了相关计算具体参数可以查MATLAB文档中关于FAP/Pd的选项。核心思想是如果某峰值的高度远高于噪声基底那它偶然出现的概率就极低对应FAP很小说明周期可信。不要只看峰的高度差要结合统计显著性判断尤其是低信噪比数据。4. 自己编写lombscargle.m时的关键细节4.1 去均值、去趋势与归一化虽然直接用plomb很方便但有时候我们得在旧版MATLAB或者没有信号处理工具箱的环境下工作这时自己写一个lombscargle.m就很有必要。我在自编时踩过几个坑第一个就是数据中心化。如果数据有一个较大的直流偏置而公式中直接用原始(x_i)去计算正弦拟合的分子计算结果会在低频段出现巨大峰值因为直流分量对任意低频正弦都有很强的投影。正确的做法是先把均值减掉变成(X_i x_i - \bar{x})。更进一步如果数据存在缓慢漂移最好先做一次线性去趋势否则周期图低频部分会被趋势主导。我习惯把函数入口设计成function [p, f] lombscargle(t, x, fgrid) % LOMBSCARGLE 计算非均匀采样时间序列的Lomb-Scargle周期图 % 输入 % t - 时间戳列向量 % x - 观测值列向量长度与t相同 % fgrid - 频率网格用于计算的频率点 % 输出 % p - 周期图功率 % f - 与p对应的频率4.2 频率网格过采样因子和最高频率的选择自编时最影响结果质量的是频率网格的选择。很多文献里提到两个参数ofacoversampling factor和hifachighest frequency factor。ofac用来控制频率网格的过采样程度通常取值4~20hifac用来控制最大扫描频率通常取平均采样间隔倒数的一半乘以某个系数。没有现成工具箱的话可以按下面逻辑设定频率上限% 平均采样间隔 dt_mean mean(diff(t)); % 一个保守的频率上限 fmax 1 / (2 * dt_mean); % 频率网格步长取单个频点对应的扫描宽度 df 1 / (range(t) * 4); fgrid (0:df:fmax).;这里range(t)是数据的总时间跨度df相当于在傅里叶分析中“频率分辨率”的四分之一目的是让谱峰更平滑便于找峰。如果数据里有很短的采样间隔用1/(2*dt_min)算出来的上限会非常高导致计算量爆炸。实际操作时我通常会结合信号物理特征把fmax限制在一个合理范围内比如机械故障诊断只看0~500 Hz天文周期搜索只扫0~1/天没必要扫描所有理论可用的频率。4.3 自编实现的计算效率窍门最朴素的实现是双层循环外层扫频率内层对每个频率做正弦拟合。这在数据量只有几百、频率点只有几百时没问题但数据量一大就非常慢。推荐的做法是向量化把时间向量写成行向量频率网格写成列向量然后用矩阵运算一次性算出所有频率下的正弦和余弦值。以核心的分子计算为例可以这样向量化omega 2 * pi * fgrid; % 频率网格角频率列向量 T t(:).; % 时间行向量 cos_wt cos(omega .* T); % 矩阵每个频率×每个时间点 sin_wt sin(omega .* T); % 重心化后的数据 xw x(:) - mean(x); % 计算每个频率下的余弦基投影 cos_sum cos_wt * xw; % 相当于 sum(xw * cos(omega_i * t_k))不过要注意Lomb-Scargle公式中还有一个相位偏移(\tau)也就是说直接用cos(omega .* T)和sin(omega .* T)做投影在非均匀采样时存在基函数不正交的问题。简单实现可以忽略(\tau)直接算得到的结果和标准版有细微差别严格实现需要先对每个频率求出(\tau)再做一次旋转。我建议初学者先用plomb对比验证自编结果如果差异在可接受范围内可以接受简化版毕竟很多实际判断只看峰位不看绝对功率。5. 常见问题与避坑经验5.1 伪峰从哪里来非均匀采样下最大的坑是“混叠”。由于采样间隔不均匀某些频率的真实信号可能会在另一个频率处产生虚假的响应。这有点像街头摄影中飞驰而过的车轮拍到的轮辐看起来转得很慢甚至反向转因为你采样时刻刚好“欺骗”了视觉。要降低伪峰风险有几个实操手段。一是尽量保证数据中缺失的块不要太大不要出现长时间完全无观测的“空洞”二是多个频率网格参数交叉验证在ofac4和ofac16下都算一遍真实峰的位置应该稳定不动伪峰则经常随网格变化三是利用FAP筛选峰值。三者结合基本能避免把随机波动当成周期。5.2 数据长度与信噪比阈值Lomb-Scargle统计量在数据量较小时有较高的不确定性。如果总观测时长只覆盖了目标周期的两三个周期即使数据确实存在周期性谱峰也会很宽、很难精确定位。我的经验是至少要让时间跨度覆盖目标周期的5倍以上峰位才能可靠。假如你的数据只有60秒就别指望能从里面可靠提取周期大于30秒的信号。另外噪声不能太大。一个粗略的参考是如果数据信噪比低于3 dB周期图上的峰就可能淹没在噪声中肉眼几乎看不出。这时候不要盲目相信plomb的输出而是尝试对数据做带通滤波预处理或者分段后重新计算。5.3 工具箱依赖与版本问题plomb函数需要Signal Processing Toolbox如果用的是轻量级MATLAB授权或者手头只有MATLAB Online基础版有可能调用时报错“未定义函数或变量”。这时候自写的lombscargle.m就派上用场了。还有一个坑是频率单位如果时间戳t的单位是小时那么plomb返回的频率单位是“周/小时”不是默认的Hz。我一度因为时间单位用错把算出来的周期读成了完全不对的秒数。建议在代码开头就把时间统一换算成秒后面所有频率解读都基于单位Hz。最后再分享一个我自己调试的小习惯拿到一批非均匀采样数据后我会先用一个已知周期的人工合成信号做一遍全流程测试确认代码、参数、解读方式没问题再处理真实数据。这样做的好处是你对结果长什么样心里有数真数据跑出奇怪结果时能更快判断是算法问题、参数问题还是数据本身的问题。非均匀采样的周期分析并不神秘关键的思路就是尊重每个采样点的真实时刻用拟合代替插值这也是lombscargle.m这类代码能站住脚的根本原因。本文还有配套的精品资源点击获取