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;