function [s_hat,Count]=SBL(H,y,thresh)
%Sparse Bayesian learning method by D.P. Wipf and B.D. Rao.
% References
%   [1]  D.P. Wipf and B.D. Rao, “Sparse Bayesian learning for basis selection,” 
%       IEEE Trans. Signal Processing, vol. 52, pp. 2153–2164, Aug. 2004.
%   [2]  D.P. Wipf and B.D. Rao, “Comparing the effects of different weight distributions on finding sparse representations, ” 
%       in Advances in Neural Information Processing Systems, Y. Weiss, B. Scholkopf and J. Platt (Eds.), 
%       Cambridge MA: MIT Press, vol. 18, 2006, pp. 1521–1528.

%Coded by Kun Qiu, Aug. 2009
%Last updated, Dec. 2009

%%%%%%%%%Function Specification%%%%%%%%%%%%%
%INPUT:
%H:                    the sensing matrix
%y:                     the measurement column vector
%thresh:            threshold for quitting the algorithm (default 1e-9)
%OUTPUT:
%s_hat:              the signal esitmate
%Count:             the number of iterations

if nargin<2
    disp('Error in calling: not enough inputs');
    return;
end
if nargin<3
    thresh=1e-9;
end

[N,m]=size(H);

% Initialization
gamma=ones(m,1);
gamma_half=gamma.^(0.5);
pinv_HGamma=pinv(H*diag(gamma_half));
s_pre=pinv(H)*y;

% Main routine
for Count=1:1000000
    gamma=diag(s_pre*s_pre'+(eye(m)-diag(gamma_half)*pinv_HGamma*H)*diag(gamma));
    gamma_half=gamma.^(0.5);
    pinv_HGamma=pinv(H*diag(gamma_half));
    s_new=gamma_half.*(pinv_HGamma*y);
    if (s_new-s_pre)'*(s_new-s_pre)/m<thresh
        s_hat=s_new;
        break;
    else
        s_pre=s_new;
    end
end

s_hat=s_new;
