burg算法估计功率谱 matlab 完全自编

burg算法估计功率谱 ，完全自编 ，没用matlab自带函数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%题目：                                                                    %
%已知AR（4）过程：                                                           %
%x（n）=2.760x（n-1）-3.809x（n-2）+2.654x（n-3）-0.924x（n-4）+w（n）                %
% 试用Burg算法利用Matlab估计模型参数、阶数和功率谱。                         %
%实验报告要求：实验仿真100次，但不要求给出所有100次结果，                     %
%只给出估计结果的方差和偏差（用估计结果均值与真值的差来估计偏差）。            %
%实验报告要求写清楚基本原理，提交程序并要有关键注释。                         %
%电子通信工程 研1205 于飞 2012020148                                        %
%2012年10月11日%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear allclose allclc;

sum=zeros（5001）;%psd估计值寄存器
for o=1:1:100
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %产生均值为零、方差为1、服从高斯（正态）分布的白噪声信号u（n）%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    w=randn（150000）;
    w=w-mean（w）;
    power_u=var（w）;
    m_u=mean（w）;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %使u（n）经过系统 x（n）=2.760x（n-1）-3.809x（n-2）+2.654x（n-3）-0.924x（n-4）+w（n）%
    %得到所要研究的随机过程x（n）%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    c=1;

    %d=[1-1.3521.338-0.6620.240];%宽带过程
    d=[1-2.7603.809-2.6540.924];%窄带过程
    y=filter（cdw）;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %burg算法计算参数、阶数、功率谱                                             %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    N=256;
    x=y（25000:25000+N-1）;%截取长度为256的样本
    f=x;%初始条件，前向预测值
    g=x;%初始条件，后向预测值
    p=（1/N）*x*x‘;%初始误差功率
    a=zeros（110）;%产生用来估计系数的向量a
    b=zeros（110）;%产生用来估计系数的向量b
    a（1）=1;%初始条件，初始系数a0=1;

    for m=1:1:10%循环递推，最多10阶
        num=-2*f（m+1:N）*g（m:N-1）‘;%计算反馈系数k的分子
        den=f（m+1:N）*f（m+1:N）‘+g（m:N-1）*g（m:N-1）‘;%计算反射系数k的分母
        k=num/den;%计算反射系数
    %%%%%%%%%%%%%%%%%%%%%%-----------计算系数a------------%%%%%%%%%%%%%%%%%%%%%%
        for n=1:1:m-1
           b（n）=a（m-n）;
        end;
        if m==1
            a（1）=k;
        else
            b（m）=1;
        end;

        a=a+k*b;
    %%%%%%%%%%%%%%%%%%%%%%------------------------------%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%-----------计算前、后项预测值-----------%%%%%%%%%%%%%%%
        ft=f;
        gt=g;
        f=ft+k*[0gt（1:N-1）];
        g=k*ft+[0gt（1:N-1）];
    %%%%%%%%%%%%%%%%%%%%%%------------------------------%%%%%%%%%%%%%%%%%%%
        pn=（1-k^2）*p;%计算误差功率
        if p-pn < 0.001%误差功率不再明显减小时，迭代结束
            break;
        end;
        p=pn;
    end;
    figure（1）
    a=[1a（1:m）];
    [HW]=freqz（1a500）;
    psd=power_u*（abs（H）.^2）;
    sum =sum +psd;
    %plot（（1/（2*pi））*W10*log10（psd））;%用估计的参数估计功率谱并绘出蓝色
    hold on;
    if o==1
        [HW]=freqz（cd500）;
        psdn=power_u*（abs（H）.^2）;
        plot（（

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2012-11-10 10:51  burgyf\
     文件        3493  2012-10-14 14:24  burgyf\ex2.asv
     文件        3537  2012-10-14 14:30  burgyf\ex2.m
     文件       69120  2012-11-08 22:25  burgyf\实验报告.doc