块稀疏贝叶斯重构算法

function res = SBI（paras）

% res = SBI（paras）
% 
% SBI（paras） performs DOA estimation using Sparse Bayesian Inference
% 
% Input:
% paras.Y: M * T matrix sensor measurements at all snapshots
% paras.A: M * N matrix columns are the steering vectors for different directions
% paras.B: M * N matrix columns are derivatives of the steering vectors wrt. different directions
% paras.sigma2: initialization of noise variance
% paras.alpha: initialization of alpha
% paras.beta: initialization of beta
% paras.rho: rho
% paras.resolution: grid resolution for the directions
% paras.maxiter: maximum iteration
% paras.tol: stopping criterion
% paras.isKnownNoiseVar: true if known variance false if unknown
% paras.K: number of sources
% Output:
% res.mu: mean estimation
% res.Sigma: variance estimation
% res.sigma2: estimated noise variance
% res.sigma2seq: estimated noise variance at all iterations
% res.alpha: reconstructed alpha
% res.beta: reconstructed beta
% res.iter: iteration used in the algorithm
% res.ML: maximum likelihood function value at all iterations
% 
% Written by Zai Yang 19 Jul 2011
% reference: Zai Yang Lihua Xie and Cishen Zhang 
%    “Off-grid direction of arrival estimation using sparse Bayesian inference“

eps = 1e-16;

Y = paras.Y;
A = paras.A;
B = paras.B;
[M T] = size（Y）;
N = size（A 2）;

alpha0 = 1 / paras.sigma2;
rho = paras.rho / T;
beta = paras.beta;
alpha = paras.alpha;
r = paras.resolution;

maxiter = paras.maxiter;
tol = paras.tolerance;
if isfield（paras ‘isKnownNoiseVar‘） && ~isempty（paras.isKnownNoiseVar）
    isKnownNoiseVar = paras.isKnownNoiseVar;
else
    isKnownNoiseVar = false;
end

if isKnownNoiseVar
   a = 1;
   b = T * M * paras.knownsigma2;
else
    a = 1e-4;
    b = 1e-4;
end

  if isfield（paras ‘K‘） && ~isempty（paras.K）
       K = paras.K;
  else
      K = min（T M-1）;
    end
   
idx = [];
BHB = B‘ * B;
converged = false;
iter_beta = 1;
iter = 0;
ML = zeros（maxiter1）;
alpha0seq = zeros（maxiter1）;

while ~converged
    iter = iter + 1;
   
    Phi = A;
    Phi（:idx） = A（:idx） + B（:idx） * diag（beta（idx））;
    
    alpha_last = alpha;
    
    C = 1 / alpha0 * eye（M） + Phi * diag（alpha） * Phi‘;
%     Sigma = diag（alpha） - diag（alpha） * Phi‘ / C * Phi * di