Matlab实现OMP（匹配追踪）

这个网上资源还是挺多哒。
我的作用仅是调试出来。
老规矩，注释必须全。抱拳。

%主文件
clc;
clear;
y=[7.6030 -4.1781 3.1123 1.0586 7.8053]';
A=[0.5377 -1.3077 -1.3499 -0.2050 0.6715 1.0347 0.8884;
    1.8339 -0.4336 3.3049 -0.1241 -1.2075 0.7269 -1.1471;
    -2.2588 0.3426 0.7254 1.4897 0.7172 -0.3034 -1.0689;
    0.8622 3.5784 -0.0631 1.4090 1.6302 0.2939 0.8095;
    0.3188 2.7694 0.7147 1.4172 0.4889 -0.7873 -2.9443];
CS_OMP(y,A,7);
%函数：CS_OMP
function [theta]=CS_OMP(y,A,t)
    %size可算出行和列分别是多少
    [y_rows,y_columns]=size(y);
    if y_rows<y_columns
        y=y';
    end
    [M,N]=size(A);                            %传感矩阵A为M*N矩阵
    theta=zeros(N,1);                         %用来存储恢复的theta(列向量)
    theta_ls=zeros(N,1);
    At=zeros(M,t);                            %用来迭代过程中存储A被选择的列
    Pos_theta=zeros(1,t);                    %用来迭代过程中存储A被选择的列序号
    r_n=y;                                    %初始化残差(residual)为y
    for ii=1:t                                  %迭代t次，t为输入参数
        fprintf('迭代次数：%d\n',ii);
        product=A'*r_n;                       %传感矩阵A各列与残差的内积
        fprintf('内积:');
        disp(product');                        %显示内积
        [value,pos]=max(abs(product));          %找到最大内积绝对值，即与残差最相关的列
        At(:,ii)=A(:,pos);                     %存储这一列
        fprintf('theta：\n')
        disp(At);
        Pos_theta(ii)=pos;                     %存储这一列的序号
        fprintf('已选择序列：');
        disp(Pos_theta);
        A(:,pos)=zeros(M,1);                   %清零A的这一列，其实此行可以不要，因为它与残差正交                                                        
        theta_ls=(At(:,1:ii)'*At(:,1:ii))^(-1)*At(:,1:ii)'*y;  %正交化
        fprintf('当前系数：');
        disp(theta_ls');                                                       
        r_n=y-At(:,1:ii)*theta_ls;             %更新残差,y一直是原始值
        fprintf('残差:');
        disp(r_n');
        error=sum(r_n.^2);                     %计算误差
        fprintf('误差=%d\n',error);
        if error < 1e-6                        %误差小于一定值结束迭代
            disp(ii);
            break;  
        end
    end
    for j=1:ii
        Pos(j)=Pos_theta(j);
    end
    fprintf('最终选定向量的下标有：')
    disp(Pos_theta);
    fprintf('系数集合：');
    disp(theta_ls');
    fprintf('归位后的系数：');
    theta(Pos)=theta_ls;  
    disp(theta');
 end