% Function [ratio_pdf_mag,pdf_mag_signal,pdf_mag_noise]=MagPDFs_varT(p,s,sigma,T,Lim)
% estimates the conditional probability density functions of noisy wavelet
% coefficient magnitudes, given the presence and the absence of a "signal of interest"
%
% INPUT:
%    p, s - "shape" and "scale" parameters of the Laplacian prior ~exp(-|x|/s)^p
%            for noise-free data
%    sigma - noise standard deviation
%    T - threshold that defines the notion of significant noise-free magnitudes 
%        (i.e., threshold that defines the signal of interest)
%    Lim - maximum value for the magnitude of the wavelet coefficient in a given subband
%
% OUTPUT:
%    pdf_mag_signal, pdf_mag_noise - probability density functions of the noisy coefficient 
%    magnitudes dominated by signal and by noise, respectively;
%    ratio_pdf_mag - likelihod ratio: pdf_mag_signal/pdf_mag_noise
%
% Author: Aleksandra Pizurica/TELIN/Ghent University
% E-mail: Aleksandra.Pizurica@telin.UGent.be
% Last modified: June 2003.
% 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 [ratio_pdf_mag,pdf_mag_signal,pdf_mag_noise]=MagPDFs_varT1(p,s,sigma,T,Lim);

x=-Lim:Lim;

% pdf of pure noise
for i=1:length(x)
   pdf_purenoise(i)=exp(-x(i)^2/(2*sigma^2));
end
pdf_purenoise=pdf_purenoise/(sigma*sqrt(2*pi));


% estimated pdf of noise-free signal
for i=1:length(x)
   pdf_sigref(i)=exp(-abs(x(i))/s)^p;
   pdf_insig(i)=exp(-abs(x(i))/s)^p;
end

z=find(abs(x)<T);
pdf_sigref(z)=0;


z=find(abs(x)>T);
pdf_insig(z)=0;


if(sum(pdf_sigref)==0)
   
   ratio_pdf_mag(1:Lim+1)=0;
   pdf_mag_signal=zeros(Lim+1);
   pdf_mag_noise=zeros(Lim+1);
   
else
   
   pdf_sigref=pdf_sigref/(sum(pdf_sigref));   
   pdf_insig=pdf_insig/sum(pdf_insig);


   % pdf of significant signal + noise
   pdf_signal=conv(pdf_sigref,pdf_purenoise); 
   pdf_signal=pdf_signal(Lim+1:length(pdf_signal)-Lim);
   
   % pdf of insignificant signal + noise
   pdf_noise=conv(pdf_insig,pdf_purenoise); 
   pdf_noise=pdf_noise(Lim+1:length(pdf_noise)-Lim);


   % pdfs of coefficient magnitudes
   pdf_mag_noise=2*pdf_noise(Lim+1:2*Lim+1);
   pdf_mag_signal=2*pdf_signal(Lim+1:2*Lim+1);


   % to prevent numerical errors
   pdf_mag_noise_c=pdf_mag_noise;
   pdf_mag_signal_c=pdf_mag_signal;

   for i=1:Lim
      if(pdf_mag_noise(i)<1e-10)
        pdf_mag_noise_c(i+1:Lim+1)=pdf_mag_noise_c(i);
        pdf_mag_signal_c(i+1:Lim+1)=pdf_mag_signal_c(i);
      end
   end

   ratio_pdf_mag=pdf_mag_signal_c./pdf_mag_noise_c;

end
