报告2AR过程的线性建模与功率谱估计.doc

实验报告2——AR过程的线性建模与功率谱估计一、 实验目的1•理解AR过程的产生机理，复习实验1估计自相关序列的方法2•利用估计出的自相关序列来求解信号的功率谱，即用周期图法来估计功率谱3. 分别采用自相关法(Yule-Walker法)，协方差法，Burg法，修正协方差法来估计功 率谱，并与周期图法进行比较，分析性能孰优孰劣4. 学习matlab在数字信号处理中的应用二、 实验过程和分析AR过程的线性建模与功率谱估计考虑AR过程：x(n) = a(1)x(n -1) + a (2) x( n - 2) + a (3) x(n - 3) + a (4) x(n - 4) + b(O)v(n)v(n)是单位方差白噪声a)取 b(0)=1, a(1)=2.7607, a(2)=-3.8106, a(3)=2.6535, a(4)=-0.9238，产生 x(n)的 N=256 个 样点用randn(1,N)产生单位方差高斯白噪声v(n)，用v(n)激励滤波器产生AR(4)过程， 即用x=filter(b,a,v)产生x(n)，b是滤波器分子系数，这里为b(0)=1，a是滤波器分母系数， a=[1 -2.7607 3.8106 -2.6535 0.9238]。

⑹ 计算其自相关序列的估计r (k)，并与真实的自相关序列值相比较x150010001 1 11 1 11 1 11 1 11 1 1: 1 1 11 1 11 1 11 1 1—來估#自相关——真实自相关1 1自相关估计5000-500-1000-150050100150200250300结论：真实自相关序列与估计出的序列对比如上图所示，两者在100点之前的形状相似, 幅度有一定差异，而且估计出的自相关序列有较大的波动，这是因为估计的点数较少， 使得估计精度不够，另外，估计自相关序列的下标越大，用来估计的点数就越少，因而 后面的估计值是很不精确的因为假设n>= 500处的数据为0,所有估计自相关不精确, 误差较大c)将下(k)的DTFT作为x(n)的功率谱估计，即:xP (ej®) = £ r (k)e-jk® =x x丄丨 X (ej® )|2NOO周期图法获取的功率谱与真实功率谱的比较0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9数字频率怎o o O5 4 3OOO-k =- N+1结论：周期图法求得的功率谱误差较大，波动也较大，两个谱峰的分辨也并不清晰，并 不精确，这是由于周期图法对自相关序列加了一个矩形窗。

d)利用所估计的自相关值和Yule-Walker法(自相关法)，估计a(1), a(2), a(3), a⑷和b(0)的 值，并讨论估计的精度a(l)=2.5018 a(2)= -3.2213 a(3)=2.0794 a(4)=-0.7025 b(0)=2.4198对比实际结果，可知Yule-Walker法(自相关法)误差较大，估计并不理想e)用(d)中所估计的a(k)和b(0)来估计功率谱为:P (e.) = 1 b(0)|2x1 + 才 a(k)e-jk®k=iCQP、S3<祥o利用Yule-Walkerte估计功率谱0.1 0.2 0..3 0.4 0.5 0..6 0.7 0.8 0.9 1数字频率"o o o5 4 3o20o JI-估计出的功率谱上图所示，难以看出有二个谱峰，这是由于自相关法在短数据量的情况 下估计是不精确的⑴将(c)和(e)的两种功率谱估计与实际的功率谱进行比较，画出它们的重叠波形OO功率谱比较0.1 0.2 0.3 0.4 0.5 0..6 0.7 0.8 0.9数字频率㈤o o O4 3 2OOOJ—-结论：用Yule-Walker法估计得到的谱曲线光滑，但无法识别出有二个谱峰，且b(0)的估 计明显是错的，使得幅度有较大的偏差；相关图的曲线不光滑，不好观察谱峰，但大致 走向与真实谱是一致的。

可以发现，Yule-Walker法(参数法)比周期图法(非参数法)精 度更高但这里Yule-Walker法得到的功率谱分辨率较低，这是由于对自相关序列加了一 个矩形窗引起的g)重复上面的(d)〜(f),只是估计AR参数分别米用如下方法：(1)协方差法；(2) Burg方 法；(3)修正协方差法试比较它们的功率谱估计精度协方差法a(1)=2.6843 a(2)=-3.6442 a(3)=2.5026 a(4)=-0.8668 b(0)=1.5506mp'En」芯徂⑷」aAWQO功率谱比较0.1 0.2 0..3 0.4 0.5 0..6 0.7 0.8 0.9数字频率"与Yule-Walker法估计出参数的误差相比，有非常大的改进，这是因为协方差法对自相关 序列估计时不需要对数据加窗，因此，对短数据记录，协方差法一般可获得比自相关法 更高的分辨率由上图可知，谱估计的幅度相差较大，但从形状和极点定位的角度而言, 估计结果很好的，曲线光滑，且极点的定位非常准确，短观测数据也能对极点精确定位 是协方差法的优点Burg方法a(l)=2.7872 a(2)=-3.8728 a(3)=2.7115 a(4)=-0.9467 b(0)=1.3006O功率谱比较0.1 0.2 0..3 0.4 0.5 0..6 0.7 0.8 0.9数字频率"Burg算法最小化前向加反向预测误差的平方和来求解全极点参数，为了保证稳定性，这 一最小化是相对于反射系数序贯地进行的，其估计精度低于协方差法。

此估计方法与 Yule-Walker法估计出参数的误差相比，有非常大的改进，但与协方差法比较，精度没有 其理想用估计出的参数来进行谱估计，如上图所示，蓝色虚线就是采用Burg法估计出 的功率谱，红线是真实功率谱，由图可知，估计结果较好，曲线光滑，且极点的定位较 为准确，Burg法是自相关法与协方差法的一个折中修正协方差法a(l)=2.7863 a(2)=-3.8716 a(3)=2.7106 a(4)=-0.9467 b(0)=1.3006功率谱比较0.1 0.2 0 3 0.4 0.5 0..6 0.7 0.8 0.9数字频率九修正协方差法是最小化前向加反射预测误差的平方和，峰值偏移对相位的敏感度较低， 比Yule-Walker法精度高，与Burg法的精度相当，但比协方差法的精度要低一些用估 计出的参数来进行谱估计，如上图所示，蓝色虚线就是采用修正协方差法估计出的功率 谱，红线是真实功率谱，由图可知，估计结果较好，曲线光滑，且极点的定位较为准确 结论：协方差法和修正协方差法可以分辨出本题中两个频率，说明它们相比于自相关法 有更高的分辨率这是因为它在形成的自相关估计是不需要对数据加窗Burg算法也没 有对数据加窗，因此其AR参数的估计也比自相关估计所获得的估计要精确。

附代码% x(n)= a(1)x(n-1) + a(2)x(n-2) + a(3)x(n-3) + a(4)x(n-4) + b(O)v(n) % v(n)是单位方差白噪声% Pb(0)=1, a=[2.7607,-3.8106,2.6535,-0.9238]clc;clear;close all;N = 256;M=1024;B=[1]；A=[1 -2.7607 3.8106 -2.6535 0.9238];a=[2.7607,-3.8106,2.6535,-0.9238];b0 = 1;nf = 0:1/1000:1;w = pi*nf;v = randn(1,N);x=filter(B,A,v);figure;st em(x);title( 'AR过程生成序列x(n)');%估计自相关Nr = N;r1 = zeros(1,Nr);for k = 0:(Nr-1)for n = 1:(N-k)r1(k+1) = r1(k+1) + x(n)*x(n+k);endr1(k+1) = r1(k+1)/N;endfigure;stem(r1, 'rx');hold on;%计算真是自相关r = ator(A,b0);for k = length(a)+2:Nr(k) = a(1)*r(k-1) + a(2)*r(k-2) + a(3)*r(k-3) + a(4)*r(k-4);endstem(r, 'b*');title('自相关估计');legend('估计自相关','真实自相关');grid on;%利用估计的自相关求功率谱P=zeros(1,M);P2=zeros(1,M);P_t=zeros(1,M);f = 0:1/(M/2-1):1;Y = fft(x,M);P=(Y.*conj(Y))/N;Py = 10*log10(P);figure;plot(f,Py(1:length(f)),'--');hold on;%真实功率谱for h =1:length(nf) temp = 0;for k=1:5temp = temp + A(k).*exp(-j*(k-1)*w(h));endP_a(h) = 1/(temp.*conj(temp));endP_a = 10*log10(abs(P_a));plot(nf,P_a,'r');title('周期法获取的功率谱和真实功率谱的比较'); legend('周期图功率谱'，'真实功率谱');xlabel('数字频率八pi');ylabel('功率谱/dB'); grid on;%利用所估计的自相关，估计AR参数a(1),a(2),a(3),a(4),b0 [est_a,epsilon] = rtoa(r1(1:5))est_bO = sqrt(epsilon)for h =1:length(nf)temp = 0;for k=1:5temp = temp + est_a(k).*exp(-j*(k-1)*w(h));endP_acor(h) = epsilon/(temp.*conj(temp));endP_acor = 10*log10(abs(P_acor)); figure;plot(nf,P_acor,'--');title('利用 Yule-Walker法估计功率谱');xlabel('数字频率 / \pi');ylabel('功率谱/dB'); grid on;figure;plot(nf,P_a,'r');hold on;plot(nf,P_acor,'--'); title('功率谱比较');legend('真实功率谱','Yule_ Walker法');xlabel('数字频率(\pi)');ylabel('功率谱(dB)') grid on;%利用协方差求功率谱p = 4;[a_cov,err_cov] = covm(x,p)for h = 1:length(nf) temp = 0;for k = 1:5temp = temp + a_cov(k) * exp(-j*(k-1)*w(h)); endP_cov。