简介一份基于分数低阶统计量与低阶循环平稳理论的 MATLAB 实现代码面向从事信号处理、阵列信号与波达方向DOA估计研究的学生和工程师重点解决脉冲噪声环境下传统循环平稳方法性能退化的问题。包内共 4 个 m 文件压缩包仅 2KB包含主程序 FLOM-TLS-Cyclic-ESPRIT1.m 以及 csd2.m、mse.m、stable (2).m 等辅助脚本涵盖分数低阶矩计算、循环互谱密度估计、误差评估与稳定性分析可支撑从算法仿真到精度验证的完整流程。目前已有 293 人学习下载适合对非高斯信号处理、Alpha 稳定分布噪声建模及循环平稳 ESPRIT 类算法感兴趣的中高级研究者。通过运行该代码可直观理解分数低阶循环平稳统计量在强脉冲干扰下的鲁棒性并借助模块化脚本自行扩展参数配置、对比不同信噪比下的估计性能为后续算法改进或论文复现提供便利。1. 为什么脉冲噪声下 ESPRIT 会失效FLOC-ESPRIT 要解决什么实测阵列信号处理有个反直觉现象发射功率提高 10 dBESPRIT 的角度估计误差反而上升。原因不在接收机饱和而在噪声模型——大气放电、舰船辐射和人为干扰产生的是脉冲噪声用高斯拟合会严重低估尾部用对称 α 稳定分布SαS更真实但这种分布不存在有限二阶矩协方差矩阵的数学基础在样本估计层面直接崩塌。FLOC-ESPRIT 的思路是保留 ESPRIT 的旋转不变框架把二阶矩替换成分数低阶矩Fractional Lower Order Statistics再借助循环平稳信号的循环频率把目标从脉冲背景里筛出来这就是标题里“循环分数低阶”的含义。下面按数学模型、MATLAB 实现、参数调节、验证技巧的顺序给出一套能直接改参数跑通的做法。2. 从协方差到分数低阶循环平稳FLOC-ESPRIT 的数学地基2.1 为什么 α 稳定分布下二阶矩失效FLOC 的数学前提SαS 分布的形状由特征指数 α ∈ (0, 2] 决定。α 2 时就是高斯分布α 越小分布尾部越厚出现大幅度脉冲样本的概率越高。工程里常见的海洋环境噪声、大气噪声、民用频段的人为干扰实测统计常常落在 α 1.2 到 1.8 之间。对这样的分布做二阶矩运算会遇到本质问题当 α 2 时E[|X|²] 趋于无穷。也就是说理论上协方差矩阵根本不存在样本协方差矩阵等于拿一组“方差无限”的随机数求平均每次新到的脉冲样本都会把矩阵元素拉出一个极端值。这解释了开头那个反直觉现象提高发射功率同时噪声脉冲的绝对幅度也成比例变大传统 ESPRIT 的协方差估计被脉冲拖垮估计误差不降反升。SαS 分布存在有限 p 阶矩的条件是 p α。这是整个分数低阶方法的数学支点选定阶数 p 在噪声特征指数 α 之下样本矩估计在统计意义上收敛脉冲对矩的贡献被非线性压缩。FLOC 本质上是用“更低阶的矩”换“对尾部有界的估计”代价是矩阶数降低后估计方差变大需要更多快拍弥补这个取舍后面第 4 章会具体展开。2.2 分数低阶矩定义与 FLOC 矩阵把 R 换成 Φ经典 ESPRIT 的出发点是协方差矩阵 R E[x x^H]。FLOC-ESPRIT 的常见做法是把 R 的每个元素替换成分数低阶相关形式。对均匀线阵接收数据 x(t) [x₁(t), …, x_K(t)]ᵀ定义 FLOC 矩阵 Φ 的第 (i, j) 个元素为Φ(i, j) E[ x_i(t) · conj(x_j(t)) · |x_j(t)|^(p-2) ]其中 1 p α。对照协方差矩阵的定义可以看到FLOC 矩阵在结构上只改了一处把 x_j 的共轭后面乘了一个幅度加权项 |x_j|^(p-2)。当 p 2 时该项恒等于 1FLOC 矩阵就退化成协方差矩阵所以高斯环境下传统 ESPRIT 可以视为 FLOC-ESPRIT 在 p 2 时的一个特例。这个加权项是 FLOC 抑制脉冲的关键。某个阵元在某一时刻被大幅值脉冲击中时|x_j|^(p-2) 会随脉冲幅度增长而急剧衰减把该样本在矩累积中的权重压下去。换句话说FLOC 矩阵给每个样本的贡献乘了一条随幅度变化的权重曲线脉冲样本的权重大幅降低正常信号样本保留下来。用 MATLAB 计算时一个快拍矩阵 xK×N对应的权重矩阵可以写成reg_delta 1e-6; % 避免 |x|0 时出现除零 w abs(x).^(p-2); % 逐元素幅度加权形状 K×N w max(w, reg_delta); % 底部截断保证数值稳定说明这里的p-2在 1 p 2 条件下是负数所以|x|越大的样本权重越小。max(w, reg_delta)是工程处理防止某个采样点恰好接近零时0^(p-2)溢出成 Inf这在低快拍或稀疏信号场景下并不罕见。需要注意的是按上面定义直接累积得到的 Φ 不保证 Hermitian 对称。ESPRIT 后续要做特征分解非 Hermitian 阵的特征向量不一定正交信号子空间提取会不稳定。常见做法是累积完做一次对称化Φ (Φ Φᵀ)/2。2.3 循环平稳如何加第二道保险循环 FLOC 矩阵的样本估计FLOC 只解决了脉冲噪声的重尾问题如果目标信号本身功率弱平稳噪声的宽带能量依然会淹没它。循环平稳思想在这里补上了第二道筛选真实通信信号通常具有循环平稳性其统计量随时间周期性变化在特定循环频率处存在谱相关峰而平稳噪声在任何非零循环频率处都没有相关性。把 FLOC 和循环平稳结合就得到循环 FLOC 矩阵Φ_ε (1/N) Σ_{t1}^N x(t) · [conj(x(t)) ⊙ |x(t)|^(p-2)]ᵀ · exp(-j2π ε t)其中 ε 是归一化循环频率。累加项中的复指数起到了“频率选择”的作用只有统计特性以 ε 为周期变化的信号会在这个方向上累积出非零结果平稳噪声和循环频率不匹配的其他信号会被平均掉。实现这个式子时要注意一个细节循环频率的单位要与数据的时间索引一致。如果 MATLAB 代码里时间轴 t 是采样点序号 0:N-1那么循环频率 ε 的单位是弧度/采样点。BPSK、QPSK 这类调制信号的循环频率与符号率有关实际工程中通常先用循环谱估计出符号率再换算成弧度/采样点。3. 用 MATLAB 实现循环 FLOC-ESPRIT主函数与完整仿真脚本3.1 仿真参数与阵列信号模型ULA 下 M 个信源先设定一组可复现的仿真参数后续所有代码都基于这套参数运行。采用 8 元均匀线阵阵元间距取半波长两个等功率信源分别位于 -20° 和 35°噪声用特征指数 α 1.5 的 SαS 分布生成广义信噪比 GSNR 定为 5 dB。快拍数取 2000这是低阶统计量收敛的基本量级。K 8; % 阵元数 d_lambda 0.5; % 阵元间距/波长 theta_true [-20, 35]; % 真实波达角单位度 M length(theta_true); % 信源数 N 2000; % 快拍数 alpha_noise 1.5; % SαS 噪声特征指数 GSNR_dB 5; % 广义信噪比 p 1.2; % 分数低阶阶数必须 alpha_noise fc_over_fs 0.1; % 循环频率单位 cycles/sample eps0 2 * pi * fc_over_fs; % 换算成弧度/采样点说明fc_over_fs模拟的是数字基带信号里残余载波或符号率分量的归一化频率。复基带数据乘以exp(j2π·0.1·t)后其二阶统计量会在 0.1 cycles/sample 处出现循环平稳特征EPS0 对应 FLOC 矩阵累加时用的复指数频率。实际应用中这个值应该由循环谱估计得到而不是手工指定。阵列流型矩阵 A 的第 k 行第 m 个元素为 exp(j2π·d_lambda·k·sin(θ_m))k 从 0 开始。用 MATLAB 生成A exp(1j * 2 * pi * d_lambda * (0:K-1). * sind(theta_true));这里(0:K-1).是列向量sind(theta_true)是行向量外积得到 K×M 矩阵。上式假设了阵列第一个阵元在原点相位参考点选在阵元 0这个约定在后面对 ESPRIT 的旋转不变关系做推导时会用到。3.2 生成 SαS 脉冲噪声和循环平稳信号sas_noise 函数统计概率分布 matlab 里的标准做法是自行实现 Chambers-Mallows-Stuck 算法不依赖 Communications Toolbox。function z sas_noise(alpha, gamma, dim) % z sas_noise(alpha, gamma, dim) % 生成实值对称 α 稳定分布随机数 % alpha: 特征指数 (0,2] % gamma: 尺度参数对应分散度 if alpha 2 z sqrt(2 * gamma) * randn(dim); % 高斯退化情形 else V pi * (rand(dim) - 0.5); % 在 (-π/2, π/2) 均匀分布 W -log(rand(dim)); % 指数分布均值 1 z sin(alpha * V) ./ (cos(V).^(1/alpha)) ... .* (cos((1-alpha) * V) ./ W) .^ ((1-alpha) / alpha); z gamma^(1/alpha) * z; end end说明V和W是两个服从标准分布的中间随机变量组合后得到标准 SαS 随机数。gamma是尺度参数gamma^(1/alpha)的作用是把标准分布拉伸到目标分散度。当 alpha 接近 1 时指数(1-alpha)/alpha接近 0数值上要注意避免负数的分数次幂好在 matlab 里cos((1-alpha)*V) ./ W在 alpha1 时可能出现负数实际使用中建议把 alpha 下限限制在 1.05 以上或者直接取绝对值加重采样处理。信号侧用一个复基带随机序列乘上残留载波构成循环平稳信号。之所以不直接用 randn 生成的平稳序列是因为纯平稳信号在非零循环频率处没有谱相关峰FLOC-ESPRIT 的循环项会把它当成干扰滤掉这会让算法“看不见”目标。base sign(randn(M, N)) 1j * sign(randn(M, N)); % 复随机序列 t 0:N-1; s base .* exp(1j * 2 * pi * fc_over_fs * t); % 乘残留载波 sig A * s; % 阵列输出参数说明base每行是一个信源的复包络sign(randn(...))生成 ±1±1j 的离散取值可以理解为简化的 BPSK/QPSK 符号。乘上复指数后信号的循环频率精确等于fc_over_fs方便后续比对。3.3 循环 FLOC-ESPRIT 主函数核心实现与角度提取主函数把 FLOC 矩累积、循环频率匹配、SVD 和 ESPRIT 旋转不变关系串起来。下面是可直接复制的函数核心代码。function theta_est floc_esprit(x, M, p, eps0, d_lambda) % FLOC-ESPRIT: 低阶循环平稳阵列角度估计 % 输入: % x K×N 复基带数据 % M 信源数 % p 分数低阶阶数1 p alpha_noise % eps0 循环频率弧度/采样点 % d_lambda 阵元间距/波长 % 输出: % theta_est 1×M 角度估计值单位度 [K, N] size(x); reg 1e-6; % 数值保护项 % 1. 幅度权重矩阵冲击样本被自动降权 w max(abs(x).^(p-2), reg); % 2. 循环频率相位项 phase exp(-1j * eps0 * (0:N-1)); % 3. 累积循环 FLOC 矩阵 C zeros(K, K); for n 1:N xn x(:, n); wn w(:, n); C C phase(n) * (xn * (conj(xn) .* wn).); end C C / N; C 0.5 * (C C); % 强制 Hermitian稳定信号子空间 % 4. SVD 取信号子空间 [U, S, ~] svd(C); Us U(:, 1:M); % 5. ESPRIT 旋转不变关系 U1 Us(1:end-1, :); U2 Us(2:end, :); psi pinv(U1) * U2; [~, D] eig(psi); % 6. 相位转角度 theta_est asin(-angle(diag(D)) / (2 * pi * d_lambda)); theta_est sort(theta_est. * 180 / pi); end代码逻辑说明第 1 步算出的权重矩阵w形状与x完全相同conj(xn) .* wn对应公式里的 conj(x_j)·|x_j|^(p-2)xn * (…).完成向量外积。theta_est 里的负号取决于阵元编号方向和 ESPRIT 投影方向如果你在自己数据上发现角度符号整体相反直接删掉负号即可这是最常见的移植坑。单次运行不一定能精确复现角度因为 SαS 噪声方差无限、单次实现的随机性很大。正确的验证方式是固定角度和噪声特征指数在多个 GSNR 点做 Monte Carlo 仿真看统计平均下的 RMSE这部分第 5 章给出脚本。3.4 顶层仿真脚本GSNR 与噪声尺度对齐把上面所有函数串起来的顶层脚本如下其中混合信噪比需要做显式对齐否则 GSNR 设定没有意义。rng(20240601); % 生成信号和阵列流型 A exp(1j * 2 * pi * d_lambda * (0:K-1). * sind(theta_true)); base sign(randn(M, N)) 1j * sign(randn(M, N)); s base .* exp(1j * 2 * pi * fc_over_fs * t); sig A * s; % 按 GSNR 调整信号幅度 gamma_noise 1; % 噪声分散度 P_sig mean(abs(sig(:)).^2); s s * sqrt(10^(GSNR_dB/10) * gamma_noise / P_sig); x A * s sas_noise(alpha_noise, gamma_noise, [K, N]); % 调用 FLOC-ESPRIT theta_hat floc_esprit(x, M, p, eps0, d_lambda); fprintf(真实角度: %s\n, mat2str(theta_true)); fprintf(估计角度: %s\n, mat2str(theta_hat, 4));GSNR 与高斯噪声下的 SNR 含义不同它把噪声功率用分散度 γ 替代因为 SαS 噪声没有有限方差不能用“噪声方差”做归一化。这里sqrt(10^(GSNR_dB/10)·γ/P_sig)保证信号平均功率与 γ 的比例关系满足设定值。4. FLOC-ESPRIT 参数怎么调阶数 p、循环频率与快拍数的取舍4.1 分数低阶阶数 p 的选取1 p α 是硬约束p 的取值同时受理论约束和工程权衡限制。理论上必须满足 p α_noise否则矩本身发散样本估计不收敛。工程上 p 越小对脉冲的压制越强但加权过程丢失的有用信息也越多估计方差变大p 越接近 α越接近二阶矩行为方差变小但脉冲残余变大偏差和抖动同时上升。我习惯的做法是取 α 的 60% 到 80% 作为 p 的初始值再结合 GSNR 微调。低 GSNR 场景脉冲占主导p 取下限高 GSNR 场景信号质量好p 可以靠近 α换取更小的角度方差。噪声特征指数 α推荐 p 区间GSNR 5 dB 建议GSNR 10 dB 建议1.21.05 ~ 1.151.05 ~ 1.101.10 ~ 1.151.51.10 ~ 1.401.10 ~ 1.201.30 ~ 1.401.81.30 ~ 1.701.30 ~ 1.451.55 ~ 1.70p 选错时最典型的症状是角度谱上出现大量虚假峰值或者特征是 SVD 的特征值差距拉不开。遇到这种情况先把 p 降到 1.1 跑一遍确认“能估计”之后再逐步增大 p 观察 RMSE 是否改善而不是一开始就把 p 定在接近 α 的位置。4.2 循环频率 eps0 的匹配误差与粗估计方法循环频率不匹配时FLOC 矩阵累加项里的复指数会把目标信号旋转平均掉等效于目标功率被|sinc(Δε·N/2)|衰减。Δε 越大衰减越剧烈。工程上常见的错误是把符号率直接当作循环频率但实际上 BPSK 的谱相关峰出现在符号率的两倍处QPSK 则为符号率整数倍位置换算时必须先确认调制类型。另一种常见情况是接收机存在残余载波频偏。此时循环频率等于残余频偏除以采样率可以从单阵元数据的幅度平方谱粗估。下面给出一种不依赖循环谱工具箱的粗估做法y abs(x(1, :)).^2; % 单阵元幅度平方 Y fftshift(fft(y)); f_axis (-N/2:N/2-1) / N * 2 * pi; % 弧度/采样点 [~, idx] max(abs(Y)); eps0_est abs(f_axis(idx));注意幅度平方处理会产生二倍频效应因此估计值可能是真实循环频率的两倍。粗估结果建议在eps0_est/2、eps0_est、2*eps0_est三个候选值上分别跑一次 FLOC-ESPRIT让计算机来选择 RMSE 最小的那个。这个方法虽然笨但在工程现场比反复调频率参数快得多。4.3 快拍数 N 与阵元数 K 的下限选择低阶统计量的样本估计收敛速度比二阶矩慢。高斯环境下 200 个快拍可以工作的场景切到 SαS 噪声后通常需要 1000 到 2000 个快拍才能稳定。快拍量不足时FLOC 矩阵的噪声子空间特征值不会出现明显跌落M 选得再准确也会在角度谱上看到毛刺。阵元数的约束来自 ESPRIT 的旋转不变构造。信号子空间取前 M 列后U1的尺寸是 (K-1)×M要求 K-1 ≥ M为了给噪声子空间留余量实际工程至少取 K ≥ 2M。快拍数和阵元数的匹配建议如下信源数 M最小阵元数 K推荐快拍数 N261000 ~ 20004102000 ~ 40008164000 ~ 8000当快拍数明显不足时与其继续增大 p不如先用分块滑窗累积 FLOC 矩阵再做特征分解。这部分在第 5.3 节展开。5. 验证 FLOC-ESPRIT 的边界特征值谱、Monte Carlo 与低复杂度实现5.1 特征值谱判断信源数是否给对循环 FLOC 矩阵的 SVD 特征值谱与经典协方差的特征值谱形态类似但脉冲噪声下高低特征值之间的差距会被压缩。建议先做一次 SVD 并打印特征值衰减曲线确认信号子空间和噪声子空间的分界。[~, S, ~] svd(C_eval); sv_dB 10 * log10(diag(S)); stem(sv_dB, filled); grid on; xlabel(特征值序号); ylabel(幅度 (dB)); title(循环 FLOC 矩阵特征值谱);操作要点找到曲线斜率由陡变缓的拐点拐点前的特征值个数就是可信的信源数。脉冲噪声较强时拐点不明显可结合 MDL 准则在特征值序列上做粗判但千万不要直接套用高斯假设下的阈值因为 FLOC 矩阵的特征值分布没有解析闭式MDL 只能指示数量级。5.2 Monte Carlo 验证RMSE vs GSNR 曲线要说服自己算法“真的有效”必须做多次独立实验统计 RMSE。下面脚本在多个 GSNR 点上重复 200 次统计估计角度与真实角度的均方根误差。MC 200; GSNR_set -5:2:15; rmse zeros(size(GSNR_set)); for gi 1:length(GSNR_set) err_acc 0; for trial 1:MC % 重新生成信号、噪声和接收数据代码同 3.4 节 x A * s sas_noise(alpha_noise, gamma_noise, [K, N]); th floc_esprit(x, M, p, eps0, d_lambda); % 角度误差需做最小距离匹配避免顺序颠倒和 ±180° 模糊 diff1 sum((th - theta_true).^2); diff2 sum((th - (theta_true 180)).^2); err_acc err_acc min(diff1, diff2); end rmse(gi) sqrt(err_acc / (MC * M)); end plot(GSNR_set, rmse, -o); xlabel(GSNR (dB)); ylabel(RMSE (度)); grid on;角度匹配是这里最容易被忽略的细节。ESPRIT 的特征值分解结果不会自动按信源顺序重排低快拍时估计角度可能发生两源互换直接相减会造成 RMSE 虚高。最小距离匹配虽然不是最优的指派策略但对于两个信源已经足够。若扩展到 M 2应该改成匈牙利算法处理排序问题。5.3 低复杂度实现滑窗累积与截断 SVDFLOC 矩阵累积是 O(NK²) 的计算快拍数大到 10⁵ 时逐样本循环会明显拖慢实时处理。常见做法是把累加改成滑窗形式每次只更新当前帧mu 0.01; % 遗忘因子 C zeros(K, K); for n 1:N xn x(:, n); wn max(abs(xn).^(p-2), reg); C (1 - mu) * C mu * phase(n) * (xn * (conj(xn) .* wn).); end这个结构可以直接移植到 FPGA 或者用 MATLAB Coder 生成 C 代码每个采样点只需要一次矩阵外积累加不再需要保存整段数据。遗忘因子 mu 对应等效快拍数约 1/mumu 太大则收敛快但方差大mu 太小则跟不上信号循环频率的慢漂移一般取 0.01 到 0.001。最后还有一个实用的回归测试技巧把 p 设成 2噪声换成高斯分布此时循环 FLOC 矩阵退化为循环协方差矩阵FLOC-ESPRIT 的结果应该与经典 ESPRIT 对齐。如果这一步对不上说明代码里的相位符号或角度映射可能出错。确认 p 2 的行为正确之后再把噪声切回 SαS、把 p 调到 0.8α 左右对比角度 RMSE 的改善幅度这才算完整验证了“低阶循环平稳”四个字的实际贡献。本文还有配套的精品资源点击获取