1. 作业背景与核心任务拆解这周的作业题目是“随机信号AR模型及MATLAB实现”看起来是数据压缩课程里的一次实践。很多同学拿到这种题目第一反应可能就是去网上找段代码改改参数把图跑出来交差。但说实话这样学完你可能只知道“AR模型在MATLAB里用aryule或者arburg函数”至于它为什么能用来“压缩”数据它的参数到底代表了什么估计还是一头雾水。我自己当年学信号处理的时候也这么干过后来在实际项目中用到AR模型做预测和特征提取才发现基础没打牢的亏有多大。所以咱们这次不光是完成作业更要把AR模型从原理到实操再到它和数据压缩的关联彻底捋清楚。AR模型全称自回归模型它的核心思想特别直观当前时刻的信号值可以用过去若干个时刻的信号值的线性组合再加上一个随机扰动白噪声来表示。这个“过去若干个时刻”就是模型的阶数。你想想看如果我能用一个简单的线性公式几个系数和一点随机噪声就很好地描述或预测一个信号那我是不是就不用存储原始信号的全部数据点了我只需要存储这几个系数和残差预测误差这不就是一种“压缩”或“简洁表示”吗这就是AR模型在数据压缩尤其是语音、音频编码中的理论基础。本次作业的核心我认为可以拆解为三个层次第一理解AR模型的数学表达和参数意义第二掌握在MATLAB中估计AR模型参数的经典方法比如Yule-Walker方程法、Burg算法第三也是最关键的如何评估你建立的AR模型好不好以及如何直观地展示这个过程。很多教程只教到第二步但第三步才是体现你理解深度的部分。2. AR模型的数学原理不只是几个公式我们先从最根本的数学表达式开始。对于一个p阶的AR模型记作AR(p)它的定义式是[ x[n] -\sum_{k1}^{p} a_k x[n-k] w[n] ]这里x[n]是我们在n时刻观测到的信号值。a_1, a_2, ..., a_p就是我们要求解的AR模型系数也叫反射系数或预测系数。注意公式前面的负号这是信号处理领域的习惯写法MATLAB的相关函数也遵循这个约定。w[n]是驱动白噪声通常假设它均值为0方差为σ²。这个噪声代表了模型无法用过去信号线性解释的部分是信号的“创新”或“随机”成分。这个公式揭示了AR模型的本质它试图用信号自身的过去来预测现在。系数a_k的大小和正负决定了过去第k个时刻的值对当前值的影响程度和方向。举个例子在语音信号中AR模型能非常有效地模拟声道的共振特性那些系数就对应着声道形状的参数。那么如何从一段观测信号x[1], x[2], ..., x[N]中估计出这p个系数a_k和噪声方差 σ² 呢这就引出了几种经典的估计算法它们在MATLAB中都有对应的函数实现。2.1 Yule-Walker方程法这是最经典的方法基于信号的自相关函数。思路是将模型定义式两边同时乘以x[n-m](m0,1,...,p)然后求数学期望。经过一系列推导这里不展开复杂的公式我们可以得到著名的Yule-Walker方程[ R_{xx}[m] -\sum_{k1}^{p} a_k R_{xx}[m-k], \quad m 1,2,...,p ]以及噪声方差[ \sigma^2 R_{xx}[0] \sum_{k1}^{p} a_k R_{xx}[k] ]其中R_xx[m]是信号x[n]的自相关函数在滞后m处的估计值。你看原来求解系数a_k的问题转化为了求解一个以自相关函数值为系数的线性方程组。这个方程组是Toeplitz结构的可以用高效的Levinson-Durbin递归算法来求解。在MATLAB中aryule函数就是基于Yule-Walker方程来估计AR参数的。它的优点是计算稳定理论完备。但缺点是对自相关函数的估计误差比较敏感特别是当数据长度较短时估计的偏差可能会较大。2.2 Burg算法最大熵谱估计Burg算法是一种更现代、通常表现更好的方法。它不直接估计自相关函数而是以前后向预测误差功率最小为准则递推地求解反射系数它和AR系数有确定转换关系。Burg算法能保证产生的AR模型是稳定的即所有极点都在单位圆内并且对于短数据序列其谱估计的分辨率和精度往往优于Yule-Walker法。在MATLAB中arburg函数实现了Burg算法。对于大多数实际应用尤其是本次作业我强烈推荐优先使用arburg。它的估计通常更准确特别是当你处理的信号长度有限时。理解这两种方法的区别很重要它不是简单的“哪个函数更好”而是背后优化准则和计算路径的不同。Yule-Walker法最小化前向预测误差但依赖于自相关函数的估计Burg法则同时最小化前向和后向预测误差直接操作数据。在作业中你可以尝试用两种方法对同一信号建模对比它们的结果这会是一个很好的加分点。3. MATLAB实战从信号生成到模型拟合光说不练假把式我们现在就用MATLAB来走一遍完整的流程。假设我们要分析一个由已知AR模型生成的测试信号这样我们可以将估计出的参数与真实值对比验证方法的有效性。3.1 生成一个测试用的AR信号首先我们得有一个信号。我们可以自己定义一个AR(4)模型然后用它来生成一段数据。这样真实的系数a_k是已知的便于后续对比。% 1. 定义真实的AR(4)模型参数 true_a [1.0, -0.9, 0.5, -0.3]; % 注意这里我们按MATLAB习惯给出的是A(z)1 a1*z^{-1}...的系数 % 即 A [1, a1, a2, a3, a4] true_A [1, true_a]; % 完整的AR多项式系数向量 % 2. 生成驱动白噪声 N 1000; % 信号长度 sigma2 0.1; % 噪声方差 w sqrt(sigma2) * randn(N, 1); % 生成N点的高斯白噪声 % 3. 使用filter函数生成AR信号 % filter(B, A, X) 是用传递函数H(z)B(z)/A(z)对X滤波。 % 对于AR模型B(z)1 A(z)就是我们的AR多项式。 x filter(1, true_A, w); % 生成的信号x现在我们手头有了一段长度为1000点的信号x它是由系数为true_A、噪声方差为0.1的AR(4)模型产生的。我们的任务就是假装不知道true_A仅通过观测信号x来估计出AR模型的阶数p和系数a_k。3.2 使用Burg算法估计AR参数我们先使用推荐的arburg函数。这里有一个关键步骤如何确定阶数p对于测试信号我们知道是4阶。但对于未知信号我们需要定阶。常用方法有AIC赤池信息准则或BIC贝叶斯信息准则它们会在模型拟合优度和复杂度之间做一个权衡。为了简化本次我们先假设通过观察信号的偏相关函数PACF或经验确定了p4。% 使用arburg函数估计AR(4)模型参数 p 4; % 假设我们已经知道或确定了阶数为4 [est_a_burg, sigma2_burg] arburg(x, p); % 注意arburg返回的est_a_burg是[A(2), A(3), ..., A(p1)]即a1, a2, ..., ap % 为了与我们的true_a对比需要构建完整的A向量 est_A_burg [1, est_a_burg];3.3 使用Yule-Walker方法估计AR参数我们也用aryule函数做一次用于对比。% 使用aryule函数估计AR(4)模型参数 [est_a_yule, sigma2_yule] aryule(x, p); est_A_yule [1, est_a_yule];3.4 结果对比与分析现在让我们把真实值和两种方法的估计值放在一起比较。fprintf(真实AR系数: [1, %.4f, %.4f, %.4f, %.4f]\n, true_a); fprintf(Burg估计系数: [1, %.4f, %.4f, %.4f, %.4f]\n, est_a_burg); fprintf(Yule-Walker估计系数: [1, %.4f, %.4f, %.4f, %.4f]\n, est_a_yule); fprintf(\n); fprintf(真实噪声方差: %.4f\n, sigma2); fprintf(Burg估计噪声方差: %.4f\n, sigma2_burg); fprintf(Yule-Walker估计噪声方差: %.4f\n, sigma2_yule);运行这段代码你可能会看到Burg算法估计的系数通常更接近真实值噪声方差的估计也更准。这就是为什么在实践里Burg算法更受青睐。但Yule-Walker法并非一无是处它在理论分析上更直观且对于某些特定类型的信号或足够长的数据表现也不错。4. 模型评估与可视化让结果“说话”把参数估计出来只是第一步。一个负责任的建模过程必须包含模型评估。我们需要回答这个AR(4)模型拟合得好不好有没有可能用更低或更高的阶数4.1 绘制功率谱密度PSD进行对比AR模型的一个重要应用就是谱估计。我们可以用估计出的AR参数来计算信号的功率谱密度并与经典的非参数化方法如周期图法进行对比。% 计算频率向量 NFFT 1024; Fs 1; % 假设采样频率为1Hz频率轴归一化到(0, 0.5) f (0:NFFT/2)*(Fs/NFFT); % 正频率部分 % 1. 使用Burg算法估计的AR参数计算PSD [H_burg, freq_burg] freqz(1, est_A_burg, NFFT/21, Fs); % 计算AR模型的频率响应 PSD_burg sigma2_burg * abs(H_burg).^2; % AR模型的PSD公式σ² / |A(e^{jω})|² % 2. 使用Yule-Walker算法估计的AR参数计算PSD [H_yule, freq_yule] freqz(1, est_A_yule, NFFT/21, Fs); PSD_yule sigma2_yule * abs(H_yule).^2; % 3. 使用周期图法非参数化方法作为参考 [Pxx_periodogram, f_periodogram] periodogram(x, hamming(length(x)), NFFT, Fs); % 绘制在同一张图上进行对比 figure(Position, [100, 100, 900, 600]); subplot(2,1,1); plot(f, 10*log10(PSD_burg), b-, LineWidth, 1.5); hold on; plot(f, 10*log10(PSD_yule), r--, LineWidth, 1.5); plot(f_periodogram, 10*log10(Pxx_periodogram), k:, LineWidth, 0.8); hold off; grid on; xlabel(归一化频率); ylabel(功率谱密度 (dB)); title(AR模型谱估计 vs. 周期图法); legend(Burg算法 (AR), Yule-Walker算法 (AR), 周期图法, Location, Best); xlim([0, 0.5]); % 绘制残差分析图 subplot(2,1,2); % 使用Burg模型对原信号进行预测滤波 x_hat_burg filter([0 -est_a_burg], 1, x); % 注意系数符号这是根据AR定义式做的预测 residual_burg x - x_hat_burg; % 残差 真实值 - 预测值 plot(residual_burg); grid on; xlabel(样本点); ylabel(幅值); title(Burg算法模型残差序列);通过这张图你可以直观地看到谱估计对比AR模型谱特别是Burg算法通常比周期图法平滑得多分辨率也更高。周期图法起伏剧烈方差大而AR谱能清晰地显示出信号的谱峰这对应于AR模型的极点位置。如果两种AR方法估计的谱形状差异很大说明模型估计可能不稳定或者阶数选择不当。残差分析一个拟合良好的AR模型其残差序列应该近似为白噪声即没有明显的自相关结构。你可以通过观察残差序列的波形如上图初步判断更严谨的做法是计算残差的自相关函数ACF并检验其是否在零滞后外接近零。4.2 定量评估残差白噪声检验我们可以用Ljung-Box检验来定量判断残差是否为白噪声。% 对Burg算法产生的残差进行Ljung-Box检验 [h_burg, pValue_burg] lbqtest(residual_burg, Lags, [5, 10, 20], Alpha, 0.05); fprintf(Ljung-Box检验结果 (Burg模型残差):\n); for i 1:length(h_burg) if h_burg(i) 0 fprintf( 滞后%d: 无法拒绝原假设 (p%.4f)残差可能为白噪声。\n, [5,10,20](i), pValue_burg(i)); else fprintf( 滞后%d: 拒绝原假设 (p%.4f)残差不是白噪声模型可能拟合不足。\n, [5,10,20](i), pValue_burg(i)); end end如果p值大于显著性水平如0.05则不能拒绝“残差是白噪声”的原假设说明模型拟合得较好已经提取了信号中主要的线性依赖关系。5. 阶数选择关键的模型复杂度权衡之前我们假设了阶数p4。但对于一个未知信号如何选择p这是一个偏差-方差权衡问题阶数太低模型太简单无法捕捉信号全部特征欠拟合阶数太高模型会开始拟合噪声导致过拟合泛化能力变差。5.1 信息准则法AIC与BIC最常用的方法是计算不同阶数p下的AIC或BIC值选择使该值最小的p。max_order 30; % 搜索最大阶数 AIC zeros(max_order, 1); BIC zeros(max_order, 1); for p_test 1:max_order [a, sigma2_est] arburg(x, p_test); N length(x); % AIC N * ln(σ²) 2 * p AIC(p_test) N * log(sigma2_est) 2 * p_test; % BIC N * ln(σ²) p * ln(N) BIC(p_test) N * log(sigma2_est) p_test * log(N); end % 找到最小AIC/BIC对应的阶数 [~, p_aic] min(AIC); [~, p_bic] min(BIC); figure(Position, [100, 100, 900, 400]); subplot(1,2,1); plot(1:max_order, AIC, o-, LineWidth, 1.5); grid on; xlabel(模型阶数 p); ylabel(AIC值); title(AIC准则定阶); hold on; plot(p_aic, AIC(p_aic), r*, MarkerSize, 15); hold off; legend(AIC, sprintf(最优阶数 p%d, p_aic)); subplot(1,2,2); plot(1:max_order, BIC, s-, LineWidth, 1.5); grid on; xlabel(模型阶数 p); ylabel(BIC值); title(BIC准则定阶); hold on; plot(p_bic, BIC(p_bic), r*, MarkerSize, 15); hold off; legend(BIC, sprintf(最优阶数 p%d, p_bic));BIC准则相比AIC对模型复杂度的惩罚更重因此BIC选择的阶数通常不会高于AIC选择的阶数有时会更低。对于我们的测试信号真实阶数4AIC和BIC曲线应该在p4附近出现一个明显的“肘点”或最小值。如果曲线一直下降或没有明显转折可能意味着信号不适合用低阶AR模型描述或者含有较强噪声。5.2 观察偏相关函数PACF对于AR模型其偏相关函数在滞后超过模型真实阶数p后理论上应截尾接近0。因此观察样本偏相关函数也是定阶的辅助手段。figure; parcorr(x, 50); % 计算并绘制滞后50以内的样本偏自相关函数并给出置信区间 title(样本偏自相关函数 (PACF));在生成的PACF图中你会看到蓝色的条形图代表不同滞后下的偏相关系数值红色的虚线是置信区间。对于一个AR(p)过程前p个偏相关系数通常显著不为零超出红色虚线范围而p阶之后的偏相关系数则应落在置信区间内呈现“截尾”现象。你可以结合此图和信息准则共同确定阶数。6. 回到数据压缩AR模型如何实现压缩最后我们回到作业标题中的“数据压缩”。AR模型用于压缩的核心在于其参数化表示。假设我们有一段语音信号采样率8kHz一帧20ms即160个样本点。原始存储需要存储160个浮点数假设每个4字节共640字节。AR模型压缩我们用一个AR(10)模型去拟合这一帧信号。通过arburg函数我们得到10个AR系数a1到a10和1个增益噪声方差σ²的平方根或叫激励增益。我们只需要存储这11个参数。如果每个参数也用4字节浮点数存储只需44字节。此外我们还需要存储残差信号原始信号与AR模型预测信号的差但残差信号的熵通常比原始信号低可以用更少的比特进行量化编码如标量量化或矢量量化。在实际的语音编码器如线性预测编码LPC中正是利用了这一原理。发送端分析语音帧得到AR参数通常转换为更稳定的线谱对LSP或反射系数和残差激励量化后传输接收端用AR模型合成滤波器和激励信号重建出语音。压缩比就来自于用少量参数替代了大量的时域样本。在MATLAB里你可以做一个简单的模拟演示% 模拟一帧信号 frame x(1:160); % 取前160点作为一帧 p_comp 10; % 使用10阶AR模型 [a_comp, var_comp] arburg(frame, p_comp); % 使用估计的模型合成预测该帧信号 % 注意这里为了演示我们使用“零激励”来合成这相当于只用了模型的“骨架” % 实际压缩中需要用编码后的残差作为激励 frame_synth filter(1, [1, a_comp], sqrt(var_comp) * randn(160, 1)); % 用白噪声激励 % 计算原始信号和合成信号的误差 mse mean((frame - frame_synth).^2); fprintf(使用AR(%d)模型合成信号与原帧信号的均方误差(MSE)为: %.6f\n, p_comp, mse); % 绘制对比 figure; subplot(2,1,1); plot(frame); hold on; plot(frame_synth, r--); hold off; legend(原始信号, AR模型合成信号); title(原始信号 vs. AR模型合成信号); xlabel(样本点); ylabel(幅值); grid on; subplot(2,1,2); plot(frame - frame_synth); title(合成误差信号); xlabel(样本点); ylabel(幅值); grid on;这个演示展示了用AR模型参数10个系数1个增益和一段随机白噪声就能大致恢复出原始信号的“轮廓”和频谱特性。误差主要来自于我们用了完全不同的随机激励。在真正的编码器中会对残差进行精细的量化编码使得用更少的数据量实现高质量的重建。完成这次作业你收获的不仅仅是一段能跑通的MATLAB代码更是一套完整的信号建模、参数估计、模型评估和应用的思维框架。下次当你听到“LPC”、“线性预测”这些词时你就能立刻联想到今天做的这一切它就是用一组AR系数去捕捉信号最本质的相关结构。 AR模型原理与MATLAB实现：从参数估计到数据压缩应用
1. 作业背景与核心任务拆解这周的作业题目是“随机信号AR模型及MATLAB实现”看起来是数据压缩课程里的一次实践。很多同学拿到这种题目第一反应可能就是去网上找段代码改改参数把图跑出来交差。但说实话这样学完你可能只知道“AR模型在MATLAB里用aryule或者arburg函数”至于它为什么能用来“压缩”数据它的参数到底代表了什么估计还是一头雾水。我自己当年学信号处理的时候也这么干过后来在实际项目中用到AR模型做预测和特征提取才发现基础没打牢的亏有多大。所以咱们这次不光是完成作业更要把AR模型从原理到实操再到它和数据压缩的关联彻底捋清楚。AR模型全称自回归模型它的核心思想特别直观当前时刻的信号值可以用过去若干个时刻的信号值的线性组合再加上一个随机扰动白噪声来表示。这个“过去若干个时刻”就是模型的阶数。你想想看如果我能用一个简单的线性公式几个系数和一点随机噪声就很好地描述或预测一个信号那我是不是就不用存储原始信号的全部数据点了我只需要存储这几个系数和残差预测误差这不就是一种“压缩”或“简洁表示”吗这就是AR模型在数据压缩尤其是语音、音频编码中的理论基础。本次作业的核心我认为可以拆解为三个层次第一理解AR模型的数学表达和参数意义第二掌握在MATLAB中估计AR模型参数的经典方法比如Yule-Walker方程法、Burg算法第三也是最关键的如何评估你建立的AR模型好不好以及如何直观地展示这个过程。很多教程只教到第二步但第三步才是体现你理解深度的部分。2. AR模型的数学原理不只是几个公式我们先从最根本的数学表达式开始。对于一个p阶的AR模型记作AR(p)它的定义式是[ x[n] -\sum_{k1}^{p} a_k x[n-k] w[n] ]这里x[n]是我们在n时刻观测到的信号值。a_1, a_2, ..., a_p就是我们要求解的AR模型系数也叫反射系数或预测系数。注意公式前面的负号这是信号处理领域的习惯写法MATLAB的相关函数也遵循这个约定。w[n]是驱动白噪声通常假设它均值为0方差为σ²。这个噪声代表了模型无法用过去信号线性解释的部分是信号的“创新”或“随机”成分。这个公式揭示了AR模型的本质它试图用信号自身的过去来预测现在。系数a_k的大小和正负决定了过去第k个时刻的值对当前值的影响程度和方向。举个例子在语音信号中AR模型能非常有效地模拟声道的共振特性那些系数就对应着声道形状的参数。那么如何从一段观测信号x[1], x[2], ..., x[N]中估计出这p个系数a_k和噪声方差 σ² 呢这就引出了几种经典的估计算法它们在MATLAB中都有对应的函数实现。2.1 Yule-Walker方程法这是最经典的方法基于信号的自相关函数。思路是将模型定义式两边同时乘以x[n-m](m0,1,...,p)然后求数学期望。经过一系列推导这里不展开复杂的公式我们可以得到著名的Yule-Walker方程[ R_{xx}[m] -\sum_{k1}^{p} a_k R_{xx}[m-k], \quad m 1,2,...,p ]以及噪声方差[ \sigma^2 R_{xx}[0] \sum_{k1}^{p} a_k R_{xx}[k] ]其中R_xx[m]是信号x[n]的自相关函数在滞后m处的估计值。你看原来求解系数a_k的问题转化为了求解一个以自相关函数值为系数的线性方程组。这个方程组是Toeplitz结构的可以用高效的Levinson-Durbin递归算法来求解。在MATLAB中aryule函数就是基于Yule-Walker方程来估计AR参数的。它的优点是计算稳定理论完备。但缺点是对自相关函数的估计误差比较敏感特别是当数据长度较短时估计的偏差可能会较大。2.2 Burg算法最大熵谱估计Burg算法是一种更现代、通常表现更好的方法。它不直接估计自相关函数而是以前后向预测误差功率最小为准则递推地求解反射系数它和AR系数有确定转换关系。Burg算法能保证产生的AR模型是稳定的即所有极点都在单位圆内并且对于短数据序列其谱估计的分辨率和精度往往优于Yule-Walker法。在MATLAB中arburg函数实现了Burg算法。对于大多数实际应用尤其是本次作业我强烈推荐优先使用arburg。它的估计通常更准确特别是当你处理的信号长度有限时。理解这两种方法的区别很重要它不是简单的“哪个函数更好”而是背后优化准则和计算路径的不同。Yule-Walker法最小化前向预测误差但依赖于自相关函数的估计Burg法则同时最小化前向和后向预测误差直接操作数据。在作业中你可以尝试用两种方法对同一信号建模对比它们的结果这会是一个很好的加分点。3. MATLAB实战从信号生成到模型拟合光说不练假把式我们现在就用MATLAB来走一遍完整的流程。假设我们要分析一个由已知AR模型生成的测试信号这样我们可以将估计出的参数与真实值对比验证方法的有效性。3.1 生成一个测试用的AR信号首先我们得有一个信号。我们可以自己定义一个AR(4)模型然后用它来生成一段数据。这样真实的系数a_k是已知的便于后续对比。% 1. 定义真实的AR(4)模型参数 true_a [1.0, -0.9, 0.5, -0.3]; % 注意这里我们按MATLAB习惯给出的是A(z)1 a1*z^{-1}...的系数 % 即 A [1, a1, a2, a3, a4] true_A [1, true_a]; % 完整的AR多项式系数向量 % 2. 生成驱动白噪声 N 1000; % 信号长度 sigma2 0.1; % 噪声方差 w sqrt(sigma2) * randn(N, 1); % 生成N点的高斯白噪声 % 3. 使用filter函数生成AR信号 % filter(B, A, X) 是用传递函数H(z)B(z)/A(z)对X滤波。 % 对于AR模型B(z)1 A(z)就是我们的AR多项式。 x filter(1, true_A, w); % 生成的信号x现在我们手头有了一段长度为1000点的信号x它是由系数为true_A、噪声方差为0.1的AR(4)模型产生的。我们的任务就是假装不知道true_A仅通过观测信号x来估计出AR模型的阶数p和系数a_k。3.2 使用Burg算法估计AR参数我们先使用推荐的arburg函数。这里有一个关键步骤如何确定阶数p对于测试信号我们知道是4阶。但对于未知信号我们需要定阶。常用方法有AIC赤池信息准则或BIC贝叶斯信息准则它们会在模型拟合优度和复杂度之间做一个权衡。为了简化本次我们先假设通过观察信号的偏相关函数PACF或经验确定了p4。% 使用arburg函数估计AR(4)模型参数 p 4; % 假设我们已经知道或确定了阶数为4 [est_a_burg, sigma2_burg] arburg(x, p); % 注意arburg返回的est_a_burg是[A(2), A(3), ..., A(p1)]即a1, a2, ..., ap % 为了与我们的true_a对比需要构建完整的A向量 est_A_burg [1, est_a_burg];3.3 使用Yule-Walker方法估计AR参数我们也用aryule函数做一次用于对比。% 使用aryule函数估计AR(4)模型参数 [est_a_yule, sigma2_yule] aryule(x, p); est_A_yule [1, est_a_yule];3.4 结果对比与分析现在让我们把真实值和两种方法的估计值放在一起比较。fprintf(真实AR系数: [1, %.4f, %.4f, %.4f, %.4f]\n, true_a); fprintf(Burg估计系数: [1, %.4f, %.4f, %.4f, %.4f]\n, est_a_burg); fprintf(Yule-Walker估计系数: [1, %.4f, %.4f, %.4f, %.4f]\n, est_a_yule); fprintf(\n); fprintf(真实噪声方差: %.4f\n, sigma2); fprintf(Burg估计噪声方差: %.4f\n, sigma2_burg); fprintf(Yule-Walker估计噪声方差: %.4f\n, sigma2_yule);运行这段代码你可能会看到Burg算法估计的系数通常更接近真实值噪声方差的估计也更准。这就是为什么在实践里Burg算法更受青睐。但Yule-Walker法并非一无是处它在理论分析上更直观且对于某些特定类型的信号或足够长的数据表现也不错。4. 模型评估与可视化让结果“说话”把参数估计出来只是第一步。一个负责任的建模过程必须包含模型评估。我们需要回答这个AR(4)模型拟合得好不好有没有可能用更低或更高的阶数4.1 绘制功率谱密度PSD进行对比AR模型的一个重要应用就是谱估计。我们可以用估计出的AR参数来计算信号的功率谱密度并与经典的非参数化方法如周期图法进行对比。% 计算频率向量 NFFT 1024; Fs 1; % 假设采样频率为1Hz频率轴归一化到(0, 0.5) f (0:NFFT/2)*(Fs/NFFT); % 正频率部分 % 1. 使用Burg算法估计的AR参数计算PSD [H_burg, freq_burg] freqz(1, est_A_burg, NFFT/21, Fs); % 计算AR模型的频率响应 PSD_burg sigma2_burg * abs(H_burg).^2; % AR模型的PSD公式σ² / |A(e^{jω})|² % 2. 使用Yule-Walker算法估计的AR参数计算PSD [H_yule, freq_yule] freqz(1, est_A_yule, NFFT/21, Fs); PSD_yule sigma2_yule * abs(H_yule).^2; % 3. 使用周期图法非参数化方法作为参考 [Pxx_periodogram, f_periodogram] periodogram(x, hamming(length(x)), NFFT, Fs); % 绘制在同一张图上进行对比 figure(Position, [100, 100, 900, 600]); subplot(2,1,1); plot(f, 10*log10(PSD_burg), b-, LineWidth, 1.5); hold on; plot(f, 10*log10(PSD_yule), r--, LineWidth, 1.5); plot(f_periodogram, 10*log10(Pxx_periodogram), k:, LineWidth, 0.8); hold off; grid on; xlabel(归一化频率); ylabel(功率谱密度 (dB)); title(AR模型谱估计 vs. 周期图法); legend(Burg算法 (AR), Yule-Walker算法 (AR), 周期图法, Location, Best); xlim([0, 0.5]); % 绘制残差分析图 subplot(2,1,2); % 使用Burg模型对原信号进行预测滤波 x_hat_burg filter([0 -est_a_burg], 1, x); % 注意系数符号这是根据AR定义式做的预测 residual_burg x - x_hat_burg; % 残差 真实值 - 预测值 plot(residual_burg); grid on; xlabel(样本点); ylabel(幅值); title(Burg算法模型残差序列);通过这张图你可以直观地看到谱估计对比AR模型谱特别是Burg算法通常比周期图法平滑得多分辨率也更高。周期图法起伏剧烈方差大而AR谱能清晰地显示出信号的谱峰这对应于AR模型的极点位置。如果两种AR方法估计的谱形状差异很大说明模型估计可能不稳定或者阶数选择不当。残差分析一个拟合良好的AR模型其残差序列应该近似为白噪声即没有明显的自相关结构。你可以通过观察残差序列的波形如上图初步判断更严谨的做法是计算残差的自相关函数ACF并检验其是否在零滞后外接近零。4.2 定量评估残差白噪声检验我们可以用Ljung-Box检验来定量判断残差是否为白噪声。% 对Burg算法产生的残差进行Ljung-Box检验 [h_burg, pValue_burg] lbqtest(residual_burg, Lags, [5, 10, 20], Alpha, 0.05); fprintf(Ljung-Box检验结果 (Burg模型残差):\n); for i 1:length(h_burg) if h_burg(i) 0 fprintf( 滞后%d: 无法拒绝原假设 (p%.4f)残差可能为白噪声。\n, [5,10,20](i), pValue_burg(i)); else fprintf( 滞后%d: 拒绝原假设 (p%.4f)残差不是白噪声模型可能拟合不足。\n, [5,10,20](i), pValue_burg(i)); end end如果p值大于显著性水平如0.05则不能拒绝“残差是白噪声”的原假设说明模型拟合得较好已经提取了信号中主要的线性依赖关系。5. 阶数选择关键的模型复杂度权衡之前我们假设了阶数p4。但对于一个未知信号如何选择p这是一个偏差-方差权衡问题阶数太低模型太简单无法捕捉信号全部特征欠拟合阶数太高模型会开始拟合噪声导致过拟合泛化能力变差。5.1 信息准则法AIC与BIC最常用的方法是计算不同阶数p下的AIC或BIC值选择使该值最小的p。max_order 30; % 搜索最大阶数 AIC zeros(max_order, 1); BIC zeros(max_order, 1); for p_test 1:max_order [a, sigma2_est] arburg(x, p_test); N length(x); % AIC N * ln(σ²) 2 * p AIC(p_test) N * log(sigma2_est) 2 * p_test; % BIC N * ln(σ²) p * ln(N) BIC(p_test) N * log(sigma2_est) p_test * log(N); end % 找到最小AIC/BIC对应的阶数 [~, p_aic] min(AIC); [~, p_bic] min(BIC); figure(Position, [100, 100, 900, 400]); subplot(1,2,1); plot(1:max_order, AIC, o-, LineWidth, 1.5); grid on; xlabel(模型阶数 p); ylabel(AIC值); title(AIC准则定阶); hold on; plot(p_aic, AIC(p_aic), r*, MarkerSize, 15); hold off; legend(AIC, sprintf(最优阶数 p%d, p_aic)); subplot(1,2,2); plot(1:max_order, BIC, s-, LineWidth, 1.5); grid on; xlabel(模型阶数 p); ylabel(BIC值); title(BIC准则定阶); hold on; plot(p_bic, BIC(p_bic), r*, MarkerSize, 15); hold off; legend(BIC, sprintf(最优阶数 p%d, p_bic));BIC准则相比AIC对模型复杂度的惩罚更重因此BIC选择的阶数通常不会高于AIC选择的阶数有时会更低。对于我们的测试信号真实阶数4AIC和BIC曲线应该在p4附近出现一个明显的“肘点”或最小值。如果曲线一直下降或没有明显转折可能意味着信号不适合用低阶AR模型描述或者含有较强噪声。5.2 观察偏相关函数PACF对于AR模型其偏相关函数在滞后超过模型真实阶数p后理论上应截尾接近0。因此观察样本偏相关函数也是定阶的辅助手段。figure; parcorr(x, 50); % 计算并绘制滞后50以内的样本偏自相关函数并给出置信区间 title(样本偏自相关函数 (PACF));在生成的PACF图中你会看到蓝色的条形图代表不同滞后下的偏相关系数值红色的虚线是置信区间。对于一个AR(p)过程前p个偏相关系数通常显著不为零超出红色虚线范围而p阶之后的偏相关系数则应落在置信区间内呈现“截尾”现象。你可以结合此图和信息准则共同确定阶数。6. 回到数据压缩AR模型如何实现压缩最后我们回到作业标题中的“数据压缩”。AR模型用于压缩的核心在于其参数化表示。假设我们有一段语音信号采样率8kHz一帧20ms即160个样本点。原始存储需要存储160个浮点数假设每个4字节共640字节。AR模型压缩我们用一个AR(10)模型去拟合这一帧信号。通过arburg函数我们得到10个AR系数a1到a10和1个增益噪声方差σ²的平方根或叫激励增益。我们只需要存储这11个参数。如果每个参数也用4字节浮点数存储只需44字节。此外我们还需要存储残差信号原始信号与AR模型预测信号的差但残差信号的熵通常比原始信号低可以用更少的比特进行量化编码如标量量化或矢量量化。在实际的语音编码器如线性预测编码LPC中正是利用了这一原理。发送端分析语音帧得到AR参数通常转换为更稳定的线谱对LSP或反射系数和残差激励量化后传输接收端用AR模型合成滤波器和激励信号重建出语音。压缩比就来自于用少量参数替代了大量的时域样本。在MATLAB里你可以做一个简单的模拟演示% 模拟一帧信号 frame x(1:160); % 取前160点作为一帧 p_comp 10; % 使用10阶AR模型 [a_comp, var_comp] arburg(frame, p_comp); % 使用估计的模型合成预测该帧信号 % 注意这里为了演示我们使用“零激励”来合成这相当于只用了模型的“骨架” % 实际压缩中需要用编码后的残差作为激励 frame_synth filter(1, [1, a_comp], sqrt(var_comp) * randn(160, 1)); % 用白噪声激励 % 计算原始信号和合成信号的误差 mse mean((frame - frame_synth).^2); fprintf(使用AR(%d)模型合成信号与原帧信号的均方误差(MSE)为: %.6f\n, p_comp, mse); % 绘制对比 figure; subplot(2,1,1); plot(frame); hold on; plot(frame_synth, r--); hold off; legend(原始信号, AR模型合成信号); title(原始信号 vs. AR模型合成信号); xlabel(样本点); ylabel(幅值); grid on; subplot(2,1,2); plot(frame - frame_synth); title(合成误差信号); xlabel(样本点); ylabel(幅值); grid on;这个演示展示了用AR模型参数10个系数1个增益和一段随机白噪声就能大致恢复出原始信号的“轮廓”和频谱特性。误差主要来自于我们用了完全不同的随机激励。在真正的编码器中会对残差进行精细的量化编码使得用更少的数据量实现高质量的重建。完成这次作业你收获的不仅仅是一段能跑通的MATLAB代码更是一套完整的信号建模、参数估计、模型评估和应用的思维框架。下次当你听到“LPC”、“线性预测”这些词时你就能立刻联想到今天做的这一切它就是用一组AR系数去捕捉信号最本质的相关结构。