基于ARMA模型的高分辨率谱分析工具箱

资 源 简 介

该项目是一个专为MATLAB环境开发的自回归滑动平均（ARMA）谱分析工具箱，旨在通过参数化模型方法对时间序列数据进行高分辨率的功率谱密度（PSD）估计。主要功能模块包括：1. 参数估计算法实现，利用扩展尤尔-沃克（Modified Yule-Walker）方程、最小二乘法或最大似然估计等技术，从观测数据中精确提取AR（自回归）和MA（滑动平均）模型系数；2. 谱密度计算，基于计算出的模型参数和方差估计，生成平滑且高分辨率的功率谱，特别适用于数据长度较短且频率特征复杂的信号分析，弥补了传统FFT周期图法在短数据记录下分辨率不足的问题；3. 模型阶数选择，内置赤池信息量准则（AIC）、贝叶斯信息量准则（BIC）等指标，自动辅助用户确定最佳的AR和MA阶数，以防止过拟合或欠拟合；4. 极零点分析，提供系统传递函数的极点和零点分布图，帮助分析信号的稳定性和频率响应特性；5. 可视化接口，提供一键式绘图功能，支持不同阶数模型谱估计结果的对比展示。该工具箱广泛适用于生物医学工程（如脑电信号分析）、地震数据处理、机械故障诊断以及金融时间序列分析等领域，能够有效提取淹没在宽带噪声中的窄带信号分量。

详 情 说 明

AutoRegressive Moving Average Spectral Analysis Toolbox

项目简介

本项目是一个专为MATLAB环境开发的自回归滑动平均（ARMA）谱分析工具箱。旨在解决短数据记录下传统FFT周期图法分辨率不足的问题，通过参数化模型方法实现高分辨率的功率谱密度（PSD）估计。该工具箱特别针对淹没在宽带噪声中的窄带信号（如靠得很近的正弦波）进行了优化，集成了信号生成、模型阶数自动选择、参数估计、谱计算及多维度可视化分析功能。

功能特性

• 高分辨率谱估计：基于ARMA模型参数计算功率谱，有效分辨频率间隔极小的信号分量（如演示中的150Hz与155Hz）。
• 混合参数估计算法：

* AR部分：采用扩展尤尔-沃克（Modified Yule-Walker）方程，利用高阶自相关函数避免MA参数的干扰。 * MA部分：采用基于残差的迭代最小二乘法（Least Squares）进行拟合。

• 智能阶数选择：内置模型定阶扫描功能，计算赤池信息量准则（AIC）和贝叶斯信息量准则（BIC），通过热力图辅助用户选择最优模型阶数。
• 多维度可视化：

* 时域波形观察。 * ARMA谱与经典周期图法的叠加对比（含局部放大）。 * 系统极零点分布图（P-Z Plot）。 * BIC定阶准则热力图。

• 合成信号仿真：内置测试信号生成器，可模拟含高斯白噪声的双频窄带信号。

系统要求

• MATLAB R2016b 或更高版本。
• 无需额外的信号处理工具箱（核心算法均为原生实现），但绘图部分依赖基本MATLAB图形库。

使用方法

1. 确保工作路径包含所有相关脚本。
2. 直接运行主程序入口函数。
3. 程序将自动执行以下流程：

* 生成模拟信号。 * 执行AIC/BIC阶数扫描。 * 基于指定阶数（默认为AR=12, MA=8）进行参数估计。 * 弹出包含四个子图的综合分析窗口。

1. 观察控制台输出的推荐阶数以及图形窗口的谱估计结果。

---

核心代码逻辑与算法分析

以下基于代码实际实现内容，对各个模块的技术细节进行深入解析：

1. 信号生成与预处理

程序首先构建了一个具有挑战性的测试环境。

• 采样参数：设定采样率Fs为1000Hz，时长2秒。
• 信号成分：合成两个频率极为接近的正弦波（150Hz和155Hz），模拟难以分辨的窄带信号。
• 噪声环境：叠加了高强度的高斯白噪声（噪声水平系数1.5），测试算法在低信噪比下的鲁棒性。
• 去均值处理：在分析前移除了信号的直流分量，防止零频分量干扰谱估计。

2. 模型阶数选择机制

代码通过双重循环遍历AR（p）和MA（q）的阶数范围。

• 扫描范围：AR阶数从1到15，MA阶数从1到10。
• 评价准则：

* AIC（赤池信息量准则）：侧重于预测精度，计算公式体现为 N*ln(v) + 2(p+q)。 * BIC（贝叶斯信息量准则）：对模型复杂度惩罚更重，计算公式体现为 N*ln(v) + (p+q)ln(N)。

• 热力图输出：最终将BIC矩阵以热力图形式展示，颜色越深（数值越小）代表模型阶数匹配度越高。

3. ARMA 参数估计算法

这是工具箱的核心数学引擎，实现了分步估计策略：

第一步：AR参数估计（Modified Yule-Walker）

• 为了解决ARMA模型中AR和MA参数耦合的问题，代码利用了高阶自相关函数的特性。
• 原理：对于ARMA(p,q)过程，其自相关函数在滞后数大于q时，主要由AR部分决定。
• 实现：构造扩展的Toeplitz矩阵（R matrix）和相关向量（r vector），利用滞后范围在 [q+1, q+p] 的样本自相关函数建立方程组 R * a = -r。
• 求解：使用伪逆（pinv）求解线性方程组，提高了矩阵病态时的数值稳定性。

第二步：MA参数估计（迭代最小二乘法）

• 残差获取：利用第一步估计得到的AR系数对原始信号进行逆滤波，获得近似的MA过程残差序列。
• 参数拟合：代码实现了一个简化的迭代高斯-牛顿（Gauss-Newton）风格的求解器。

* 它构建回归矩阵M，尝试寻找MA系数向量b，使得 M * b 逼近残差序列。 * 通过少量迭代（默认为5次）更新误差估计及系数，从而提取MA部分的特征。

• 方差估计：最终计算拟合残差的方差作为白噪声激励的方差估计值。

4. 功率谱密度（PSD）计算

代码没有仅依赖FFT对数据进行变换，而是利用估计出的模型参数计算解析谱。

• 传递函数：构建系统函数 H(z) = B(z) / A(z)，其中B(z)为MA部分多项式，A(z)为AR部分多项式。
• 频率响应：分别计算分子和分母在单位圆上的频率响应。
• 谱与频率变换：计算公式为 sigma^2 * |B(w)|^2 / |A(w)|^2。为了便于对比，将归一化频率映射回实际物理频率（Hz），并进行了相应的幅度缩放。

5. 结果可视化

绘图模块生成一个包含四个子图的窗口：

• 时域图：展示前300个采样点，直观显示含噪信号形态。
• PSD对比图：将经典的周期图法（灰色曲线）与ARMA参数法（红色曲线）叠加。图中包含一个局部放大插图（Inset），专门展示140Hz-165Hz频段，证明ARMA方法能清晰分辨出周期图法无法区分的150Hz和155Hz双峰。
• 极零点图：绘制单位圆内的极点和零点位置，用于分析系统的稳定性及频率选择特性。
• 定阶热力图：直观展示不同(p, q)组合下的BIC值分布，原点置于左下角。