matlab利用levinson-durbin算法进行功率谱估计仿真实验

Levinson-Durbinalgorithm利用Levinson-Durbin算法计算序列x(n)的功率谱，cos0.3cos0.32xnnnwn其中w(n)为高斯白噪声序列（1）假设x(n)的采样点数为128，AR模型阶数为80，（2）假设x(n)的采样点数为512，AR模型阶数为80。1、数学模型Levinson算法主要是从AR(k)模型参数出发，得到AR(k+1)阶模型参数。目前已经得到k阶Yule-Walker方程参数{𝑎𝑘,1,𝑎𝑘,2,⋯,𝑎𝑘,𝑘,𝜎𝑘2}，接下来以此求解k+1阶Yule-Walker方程，k阶Yule-Walker方程为：2,1,1011010100kkkkRRRkaRRRkaRkRkRK+1阶Yule-Walker方程为：211,11,1,k11(0)(1)(k)(k1)(1)(0)(k1)(k)0(k)(k1)(0)(1)0(k1)(k)(1)(0)0kkkkkRRRRaRRRRaRRRRaRRRR将k阶方程的系数矩阵增加一行一列，成为一下的扩大方程：2,1,(0)(1)(k)(k1)1(1)(0)(k1)(k)0(k)(k1)(0)(1)0(k1)(k)(1)(0)0kkkkkRRRRRRRRaRRRRaRRRRD其中𝐷𝑘由下式定义：,,001,1kkkikiDaRkia将以上扩大方程的行倒序，列也倒序，得到预备方程：,,12(0)(1)(k)(k1)00(1)(0)(k1)(k)0(k)(k1)(0)(1)(k1)(k)(1)(0)1kkkkkDRRRRRRRRaRRRRaRRRR将待求解的k+1阶Yule-Walker方程的解表示成扩大方程的解和预备方程的解的线性组合形式：1,1,1,11,.,11,111001kkkkkkkkkkkkaaaaaaa或：1,,1,1,1,2,,kikikkkiaaaik式中𝛾𝑘+1是待定系数，称为反射系数，上式中各项都右乘k+1阶系数矩阵，得到：221120000kkkkkkDD由此可以得到：12kkkD22221111kkkkkkD由扩大方程的第一个方程可以得到：2,10kkkiiRaRi从以上的推导可以归纳出如下的由k阶模型参数求k+1阶模型参数的计算公式：2,10kkkiiRaRi,,001,1kkkikiDaRkia12kkkD222111kkk1,,1,11,k11,1,2,,;kikikkkikkaaaika对于AR(p)模型，递推计算到k+1=p为止，得到AR(p)模型的参数后，根据以下式子即可得到功率谱估计值：22,1ˆ1pjpjipiiSeae2、算法模型Levinson-Durbin算法计算流程如下：（1）初始化：2000,00,1,a1RDR（2）开始迭代计算：先由12kkkD，根据k阶的𝐷𝑘与𝜎𝑘2求出反射系数𝛾𝑘+1，之后利用222111kkk以及已经得到的𝜎𝑘2求出𝜎𝑘+12，利用1,,1,1kikikkkiaaa，求出k+1阶的AR模型参数。利用得到的𝑎𝑘+1,𝑖根据,01kkkiiDaRki可以求出𝐷𝑘+1，自此完成从k阶到k+1阶的迭代。如此依次递推计算知道达到要求的阶数时停止计算即可完成对AR(p)模型参数的估计。3、源代码clear;clc;n=[1:128];%时间序列p=80;%AR模型阶数N=length(n);%数据点数x=cos(0.3*pi*n)+1*cos(0.32*pi*n);x=awgn(x,4);%生成带噪声的数据序列%-----------自相关序列-----------x1=[xzeros(1,N-1)];fori=1:length(n)r(i)=sum(x.*x1(i:N+i-1))/N;end%-------------------------------%--------迭代k=0---------------sigma=r(1);D=r(2);gamma=D/sigma;a=zeros(1,p);%--------迭代k=1---------------a(1)=1;a(2)=-gamma;sigma=r(1)+a(2)*r(2);D=r(3)+a(2)*r(2);%--------迭代k2---------------fork=2:p-1gamma=D/sigma;ar=[a(1)fliplr(a(2:k))zeros(1,p-k)];a(2:k)=a(2:k)-gamma*ar(2:k);%公式4.25a(k+1)=-gamma;sigma=(1-gamma*gamma)*sigma;D=sum(a(1:k+1).*fliplr(r(2:k+2)));end%----根据AR模型参数画出功率谱--------------l=[1:p-1];num=1000;f=2*[1:1:num]/num;i=sqrt(-1);forjj=1:numS(jj)=sigma/((abs(1+sum(a(2:p).*exp(-i*2*pi*l*jj/num)))).^2);End%S=20*log10(S);plot(f(1:num/2),S(1:num/2))4、仿真结果与分析在数据点数为128，阶数为80时得到的功率谱为：图1数据点数为128，阶数为80，频率间隔0.02时的功率谱数据点数为512，阶数为80时得到的功率谱为：图2数据点数为512，阶数为80，频率间隔0.02时的功率谱将输入两个余弦信号的频率间隔改为0.01后，在相同信噪比下得到的功率谱为：图3数据点数为128，阶数为80，频率间隔0.01时的功率谱图4数据点数为128，阶数为80，频率间隔0.01时的功率谱可以看到，数据点数不同的两种情况下利用Levinson-Durbin算法都可以在0.3和0.32出区分出两个谱峰，对比图1与图2，可以看到，数据点数越多，能量就越多的集中在0.3与0.32上。接下来考虑减小频率间隔，如图3,4所示，512点的情况可以将两个频率分的更开。说明数据点数越多，对频谱的分辨能力就越强。其原因在于采样点数越多，取样自相关函数就越接近原本的自相关函数，以此为基础进行的谱估计也更精确。