随机数字信号处理实验报告三——Levinson和Burg递推法MATLAB实现

完整的实验报告下载连接https://download.csdn.net/download/LIsaWinLee/14884452

一、实验原理

    随机信号的功率谱密度用来描述信号的能量特征随频率的变化关系。功率谱密度简称为功率谱，是自相关函数的傅里叶变换。对功率谱密度的估计又称功率谱估计。

1.1 Levinson递推法

自相关法——列文森递推法是已知信号观测数据，估计功率谱。它的出发点是选择AR模型参数使预测误差功率最小。

假设信号xn 的数据区在0≤n≤N-1 范围，有P个预测系数，N个数据经过冲激响应为api(i=0,1…,p) 的滤波器，输出预测误差e(n) 的长度为N+P ,因此有预测误差功率为 。

1.2 Burg递推法

Levinson递推法需要由观测数据估计自相关函数，这是它的缺点。而Burg递推法则由信号观测数据直接计算AR模型参数。

Burg递推法思想：借助格型预测误差平均功率，选择kp 使其最小，求出kp 。之后再利用Levinson递推法求模型参数和输入噪声方差。

二、实验内容

信号为两个正弦信号加高斯白噪声，各正弦信号的信噪比均为10dB，长度为N，信号频率分别为 和 ，初始相位 ，取 ， 取不同数值：0.3，0.25。 为采样频率。分别用Levinson递推法和Burg法进行功率谱估计，并分析改变数据长度、模型阶数对谱估计结果的影响。

三、实验过程及结果

3.1 原始信号

原始信号为s=sin2πf1kfs+sin2πf2kfs ，观测信号是x=s+w 。取fs=1，f1=0.2，f2=0.3， 数据长度N为100，阶数p1为20。输出原始图像如图3.1所示。

图3.1 原始信号图像

3.2 Levinson递推法

3.2.1 研究数据长度对谱估计的影响

取fs=1，f1=0.2，f2=0.3， 阶数p1为20保持不变。改变不同的数据长度，比较数据长度对谱估计的影响。

1）数据长度N为35

图3.2 阶数20、数据长度35的功率谱估计

2）数据长度N为100

图3.3 阶数20、数据长度100的功率谱估计

3）数据长度N为1000

图3.4 阶数20、数据长度1000的功率谱估计

分析：由上述3个实验对比可以得出，在阶数一定时，观测数据较短时，估计误差较大，会出现谱峰频率偏移和谱线分裂；数据长度越长，估计自相关函数越准确，但计算量较大。

3.2.2 研究阶数对谱估计的影响

取fs=1，f1=0.2，f2=0.3， 数据长度为100保持不变。改变不同的阶数，比较阶数对谱估计的影响。

1）阶数p1为2

图3.5 阶数2、数据长度100的功率谱估计

2）阶数p1为8

图3.6 阶数8、数据长度100的功率谱估计

3）阶数p1为20

图3.7 阶数20、数据长度100的功率谱估计

分析：由上述3个实验对比可以得出，在数据长度一定时，当阶数较低，会使谱估计发生偏移，降低分辨率；阶数越高，分辨率越高；但阶数太高，会使估计误差加大，谱峰分裂。

3.3 Burg递推法

3.3.1 研究数据长度对谱估计的影响

取fs=1，f1=0.2，f2=0.3， 阶数M为20保持不变。改变不同的数据长度，比较数据长度对谱估计的影响。

1）数据长度N为35

图3.8 阶数20、数据长度35的功率谱估计

2）数据长度N为100

图3.9 阶数20、数据长度100的功率谱估计

3）数据长度N为1000

图3.10 阶数20、数据长度1000的功率谱估计

分析：由上述3个实验对比可以得出，在阶数一定时，观测数据较短时，估计误差较大，会出现谱峰频率偏移和谱线分裂；数据长度越长，估计自相关函数越准确，但计算量较大。频率越靠近的谱估计，需要的阶数越高。

3.3.2 研究阶数对谱估计的影响

取fs=1，f1=0.2，f2=0.3， 数据长度为900保持不变。改变不同的阶数，比较阶数对谱估计的影响。

1）阶数p1为4

图3.11 阶数4、数据长度900的功率谱估计

2）阶数p1为8

图3.12 阶数8、数据长度900的功率谱估计

3）阶数p1为20

图3.13 阶数20、数据长度900的功率谱估计

分析：由上述3个实验对比可以得出，在数据长度一定时，当阶数较低，会使谱估计发生偏移，降低分辨率；阶数越高，分辨率越高；但阶数太高，会使估计误差加大，谱峰分裂。

四、总结与展望

本次实验分别采用Levinson和Burg递推法进行功率谱估计，采用单一变量法研究。在数据长度和模型阶数中，一者不变时，不断改变另一个的大小，研究数据长度和模型阶数对谱估计结果的影响。通过实验，学习Levinson和Burg递推法的基本原理和一般流程，学会运用MATLAB编写程序，完成实验。通过实验对比得出，阶数越高，分辨率越高；数据长度越长，自相关谱估计越准确。不过不能一味的数值大，应该通过多次实验，设置最为合理的参数结果。

五、附录

5.1 Levinson递推法的程序

clc;

clear all;

fs=1;%采样频率

Ts=1/fs;

N=100;%数据长度

p1=8;%阶数

f1=0.2*fs;f2=0.3*fs;%信号频率

pha1=0;pha2=0;%初始相位

SNR=2;%设置信噪比

%产生信号

w=randn(1,N);

Am=sqrt(2*10^(SNR/10));

k=1:N;

x1=Am*sin(2*pi*f1*k*Ts+pha1);x2=Am*sin(2*pi*f2*k*Ts+pha2);

xx=x1+x2;x=w+x1+x2;

figure(1)

subplot(2,1,1);plot(xx);title('两个正弦信号波形');grid on;

subplot(2,1,2);plot(x);title('加噪信号波形');grid on;

%计算自相关函数

R=[];

for m=1:N

    s=0;

    for n=1:N-(m-1)

        s=s+x(n)*x(n+m-1);

    end

    R(m)=s/N;

end

%列文森递推法

a(1,1)=-R(2)/R(1);cov(1)=R(1)+a(1,1)*R(2);

for p=2:p1

    sum=0;

    for k=1:p-1

        sum=sum+a(p-1,k)*R(p-k+1);

    end

    a(p,p)=-(R(p+1)+sum)/cov(p-1);

    for k=1:p-1

        a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);

    end

    cov(p)=(1-(abs(a(k,k)))^2)*cov(p-1);

end

%计算功率谱，功率谱以2*pi为周期，又信号为实信号，只需要输出0到p1即可

W=0.01:0.01:pi;Z=0;

for k=1:p1

    Z=Z+(a(p1,k).*exp(-j*k*W));

end

PSD=1./((abs(1+Z)).^2);%算出功率谱

F=W*fs/(pi*2);%角频率坐标换算成频率坐标

figure(2)

plot(F,abs(PSD));xlabel('频率（Hz）');ylabel('功率谱');

title('自相关法-列文森Levenson递推法的功率谱估计');grid;

figure(3)

p=1:p1;plot(p,cov(p),'b');xlabel('模型阶数');ylabel('预测误差功率大小');

title('预测误差功率变化趋势');grid;

5.1 Burg递推法的程序

clear all;

clear

%产生信号

fs=1;%设采样频率为1

N=900;%数据长度 改变数据长度会导致分辨率的变化

f1=0.2*fs;%第一个正弦信号的频率，f1/fs=0.2

f2=0.3*fs;%第二个正弦信号的频率，f1/fs=0.3

M=20;%滤波器阶数的最大取值，超过则认为歹加太大而放弃

n=1:N;

s=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs);%s为原始信号

x=awgn(s,10);%x为观测信号，即对原始信号加入白噪声，信噪比10dB

for i=1:N

    ef(1,i)=x(i);

    eb(1,i)=x(i);

end

sum=0;

for i=1:N

    sum=sum+x(i)*x(i);

end

r(1)=sum/N;

%Burg递推

for p=2:M

    %求解第p个反射系数

    sum1=0;

    for n=p:N

        sum1=sum1+ef(p-1,n)*eb(p-1,n-1);

    end

    sum1=-2*sum1;

    sum2=0;

    for n=p:N

        sum2=sum2+ef(p-1,n)*ef(p-1,n)+eb(p-1,n-1)*eb(p-1,n-1);

    end

    k(p-1)=sum1/sum2;

    %求解预测误差平均功率

    r(p)=(1-k(p-1)*k(p-1))*r(p-1);

    %求解p阶白噪声方差

    q(p)=r(p);

    %系数a

    if p>2

        for i=1:p-2

            a(p-1,i)=a(p-2,i)+k(p-1)*a(p-2,p-1-i);

        end

    end

    a(p-1,p-1)=k(p-1);

    %求解前向预测误差

    for n=p+1:N

        ef(p,n)=ef(p-1,n)+k(p-1)*eb(p-1,n-1);

    end

    %求解后向预测误差

    for n=p:N-1

        eb(p,n)=eb(p-1,n-1)+k(p-1)*ef(p-1,n);

    end

end

%计算功率谱

for j=1:N

    sum3=0;

    sum4=0;

    for i=1:p-1

        sum3=sum3+a(p-1,i)*cos(2*pi*i*j/N);

    end

    sum3=1+sum3;

    for i=1:p-1

        sum4=sum4+a(p-1,i)*sin(2*pi*i*j/N);

    end

    pxx=sqrt(sum3*sum3+sum4*sum4);

    pxx=q(M)/pxx;

    pxx=10*log10(pxx);

    pp(j)=pxx;

end

%画出功率谱

ff=1:N;

ff=ff/N;

figure;

plot(ff,pp),axis([0 0.5 -20 10]),xlabel('频率'),ylabel('幅度（dB）'),title('自相关法-Burg递推法的功率谱P');

grid;