function [median,sigma] = plotNR(x,y,x1,x2)
% description: function plotNR plots cumulative distribution functions 
% on probability paper (normal y-axis)
%
% usage: plotNR(x,y,x1,2), x is a range of thresholds, y is a range of
% cumulative probabilities from 0 to 1, and x1 & x2 are the threshold
% limits
%
% plotNR also performs a regression fit on the data and determines the
% median and sigma (both in dB) of the best-fit Gaussian PDF

%-------------------------------------------------------------------------

for i = length(y): -1: 1                     % remove the portions of the 
                                             % CDF corresponding to pr =0
   if((y(i)>=1-1E-9)|(y(i)<=0))              % or pr = 1
       x(i) = []; y(i) = [];
   end 

end   

%--------------------------------------------------------------------------

yp=invQ(y);        % distort the "y" component of the function so that
		   % a Gaussian CDF will plot as a straight line in the 
		   % physical coordinate space of the graph

%--------------------------------------------------------------------------

c=polyfit(x,yp,1); % use regression to find the coefficients "c" of the 
                   % polynomial of order 1 that best fits the ordered 
		   % pairs (x,yp)
		   
yfit=polyval(c,x); % generate a set of points which pass through the 
                   % regression line defined by the coefficients "c"
	    
%--------------------------------------------------------------------------
		   
plot(x,yp,'-',x,yfit,'--'); % plot the data points with a solid line and
                            % the regression points with a dashed line

%--------------------------------------------------------------------------

X90=(invQ(0.9)-c(2))/c(1);  % y = mx + b  => x = (y-b)/m
X50=(invQ(0.5)-c(2))/c(1);  % 

median=X50;                 % calculate the median and sigma of the
sigma=(X90-X50)/1.29;       % best fit Gaussian distribution

%--------------------------------------------------------------------------   
% tic locations for the vertical axis

  % if the desired probability limits are [0.0001...0.99999]: 

  y1=[0.0001 0.001 0.01 0.05 0.1  0.2  0.3  0.4  0.5...
      0.6    0.7   0.8  0.9  0.95 0.99 0.999 0.9999];

  % if the desired probability limits are [0.01...0.99]:
      
  %  y1 =[0.01 0.05 0.1 0.2 0.3  0.4 0.5...
  %	0.6  0.7  0.8 0.9 0.95 0.99 ];    

z=invQ(y1);   % physical locations for the y tic marks

%--------------------------------------------------------------------------   
   
a=gca;                             % handle for the current plot

set(a,'YTick', z)                  % set locations for yticks
set(a,'YTickLabel', y1)            % set values for yticks

set(a,'YLim',[z(1) z(length(z))])  % set physical min and max for y-axis 
set(a,'XLim',[x1 x2])              % set physical min and max for x-axis

set(a,'YGrid','on') 		   % turn on y grid
set(a,'XGrid','on')                % turn on x grid

% the values for Title, Legend, XLabel, YLabel, etc. are normally set by the
% calling program

%--------------------------------------------------------------------------
