function output = blfdr(data, WaveletType, par, coarsest, gam)
% Shrinkage by Local False Discovery Rate in the Wavelet domain 
% Usage: output = blfdr(data, WaveletType, par, coarsest, gam)
% Input: 
%           data - data vector - 1-d signal;
%           WaveletType - string, 'Haar', 'Beylkin', 'Coiflet', 'Daubechies',
%                                'Symmlet', 'Vaidyanathan','Battle'
%           par  -  integer, it is a parameter related to the support and vanishing
%                         moments of the wavelets.
%           coarsest - Coarsest Level of V_0;  L << J
%           gam - parameter gamma in the formulat for determining level dependent prior
%           probabilities of H_0, gam 1.5 - 2.5 and higher work fine
%        
% Output:
%           output - 1-d denoised signal
%
%For reference please see "Local Bayesian False Discovery Rate Wavelet
%Shrinkage: Theory and Applications", 

N = length(data);  
num_lev = log2(N) - coarsest;
finest_level = log2(N) - 1;


%-----------------Wavelet decomposition
qmf = MakeONFilter(WaveletType, par);
wd = FWT_PO(data, coarsest, qmf);

  %---------------------------
    fin_ind=dyad(log2(N)-1);
    a=dyad(log2(N)-num_lev);
    sm_ind=1:a(1)-1; a=a(1);
    %--------mu----------!!!!!!!!!!!!!!!!---------
    finest_lev=wd(fin_ind);

    q1=prctile(finest_lev,25);
    q2=prctile(finest_lev,75);
    pseudos = abs(q2-q1)/1.5; 
    mu = 1/pseudos^2;      
    %----------------
    
    sigS=zeros(1,N);
    sigS(sm_ind) = wd(sm_ind);
    %------------------------------

   j = 1; 
   for i = finest_level:-1:coarsest
    
    PI(i) = 1 - 1/(i-coarsest +1).^gam;
    tau(i)=(max((var(data)-1/mu)/(2*(1-PI(i))^2),10^(-6)))^0.5;
    lev_ind=dyad(log2(N)-j);
    Bs = b_factor(wd(lev_ind), mu, tau(i));
    bi = find(Bs < 1);
    sigS(lev_ind(1)+ bi -1) =wd(lev_ind(1)+bi - 1);  
    j = j+1;
   end;
    %-------------------------------
       
   output = IWT_PO(sigS, coarsest,qmf);