function [ spectrum,peak_ind,count ] = Re_M_FOCUSS( A,y,p,lambda )
%输入参数
% A:过完备基,字典 | M*N维
% y:观测矢量 | M*L维
% p:p范数
% lambda:正则化参数
%计算维度
[M,N]=size(A);
[~,L]=size(y);
%默认参数设置
DISCARD_THRESHOLD=1e-4; %剔除0值的门限
END_THRESHOLD=1e-8; %结束迭代门限
MAX_ITERS=800; %最大迭代次数
%初始化
last_ans=ones(N,1); %初始化用于计算的x
x_sup=[1:N]'; %记录解x中非零值的位置,即x支撑集的坐标
sup_num=length(x_sup); %解中非零值个数
answear=zeros(N,L); %初始化解矩阵
count=0; %迭代次数
%迭代部分
while(1)
%剔除零值,减小计算量
if (min(last_ans)<DISCARD_THRESHOLD)%看解中有没有小于下限的值
if (~sup_num)
break;
end
index=find(last_ans>DISCARD_THRESHOLD);%把所有大于下限的值的位置拿出来,find()的返回值是第几个数,是坐标
last_ans=last_ans(index);%新的解长度变短了,仅有大于门限的元素
A=A(:,index);%同时也去掉小于门限值对应的导向矢量,仅保留新的解的对应部分
x_sup=x_sup(index);%更新支撑集坐标,这样不管我们怎么剔值,非零值在原谱上的位置是有记录的
sup_num=length(x_sup);
end
%计算结果
AA=A*diag(last_ans);
[U,S,V] = svd(AA,'econ');%econ是把多余的部分去掉,AA不是方阵时中间那个奇异值矩阵有0行或列,给去除了,奇异矢量也对应去除
%注意这里S,当sup_num>M时,S为M*M维,当sup_num<M时,S为sup_num*sup_num维(被econ给把多的去了)
diag_S=diag(S);
ans_last=answear;%answear用来放新的解,上一次的结果扔到ans_last中
answear=diag(last_ans)*V*diag(diag_S./(diag_S.^2 + lambda + 1e-16))*(U(:,1:min(M,sup_num)))'*y;
%更新权值
last_ans=sqrt((sum(abs(answear).^2,2)/L).^(1-p/2));
%检查结束条件
count=count+1;
if (count>MAX_ITERS)
break;
end
if (size(answear)==size(ans_last))%没再剔除0值
if (max(max(abs(answear-ans_last)))<END_THRESHOLD)%两次结果近似
break;
end
end
end
%输出结果
peak_ind=sort(x_sup);
spectrum=zeros(N,L);
spectrum(x_sup,:)=answear;
end
- 1
- 2
前往页