function [spectrum,cost] = MMV_ADM(A,Y,lamda,mu,pho)
% 该函数采用ADM方法解决如下问题
% min_X {||X||_2,1 + .5/lamda*||A*X-Y||_F^2}
% X为待稀疏恢复矩阵
% Y为多重观测向量
% lamda 为正则化常数
% mu为拉格朗日乘数
% pho 控制着敛散性
% 对于DOA估计，建议各参数设置如下
% lamda = 100 mu = 0.01 pho = 0.002(该参数不建议超过0.005，否则不收敛)
% 另外该方法对于相干信号效果不是很好，非相干效果比较好
% 参考论文 "A Fast Algorithm for Recovery of Jointly Sparse Vectors based on
% the Alternating Direction Methods" 2011
% 写于2023/3/27 肖武当
[M,N] = size(A);
[~,snap] = size(Y);
%初始化
X = A'*Y;
U = zeros(M,snap);
MAX_ITER = 500; %最大迭代次数
TOL_STOP = 1e-3; %误差界限
iter = 0;
old_X = X;
delta_X = inf;
% 声明一些中间需要变量，已提高运行速度
cost = zeros(MAX_ITER,1);
while (delta_X > TOL_STOP) && (iter < MAX_ITER)
iter = iter + 1;
%计算E
E = (lamda*mu/(1+lamda*mu))*(1/mu*U-(A*X-Y));
%计算梯度
G = A'*(A*X+E-Y-1/mu*U);
%更新X
X = Row_Shrink(X-pho*G,pho/mu);
%更新U
U = U-mu*(A*X+E-Y);
% 计算历史代价函数
for k = 1:N
cost(iter) = cost(iter) + sqrt(norm(X(k,:)));
end
delta_X = norm(X-old_X,'fro');
old_X = X;
end
spectrum = zeros(N,1);
for j = 1:N
spectrum(j) = norm(X(j,:));
end
end
