levinson durbin迭代算法

clc,clear
%产生零均值、方差为1的复高斯白噪声序列
N=256;
noise=(randn(1,N)+1i*randn(1,N))/sqrt(2);
%产生三个复正弦信号
f1=0.15;
f2=0.17;    %复正弦信号归一化频率
f3=0.26;
SNR1=30;
SNR2=30;    %复正弦信号信噪比
SNR3=27;
A1=10^(SNR1/20);
A2=10^(SNR2/20);     %复正弦信号幅度
A3=10^(SNR3/20);
signal1=A1*exp(1i*2*pi*f1*(0:N-1));
signal2=A2*exp(1i*2*pi*f2*(0:N-1));     %复正弦信号
signal3=A3*exp(1i*2*pi*f3*(0:N-1));
%产生观察样本
un=signal1+signal2+signal3+noise;
%计算自相关函数值
p=16;  %AR模型的阶数
r0=xcorr(un,p,'biased');  %直接计算自相关函数
r=r0(p+1:2*p+1);  %提取r(0),r(1),...,r(p)
%计算一阶AR模型的系数与输入方差
a(1,1)=-r(2)/r(1);  %1阶AR模型的系数
sigma(1)=r(1)-(abs(r(2))^2)/r(1); %1阶AR模型的输入方差
%Levinson-Durbin迭代算法的实现
for m=2:p
    k(m)=-(r(m+1)+sum(a(m-1,1:m-1).*r(m:-1:2)))/sigma(m-1); %反射系数
    a(m,m)=k(m);