压缩感知重构算法之子空间追踪(SP)

function [ theta ] = CS_SP( y,A,K )
%CS_SP Summary of this function goes here 
%Version: 1.0 written by jbb0523 @2015-05-01 
%   Detailed explanation goes here 
%   y = Phi * x 
%   x = Psi * theta 
%   y = Phi*Psi * theta 
%   令 A = Phi*Psi, 则y=A*theta 
%   K is the sparsity level 
%   现在已知y和A，求theta 
%   Reference:Dai W，Milenkovic O．Subspace pursuit for compressive sensing 
%   signal reconstruction[J]．IEEE Transactions on Information Theory， 
%   2009，55(5)：2230-2249. 
    [y_rows,y_columns] = size(y); 
    if y_rows<y_columns 
        y = y‘;%y should be a column vector 
    end 
    [M,N] = size(A);%传感矩阵A为M*N矩阵 
    theta = zeros(N,1);%用来存储恢复的theta(列向量) 
    Pos_theta = [];%用来迭代过程中存储A被选择的列序号 
    r_n = y;%初始化残差(residual)为y 
    for kk=1:K%最多迭代K次 
        %(1) Identification 
        product = A‘*r_n;%传感矩阵A各列与残差的内积 
        [val,pos]=sort(abs(product),‘descend‘); 
        Js = pos(1:K);%选出内积值最大的K列 
        %(2) Support Merger 
        Is = union(Pos_theta,Js);%Pos_theta与Js并集 
        %(3) Estimation 
        %At的行数要大于列数，此为最小二乘的基础(列线性无关) 
        if length(Is)<=M 
            At = A(:,Is);%将A的这几列组成矩阵At 
        else%At的列数大于行数，列必为线性相关的,At‘*At将不可逆 
            break;%跳出for循环 
        end 
        %y=At*theta，以下求theta的最小二乘解(Least Square) 
        theta_ls = (At‘*At)^(-1)*At‘*y;%最小二乘解 
        %(4) Pruning 
        [val,pos]=sort(abs(theta_ls),‘descend‘); 
        %(5) Sample Update 
        Pos_theta = Is(pos(1:K)); 
        theta_ls = theta_ls(pos(1:K)); 
        %At(:,pos(1:K))*theta_ls是y在At(:,pos(1:K))列空间上的正交投影 
        r_n = y - At(:,pos(1:K))*theta_ls;%更新残差  
        if norm(r_n)<1e-6%Repeat the steps until r=0 
            break;%跳出for循环 
        end 
    end 
    theta(Pos_theta)=theta_ls;%恢复出的theta 
end