高斯信道衰减-MATLAB程序.doc

clear all;
close all;

echo on;
display('Beginning');
echo off;

N=50;
SNRindB=-5;
snr=10.^(SNRindB./10);
a=[0.005 0.01 0.025 0.05 0.1 0.15 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 0.98];
CV=[4.559 3.857 3.070 2.492 1.933 1.610 1.405 1.126 0.924 0.773 0.650 0.545 0.448 0.346 0.282 0.227];
% a=[0.98 0.96 0.93 0.9 0.8 0.6 0.5 0.4 0.3 0.2 0.1 0.05 0.01 0.005 0.001];
exp_time=10000;
%%%%%%%%%%%%%%%%%M是高斯信号%%%%%%%%%%%%%%%%%%%%%%%
for i=1:length(a);
display(i);
count=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
ave_segma_MLE=0;
for k=1:exp_time
% display(k); V
%%%%产生高斯随机数，并进行处理%%%%%%%%
R1=randn(1,N);
R2=randn(1,N);
R=(sqrt(snr).*R1+R2);
Y=sort(R);
%%%%%%%%%%%%%进行极大似然估计%%%%%%%%%%%
% R3=sum(R)/N;
% for w=1:N;
% A(w)=(R(w)-R3).^2;
% end
% segma_MLE=sqrt(sum(A)/N);
segma_MLE=1;
temp2=0;
for j=1:1:N

%%%%假设的高斯分布函数F（x）
%%%Z_i
h1=@(t)((1/(sqrt(2*pi).*segma_MLE)).*exp((t.^2)/(-2*(segma_MLE.^2))));
F_Y1=quadl(h1,-10000,Y(j));
Z_1=F_Y1;

%%%%Z_(N-i+1)
h2=@(t)((1/(sqrt(2*pi).*segma_MLE)).*exp((t.^2)/(-2*(segma_MLE.^2))));
F_Y2=quadl(h2,-10000,Y(N-j+1));
Z_2=F_Y2;

%%%%计算差值
temp2=temp2+(2*j-1)*(log(Z_1)+log(1-Z_2));
% temp2=temp2+(2*j-1)*(log(Z_1))+(2*N+1-2*j)*(log(1-Z_1));
end
%A_n2=-N-temp2;
A_n2=-N-(1/N)*temp2;
%display(A_n2);

if A_n2>=CV(i)
count(i)=count(i)+1;
end
count1(i)=count(i);
end
P1(i)=count1(i)/exp_time;
display(P1(i));
% ave_segma_MLE=ave_segma_MLE/exp_time;
% display(ave_segma_MLE);
end
%%%%%%%%%%%%%%%%%%%%M是正弦信号%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:length(a);
display(i);
count=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
ave_segma_MLE=0;
for k=1:exp_time
% display(k);
%%%%产生正弦信号随机数，并进行处理%%%%%%%%
R=linspace(0,2*pi,N);
for f=1:N;
Y1(f)=sqrt(2)*sin((pi/3)*R(f)+(pi/3));
end
%%%%产生高斯分布随机数，并进行排序处理
%X=raylrnd(segma_MLE,[1 N]);
%Y=sort(X);
R1=randn(1,N);
Y2=(sqrt(snr).*Y1+R1);
Y=sort(Y2);
%%%%%%%%%%%%%进行极大似然估计%%%%%%%%%%%
% R2=sum(Y)/N;
% for w=1:N;
% A(w)=(Y(w)-R2).^2;
% end
% segma_MLE=sqrt(sum(A)/N);
segma_MLE=1;
temp2=0;
for j=1:1:N

%%%%假设的高斯分布函数F（x）
%%%Z_i
h1=@(t)((1/(sqrt(2*pi).*segma_MLE)).*exp((t.^2)/(-2*(segma_MLE.^2))));
F_Y1=quadl(h1,-10000,Y(j));
Z_1=F_Y1;

%%%%Z_(N-i+1)
h2=@(t)((1/(sqrt(2*pi).*segma_MLE)).*exp((t.^2)/(-2*(segma_MLE.^2))));
F_Y2=quadl(h2,-10000,Y(N-j+1));
Z_2=F_Y2;

%%%%计算差值
temp2=temp2+(2*j-1)*(log(Z_1)+log(1-Z_2));
% temp2=temp2+(2*j-1)*(log(Z_1))+(2*N+1-2*j)*(log(1-Z_1));
end
%A_n2=-N-temp2;
A_n2=-N-(1/N)*temp2;
%display(A_n2);

if A_n2>=CV(i)
count(i)=count(i)+1;
end
count1(i)=count(i);
end
P2(i)=count1(i)/exp_time;
display(P2(i));
% ave_segma_MLE=ave_segma_MLE/exp_time;
% display(ave_segma_MLE);
end
%%%%%%%%%%%%%%%%%%%%%%%%M是1%%%%%%%%%%%%%%%%%%%%%%%%%%%%555
for i=1:length(a);
display(i);
count=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
ave_segma_MLE=0;
for k=1:exp_time
% display(k);
%%%%产生高斯随机数，并进行处理%%%%%%%%
R1=randn(1,N);
%%%%%m=1%%%%%%%%%%%%%%%
R2=ones(1,N);
R=(sqrt(snr).*R2+R1);
Y=sort(R);
%%%%%%%%%%%%%进行极大似然估计%%%%%%%%%%%
% R3=sum(R)/N;
% for w=1:N;
% A(w)=(R(w)-R3).^2;
% end
% segma_MLE=sqrt(sum(A)/N);
segma_MLE=1;
temp2=0;
for j=1:1:N

%%%%假设的高斯分布函数F（x）
%%%Z_i
h1=@(t)((1/(sqrt(2*pi).*segma_MLE)).*exp((t.^2)/(-2*(segma_MLE.^2))));
F_Y1=quadl(h1,-10000,Y(j));
Z_1=F_Y1;

%%%%Z_(N-i+1)
h2=@(t)((1/(sqrt(2*pi).*segma_MLE)).*exp((t.^2)/(-2*(segma_MLE.^2))));
F_Y2=quadl(h2,-10000,Y(N-j+1));
Z_2=F_Y2;

%%%%计算差值
temp2=temp2+(2*j-1)*(log(Z_1)+log(1-Z_2));
% temp2=temp2+(2*j-1)*(log(Z_1))+(2*N+1-2*j)*(log(1-Z_1));
end
%A_n2=-N-temp2;
A_n2=-N-(1/N)*temp2;
%display(A_n2);

if A_n2>=CV(i)
count(i)=count(i)+1;
end
count1(i)=count(i);
end
P3(i)=count1(i)/exp_time;
display(P3(i));
% ave_segma_MLE=ave_segma_MLE/exp_time;
% display(ave_segma_MLE);
end

for i=1:length(a),
if a(i)>0.1,
aa=N/2;
else
aa=N;
end
threshold(i)=2*fzero(@(x)threshold_pf(x,N/2,a(i)),aa)/N;
thr=threshold(i)*N;
sr=snr*N;
Pd_AWGN(i)=marcumq(sqrt(sr),sqrt(thr),N/2);
end
% display(a);
% display(P1,P2,P3);
plot(a,Pd_AWGN,'b-o',a,P1,'b*-',a,P2,'ks-',a,P3,'ro-');
axis([1e-3 1 1e-3 1]);
xlabel('a:false alarm Probability');
ylabel('P:Probability of detection');
temp3=['ED N=50 SNR=-5 '];
temp4=['AD m=awgn N=50 SNR=-5 '];
temp5=['AD m=sin N=50 SNR=-5 '];
temp6=['AD m=1 N=50 SNR=-5 '];
legend(temp3,temp4,temp5,temp6);
grid;

% loglog(a,P,'b--');
% axis([1e-3 1 1e-3 1]);
% %temp0=['AWGN local detection, N=1'];
%
% %temp4=['AWGN AND rule, K=',num2str(N),', theory'];
% %temp5=['AWGN Majority rule, K=',num2str(K),', theory'];
% %temp6=['AWGN OR rule, K=1, theory'];
%
% %legend(temp4,temp5,temp6);
% xlabel('a:false alarm Probability');
% ylabel('P:Probability of detection');
% grid;