clc
clear all
close all
randn('seed',0);
lambda=5;
sigma2=1;
r=0.9;
N=100;
nreal=100;
for i=1:nreal
x0=sqrt(sigma2)*randn(1,N);
down=0;
up=0;
for n=1:N
down=down+r^(2*n);
end
A=sqrt(lambda*sigma2/down);
% A=up/down;
sn=A*r.^(1:N);
x1=x0+sn;
for n=1:N
up=up+r^n*x1(n);
end
sigma02=mean(x1.^2);
sigma12=0;
A1=up/down;
for n=1:N
sigma12=sigma12+(x1(n)-A1*r^n)^2;
end
T1(i)=N*log(N*sigma02./sigma12);
end
L=51;
for i=1:L
Pfaa(i)=(i-1)/(L-1);
u(i)=Qinv(Pfaa(i)/2);
Pda(i)=Q((u(i))-sqrt(lambda))+Q((u(i))+sqrt(lambda));
end
P=zeros(1,length(u));
u_new=sort(u,'descend');
for i=1:length(u_new)
clear Mgam
Mgam=find(T1>(u_new(i))^2);
disp(length(Mgam));
P(i)=length(Mgam)/nreal;
end
figure
subplot(1,2,1)
plot(Pfaa,Pda,Pfaa,P)
legend('理论值','蒙特卡洛结果');
xlabel('虚警概率PFA');
ylabel('检测概率PD');
title(['r=0.9,仿真次数100次']);
nreal=1000;
for i=1:nreal
x0=sqrt(sigma2)*randn(1,N);
down=0;
up=0;
for n=1:N
down=down+r^(2*n);
end
A=sqrt(lambda*sigma2/down);
% A=up/down;
sn=A*r.^(1:N);
x1=x0+sn;
for n=1:N
up=up+r^n*x1(n);
end
sigma02=mean(x1.^2);
sigma12=0;
A1=up/down;
for n=1:N
sigma12=sigma12+(x1(n)-A1*r^n)^2;
end
T1(i)=N*log(N*sigma02./sigma12);
end
L=51;
for i=1:L
Pfaa(i)=(i-1)/(L-1);
u(i)=Qinv(Pfaa(i)/2);
Pda(i)=Q((u(i))-sqrt(lambda))+Q((u(i))+sqrt(lambda));
end
P=zeros(1,length(u));
u_new=sort(u,'descend');
for i=1:length(u_new)
clear Mgam
Mgam=find(T1>(u_new(i))^2);
disp(length(Mgam));
P(i)=length(Mgam)/nreal;
end
subplot(1,2,2)
plot(Pfaa,Pda,Pfaa,P)
legend('理论值','蒙特卡洛结果');
xlabel('虚警概率PFA');
ylabel('检测概率PD');
title(['r=0.9,仿真次数1000次']);
评论2