pjh1.m

来自「PSD using BT/AR/ method,get the auto-cor」· M 代码 · 共 61 行

%auto-correlation
close all;
clear all;
N1=10000;	%总点数
N=1000;		%采样点数
a1=-0.195;
a2=0.95;
Pu=1;
var_v=(1-a2)*((1+a2)^2-a1^2)*(Pu^2)/(1+a2);
%get the noise-data accord N() distributing
V=sqrt(var_v)*randn(1,N1);
u1=zeros(1,N1);
u1(1)=0;
u1(2)=0;
%get the final 1000 data from the 10000 data
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;
%get the auto-correlation funtion
for m=1:N
Pre_r(m)=0;
for n=m:N
Pre_r(m)=u(n)*conj(u(n-m+1))+Pre_r(m);
end
Pre_r(m)=Pre_r(m)/N;
end
plot(Pre_r)
title(['Pre_r(0)=',num2str(Pre_r(1))]);

% estimate a1, a2, var_v
r0=Pre_r(1);
r=zeros(1,2);
r=Pre_r(2:3);
a1=-0.195;a2=0.95;
a1_est=-r(1)*(r0-r(2))/(r0^2-r(1)^2)
a2_est=(-r0*r(2)-r(1)^2)/(r0^2-r(1)^2)
var_v_est=r0+a1*r(1)+a2*r(2)

%PSD using BT method with estimated data
w=0:pi/200:2*pi;
Pre_Pbtsum=0;
for m=1:255
Pre_Pbtsum=Pre_r(m)*cos(w*m)+Pre_Pbtsum;
end
Pre_Pbt=Pre_Pbtsum*2+r0;      %r0=1
figure;
subplot(1,2,1),plot(w,abs(Pre_Pbt));
grid on;
title('PSD estimation of u(n) by BT method');
%--------------------------------------------------------------------------------------------------------------------
%PSD using AR method with estimated data
w=0:pi/200:2*pi;
Hw=1+a1_est*exp(-i*w)+a2_est*exp(-2*i*w); 
Hw=abs(1./Hw);
PAR=Hw.^2*var_v_est;
subplot(1,2,2),plot(w,PAR);
grid on;
title('PSD estimation of u(n) by AR method');