burg3.m

高阶频谱分析防真实验,MATLAB实现.说明详细,容易看懂.

clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   产生信号    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

FS=100;%设采样频率为100，则根据F1/FS=0.2，F2/FS=0.3 OR 0.25 ，可以得到F1，F2便于计算；
M=16;
N=155;                             %数据长度 改变数据长度会导致分辨率的变化；
n=0:1/FS:N/FS;                    
F1=20;                            %第一个SIN信号的频率；
F2=20.5;                            %第二个SIN信号的频率，取25或者30；
X=sin(2*pi*F1*n)+sin(2*pi*F2*n)+(1/20)*randn(size(n)) ; %产生信号，由信噪比为10dB推出噪声功率；信号长度从X（1）到X（N+1）；
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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(F,PXX)
grid