burg4.m

雷达信号带限外推。包括burg外推的方法

clear;
N=16;                             %数据长度 改变数据长度会导致分辨率的变化；
n=2:N+1;
k=2:N+1;
F1=0.320;
F2=0.560;
sga=0.01;
r1=randn(1,N);
r2=randn(1,N);
Z1=sga*sqrt(-2*log(r1)).*(cos(2*pi*r2));
Z2=sga*sqrt(-2*log(r1)).*(sin(2*pi*r2));
Za=Z1+j*Z2;
X=3*exp(j*F1*pi*n)+2*exp(j*F2*pi*n)+Za;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Burg算法 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S=0;
for n=1:N+1
    e(1,n)=X(n);
    b(1,n)=X(n);
    S=S+(X(n)*X(n));
end
U2(1)=(1/(N+1))*S;

W=0;S=0;
for n=1:N
    W=W+(e(1,n+1)*b(1,n));
    S=S+(e(1,n+1)*e(1,n+1)+b(1,n)*b(1,n));
end
K(1)=-2*W/S;
a(1,1)=K(1);
U2(2)=(1-K(1)^2)*U2(1);
for n=1:N
    e(2,n+1)=e(1,n+1)+K(1)*b(1,n);
    b(2,n+1)=b(1,n)+K(1)*e(1,n+1);
end

W=0;S=0;
for p=2:M
    for m=p:N
        W=W+(e(p,m+1)*b(p,m))
        S=S+(e(p,m+1)*e(p,m+1)+b(p,m)*b(p,m))
    end
    K(p)=-2*W/S;
    W=0;S=0;
    for m=p:N
        e(p+1,m+1)=e(p,m+1)+K(p)*b(p,m);
        b(p+1,m+1)=b(p,m)+K(p)*e(p,m+1);
    end
    a(p,p)=K(p);
    for n=1:p-1
        a(p,n)=a(p-1,n)+a(p,p)*a(p-1,p-n);
    end
    U2(p+1)=(1-K(p)^2)*U2(p);
end
disp('输出模型参数a');
for k=1:M
    disp(a(M,k));
end
disp('输出误差功率');
disp(U2(M+1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% AR模型参数确定后计算出功率谱 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Z=0;
w=0:0.01:pi;                     %功率谱以2*pi为周期，又信号为实信号，只需输出0到PI即可；
for k=1:M 
    Z=Z+(a(M,k)*exp(-j*k*w));
end
PXX=U2(M)./((abs(1+Z)).^2);    %得到功率谱函数；
%F=w*FS/(pi*2);                   %将角频率坐标换算成HZ坐标，便于观察；
plot(w,PXX)
grid