bt.m

ar模型 bt算法 ls rls lms music等数字信号处理的源代码

clear;
clc;
N1=10000;
N=1000;
a1=-0.195;
a2=0.95;
deltau=1;
tu=1
tv=tu*(1-a2)*((1+a2)^2-a1^2)/(1+a2)
tv=0.096525;
r=zeros(1,N);
v=sqrt(tv)*randn(1,N1);
u1=zeros(1,N1);
u1(1)=0;
u1(2)=0;
for m=1:N1-2
    u1(m+2)=-a1*u1(m+1)-a2*u1(m)+v(m);
end
u=u1(N1-N+1:N1);
clear u1 v;
for m=1:N
    for n=m:N
        r(m)=r(m)+u(n)*conj(u(n-m+1));
    end
end
r=r/N;
a1=-r(2)*(r(1)-r(3))/((r(1))^2-(r(2))^2)
a2=((r(2))^2-r(1)*r(3))/((r(1))^2-(r(2))^2)
sigma_v=r(1)+a1*r(2)+a2*r(3)
%%%%%%%%%%%%%%%   BT method    %%%%%%%%%%%%%
r1=[fliplr(r(2:256)) r(1:256)];
k=-100:100;
n1=-255:255;
P_BT=r1*exp(-j*pi/100).^(n1'*k);
plot(pi/100*k,10*log10(P_BT/max(P_BT)))
title('BT method PSD curve')
xlabel('w')
ylabel('normalized PSD(db)')
%%%%%%%%%%%%%%%   AR method    %%%%%%%%%%%%%
w=k*pi/100;
P_AR=(abs(1./(1+a1*exp(-j*w)+a2*exp(-2*j*w)))).^2*sigma_v;
figure
plot(w,10*log10(P_AR/max(P_AR)))
title('AR method PSD curve')
xlabel('w')
ylabel('normalized PSD(db)')