function output = bafdr(data, WaveletType, par, coarsest, alpha, pi_0)
% Shrinkage by Ordering of Posterior Probabilities (BaFDR)
% Usage: output = bafdr(data, WaveletType, par, coarsest, alpha, pi_0)
% 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
%           alpha -  Pr( type I error) - significance level, set alpha to be small 0.05 for example
%           pi_0  - prior probability of H_0 ( wavelet coef = 0)  set to be
%                     close to 1, pi_0  = 0.90 or 0.95 for example
%        
% 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;
pi_1 = 1- pi_0;
sig = zeros(1,N);


%-----------------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);
    %---------------------------
    sig(sm_ind) = wd(sm_ind);
    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;        
    tau=(max((var(data)-1/mu)/(2*(1-pi_0)^2),10^(-6)))^0.5;


    Bs=b_factor(wd(a:N),mu,tau);
    Ps = Bs./(pi_1/pi_0 + Bs);

    [P_sort, Ind_sort]=sort(Ps);
    EQ=0; R=1;

    while EQ<alpha
        R=R+1;
        P0ly=[1-P_sort(1), P_sort(1)];
        for i=2:R
            b=[1-P_sort(i), P_sort(i)];
            P0ly=conv(P0ly,b);
        end;
        P0ly1=P0ly(2:R+1);
        zz=1:R;
        P0ly1=zz.*P0ly1;
        EQ=sum(P0ly1)/R;
        EQQ(R)=EQ;
    end;
    R=R-1;
    s_ind=Ind_sort(1:R);
    sig(s_ind+(a-1))=wd(s_ind+(a-1));
    output=IWT_PO(sig,coarsest,qmf);