
% Function [p,s]=FindGenLapPar(D,sigma_n)
% estimates the parameters of the generalized Laplacian prior exp(-|x|/s)^p 
% for noise-free wavelet coefficients from the noisy histogram
%
% INPUT:
%    D- noisy subband
%    sigma_n - noise standard deviation in the given subband
%
% OUTPUT:
%    p, s - parameters of the generalized Laplacian prior for noise-free data
%
% 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 [p,s]=FindGenLapPar(D,sigma_n);



[H,V]=size(D);
h=reshape(D,1,H*V);

m4=mean(h.^4);
sigma_s=std(h);

Q=(m4+3*sigma_n^4-6*sigma_s^2*sigma_n^2)/(sigma_s^2-sigma_n^2)^2;

% For the generalized Laplacian distribution A*exp(-abs(x/s)^p) :
% gamma(1/p)*gamma(5/p)/(gamma(3/p)^2) = Q

load Lap_p.mat
load Lap_kurtosis.mat

i=0;
Diff=1000000;
fitting='y';
while(fitting=='y')
   i=i+1;
   K=Lap_kurtosis(i);
   Diff_old=Diff;
   Diff=abs(Q-K);  
   if (i>size(Lap_p)-1)
           fitting='n';
   end
   if((Diff>Diff_old))
      fitting='n';
   end
end

p=Lap_p(i);


s=sqrt((sigma_s^2-sigma_n^2)*gamma(1/p)/gamma(3/p));


