www.pudn.com > Levenson.rar > Levenson.m, change:2014-11-15,size:1246b

```clear all;
clear
tic;
fs=1;
N=input('the value of length:');
f1=0.2*fs;
f2=0.3*fs;
M = 16;
L = 2*N;
n=1:N;
s = sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs);
x = awgn(s,10);
subplot(2,2,1);
plot(s,'b'),xlabel('时间'),ylabel('幅度'),title('原始信号s');
grid;
subplot(2,2,2);
plot(x,'r'),xlabel('时间'),ylabel('幅度'),title('观测信号x');
grid;

rxx = xcorr(x,x,M,'biased');
r0 = rxx(M+1);
R = rxx(M+2:2*M+1);
%Levinson 递推算法
a = zeros(M,M);
FPE = zeros(1,M);
var = zeros(1,M);
a(1,1) = -R(1)/r0;
var(1) = (1-(abs(a(1,1)))^2)*r0;
FPE(1) = var(1)*(M+2)/(M);
for p=2:M
sum=0;
for k=1:p-1
sum=sum+a(p-1,k)*R(p-k);
end
a(p,p)=-(R(p)+sum)/var(p-1);
for k=1:p-1
a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);
end
var(p)=(1-a(p,p)^2)*var(p-1);
FPE(p)=var(p)*(M+1+p)/(M+1-p);
end
min=FPE(1);
p = 1;
for k=2:M
if FPE(k)<min
min=FPE(k);
p=k;
end
end
W=0.01:0.01:pi;
he=ones(1,length(W));
for k=1:p
he=he+(a(p,k).*exp(-j*k*W));
end
Pxx=var(p)./((abs(he)).^2);
F=W*fs/(pi*2);
subplot(2,1,2);
plot(F,abs(Pxx))；
xlabel('频率/Hz'),ylabel('功率谱P');
title(['AR模型的最佳阶数p=' int2str(p)]);
grid;
toc;
```