p506.m

现代功率谱估计中的参数估计AR模型 比经典谱估计有更好的分辨率

clear all
%program p506
%%%%%%%%%%%%%%%2阶AR模型u(n)+a1u(n-1)+a2u(n-2)=v(n) %%%%%%%%%%%
%%%%%%%%%%%%v(n)均值为0 方差为1 a1=-0.192    a2=0.92%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%产生随机过程的样本值%%%%%%%%%%%%%%%%%%%%%%%%%
a1=-0.192;
a2=0.92;
n=[1:1000]';
v0=randn(1,1000);
var=(1-a2)/(1+a2)*((1+a2).^2-a1.^2);
v=0.3*v0;
u0=[0 0];
num=1;
den=[1 a1 a2];
Zi=filtic(num,den,u0)
[u,Zf]=filter(num,den,v,Zi);
xlabel('n');
ylabel('V0(n)');
title('White Noise');
grid on;
legend('var=1',1)
subplot(221);
plot(n,v)
xlabel('n');
ylabel('V(n)');
title('White Noise');
grid on;
legend('var=0.3',1)
subplot(222)
plot(n,u)
xlabel('n');
ylabel('X(n)');
title('signal with white noise');
grid on;
%%%%%%%%%%%%%%%%计算r(m)的估计值,m=255.%%%%%%%%%%%%%%%%%%%%
lag=255;
[r,z]=xcorr(u,lag,'biased');
stem(z,r);
hndl=stem(z,r);
set(hndl,'Marker','o')
set(hndl,'Markersize',0.5)
xlabel('n');
ylabel('r(n)');
title('Autocorrelation');
grid on;
%%%%%%%%%%%%%画功率谱图%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=1:255;
k=-500:500;
w=(pi/500)*k;
R=k/500;
r0=0.9141;
r2=r(257:511)
X=r2*cos(pi/500*m'*k)+r0;
magX=abs(X);
subplot(223);
plot(R,magX);
xlabel('Angular frequency ,w(pi)');
ylabel('Pbt(w)');
title('Black_Turkey 功率谱曲线')
grid on;
%AR
HH=zeros(1,1001);
w=(pi/1000)*k;
a1=-0.195;
a2=0.95;
var=1;
H=1+a1*exp(-j*pi/500*k)+a2*exp(-j*2*pi/500*k);
magH=abs(H);
i=1;
while i<1001
    HH(i)=1/magH(i).^2;
    magH(i)=magH(i+1);
    i=i+1;
end
HH(1)=1/magH(1).^2;
HH(1001)=1/magH(1001).^2;
subplot(224)
plot(R,var*HH)
xlabel('Angular frequency,w(pi)');
ylabel('Par(w)');
title('AR功率谱曲线')
grid on;