% Function [r,pdf_as,pdf_an]=AverageMagPDFs(pdf_s,pdf_n,W,Lim)
% estimates the conditional probability density functions for the locally averaged 
% coefficient magnitude, given the presence and the absence of a "signal of interest"
% where signal of interest is a noise-free magnitude that exceeds a certain threshold
%
% INPUT:
%    pdf_s, pdf_n - probability density functions of the noisy wavelet coefficients dominated by 
%    signal and by noise, respectively; 
%    W - "half window size", e.g., W=1 for 3x3 window and W=2 for 5x5 window 
%    Lim - maximum value for the magnitude of the wavlet coefficient in a given subband
% OUTPUT:
%    pdf_as, pdf_an - probability density functions of the locally averaged coefficient magnitude 
%    dominated by signal and by noise, respectively; 
%    r - likelihood ratio r=pdf_as/pdf_an
%
% Author: Aleksandra Pizurica/TELIN/Ghent University
% E-mail: Aleksandra.Pizurica@telin.UGent.be
% Last modified: June 2002.
% Copyright: All software, documentation, and related files in this distribution
%            are Copyright (c) 2003 Aleksandra Pizurica. All rights reserved.
%
% This software is distributed and licensed to you on a non-exclusive 
% basis, free-of-charge. Redistribution and use in source and binary forms, 
% with or without modification, are permitted provided that the above 
% copyright notice is clearly maintained.


function [r,pdf_as,pdf_an]=AverageMagPDFs(pdf_s,pdf_n,W,Lim);


Lim=length(pdf_s)-1;

if(sum(pdf_s)==0)
   
   r(1:Lim+1)=0;
   pdf_as=0;
   pdf_an=0;
   
else


% Number of coefficient magnitudes to be averaged
% neighbors of a central pixel in a a sqare window (2W+1)x(2W+1)
Count=(2*W+1)^2-1;



pdf_as=pdf_s;
pdf_an=pdf_n;


for i=1:(Count-1)
   pdf_as=conv(pdf_as,pdf_s);
   pdf_as=2*pdf_as(1:2:length(pdf_as));
   pdf_as=pdf_as/sum(pdf_as);
   
   pdf_an=conv(pdf_an,pdf_n);
   pdf_an=2*pdf_an(1:2:length(pdf_an));
   pdf_an=pdf_an/sum(pdf_an);
end



k=1;
while(pdf_an(k)==0)
   k=k+1;
end
r(1:k)=0;

pdf_an_c=1000*pdf_an;
pdf_as_c=1000*pdf_as;


% to prevent numerical errors

[Max,ind]=max(pdf_an);

for i=ind:Lim
   if(pdf_an(i)<1e-10)
      pdf_an_c(i+1:Lim+1)=pdf_an_c(i);
      pdf_as_c(i+1:Lim+1)=pdf_as_c(i);
   end
end
r(k:(Lim+1))=pdf_as_c(k:(Lim+1))./pdf_an_c(k:(Lim+1));

end


%figure
%subplot(2,1,1),plot(0:Lim,pdf_s,'k',0:Lim,pdf_as,'r');
%subplot(2,1,2),plot(0:Lim,pdf_n,'k',0:Lim,pdf_an,'r');

