注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Mr.Right

不顾一切的去想,于是我们有了梦想。脚踏实地的去做,于是梦想成了现实。

 
 
 

日志

 
 
关于我

人生一年又一年,只要每年都有所积累,有所成长,都有那么一次自己认为满意的花开时刻就好。即使一时不顺,也要敞开胸怀。生命的荣枯并不是简单的重复,一时的得失不是成败的尺度。花开不是荣耀,而是一个美丽的结束,花谢也不是耻辱,而是一个低调的开始。

网易考拉推荐

Monte carlo simulation (蒙特卡罗模拟) 完整范例  

2012-03-17 11:30:43|  分类: 学习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
% 将下面代码保存为PdPfa.m,然后Run
% PdPfa.m
%
%  Pd/Pfa via analysis and simulation for
%  nonfluctuating target and Rayleigh interference
%
clear, hold off
close all
format compact
j = sqrt(-1);

% Form complex interference signal with Rayleigh amplitude
% and unit mean

N = input('Enter desired noise sequence length: ');     %  输入500000
snrdB = input('Enter desired SNR (dB): ');                      % 输入 13

beta = 2/sqrt(pi)*1;    %  beta是噪声的振幅 beta = 1.1284
cnoise = beta/sqrt(2)*randn(N,1)+j*beta/sqrt(2)*randn(N,1); % 复高斯噪声,均值为beta
ray = abs(cnoise);
mean_ray = mean(ray)  % just checking to make sure mean ~= 1 , 在 beta=1.1284 左右波动

% Now create another sequence with a target
% added to each sample at the specified SNR.

npow = beta^2;    % n平方,beta^2表示噪声功率
S = sqrt((10^(snrdB/10))*npow);   % 给定噪声功率npow,如何加入snrdB的信号S

snoise = cnoise + S;

sray = abs(snoise);

% The matched filter output is simply the input multiplied
% by conj(S) (but S is real here); the test statistic is the magnitude of that
% quantity (see Fig. 6-5 in text).  However, since we usually don't
% know S in advance, I will assume S = 1 for scaling the data and
% establishing the threshold.

% Now let's set up a series of threshold settings to achieve
% desired Pfa's using analytic formula from text

Pfa = [1e-7 2e-7 5e-7 1e-6 2e-6 5e-6 1e-5 4e-5 1e-4 4e-4 1e-3  1e-2]';
T = beta*sqrt(-log(Pfa));    % 注意门限 T = S*beta*sqrt(-log(Pfa)),这里由于假设 S = 1

% Now let's measure the Pfa with these thresholds and see how
% the results stack up against the design goal

for k=1:length(T)
  Pfa_sim(k) = sum( ray>T(k) )/N;
end
fprintf('\n\n Threshold    Theoretical Pfa    Simulated Pfa\n')
for k=1:length(T)
  fprintf('%9.6g          %5.3g         %5.3g\n',T(k),Pfa(k),Pfa_sim(k))
end
Monte carlo simulation (蒙特卡罗模拟) 完整范例 - 阿英 - Mr.Right
 
figure
loglog(Pfa, Pfa_sim);xlabel('Theoretical Pfa');ylabel('Simulated Pfa');
grid


% Now let's predict Pd for the same thresholds analytically,
% and again compare against simulation

% Do two versions, one using Marcum Q to predict Pd, the other
% using Albersheim with N=1 (and rearranged to compute Pd instead
% of SNR)

% Some caution needed here in using the marcum function. For the
% simulation, where I was mimicking the actual processing of data, I
% couldn't assume I knew S, so I just used S = 1 in the threshold
% calculation.  To use the Marcum function, however, I do have to plug in
% the actual SNR, which means using the actual value of S, and I need to
% use the threshold I would have computed if I had known S.  That threshold
% would be larger than my actual threshold by a factor of S, so where the
% Marcum function uses sqrt(2)*T(k)/S/beta for the second argument, I just
% want to use sqrt(2)*T(k)/beta, because my threshold is already smaller by
% the factor S.

for k=1:length(T)
   A=log(0.62/Pfa(k));
   Z=snrdB/(6.2+(4.54/sqrt(1.44)));
   B=(10^Z-A)/(1.7+0.12*A);
   Pd_Alb(k)=1/(1+exp(-B));
   
   Pd_Q(k)=marcum(sqrt(2)*S/beta,sqrt(2)*T(k)/beta,0.001);
   
   Pd_sim(k) = sum( sray>T(k) )/N;
end
fprintf('\n\n Threshold    Pd (Marcum Q)    Pd (Albersheim)     Simulated Pd\n')
for k=1:length(T)
  fprintf('%9.6g        %6.4g              %6.4g              %6.4g\n',T(k),Pd_Q(k),Pd_Alb(k),Pd_sim(k))
end

% Finally, let's plot Pd vs. Pfa in  "receiver operating characteristic"
% (ROC) curve.  Plot both predicted and simulated Pd.
Monte carlo simulation (蒙特卡罗模拟) 完整范例 - 阿英 - Mr.Right
 
figure
semilogx(Pfa,Pd_Alb,Pfa,Pd_Q,Pfa,Pd_sim); grid; xlabel('Pfa'); ylabel('Pd');
legend('Albersheim','Marcum Q','Simulated',-1)
title(['SNR=',int2str(snrdB),' dB'])

figure
hist([ray,sray],100)
Monte carlo simulation (蒙特卡罗模拟) 完整范例 - 阿英 - Mr.Right
 

% -----------------------------------------------
% 用到的函数
function q=marcum(a,t,prec)
% For (x,y) ~ N(0,a,1,1,0) marcum(a,t,prec) computes the 
%    probability that x^2+y^2 > t^2 using a method from 
%    Brennan & Reed, IEEE Trans. Info Theory, 1965, 
%    pp. 312-313
 
% Modified to avoid infinite loops
%   for very small probabiliities ( q < 1e-6), which occur for large values
%   of 'a', which occur when the SNR in the calling function
%   is large
%
% prec=precision of computation. prec=0.001 => compute Q to 0.1% accuracy.
%

ct = .5*t^2;
ca = .5*a^2;
et = exp(-ct);
g  = 1 - et;
k  = exp(-ca);
gnew = et;
q = 1 - g*k;

n=0;
while ( n < a*t/2) | ( g*k/(1-(.5*a*t/n)^2) > prec*q )
    if (q < 1e-6)
        return
    end
n = n + 1;
gnew = gnew * ct/n;
g = g - gnew;
k = k * ca/n;
q = q - g*k;
end

  评论这张
 
阅读(1045)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016