在matlab和opencv中分别实现稀疏表示

在本文中，稀疏表示的原理不再具体讲解，有需要的同学请自行百度。

本文采用OMP算法来求解稀疏系数。首先随机生成字典数据和待测试数据

字典数据：

dic =[  
   6,   7,   9,   9,   7,   0,   6,   3,   6,   9;
   1,   8,   7,   8,   5,   3,   8,   1,   7,   3;
   3,   3,   5,   4,   8,   2,   6,   1,   2,   2;
   6,   1,   0,   7,   3,   5,   0,   6,   3,   3;
   7,   5,   0,   5,   3,   0,   2,   7,   1,   7];

　　这是一个5*10的矩阵，行数代表维度，列数代表样本数。列数在字典中也叫字典原子，此处有10个原子，原子数大于维数，符合过完备要求。

信号数据：

signal=[  9;   8;   8;   3;   9];

　　为了简便，只模拟了一个信号数据，是一个5*1的矩阵，如果有多个数据，则应该是5*n的矩阵。求解的时候，可用循环求解。

一、在matlab中实现稀疏表示，求解稀疏系数

clc;close all;clear all;
dic =[  
   6,   7,   9,   9,   7,   0,   6,   3,   6,   9;
   1,   8,   7,   8,   5,   3,   8,   1,   7,   3;
   3,   3,   5,   4,   8,   2,   6,   1,   2,   2;
   6,   1,   0,   7,   3,   5,   0,   6,   3,   3;
   7,   5,   0,   5,   3,   0,   2,   7,   1,   7];    %字典
signal=[  9;   8;   8;   3;   9];    %原始信号 
dic=dic*diag(1./sqrt(sum(dic.^2)));  %字典原子单位化，即每列的norm为1
signal=signal/norm(signal);  %信号单位化
[A,res]=OMP(dic,signal,6);  %稀疏度设定为6，即非零元素最多为6个
A  %输出系数
res  %输出残差
epsilon=norm(signal-dic*A)  %验证残差  ||Y-Dx||2

其中OMP算法：

%OMP计算稀疏系数
function [A,res]=OMP(D,X,L)
% 输入参数:
%       D - 过完备字典，注意：必须字典的各列必须经过了规范化
%       X - 信号
%       L - 稀疏度，系数中非零元个数的最大值
% 输出参数:
%       A - 当前信号的系数
%       res - 残差

%%
residual=X; %初始化残差
indx=zeros(L,1);
for i=1:L,
    proj=D'*residual;%D转置与residual相乘，得到与residual与D每一列的内积值
    [~,pos]=max(abs(proj));%找到内积最大值的位置
    pos=pos(1);%若最大值不止一个，取第一个
    indx(i)=pos;%将这个位置存入索引集的第j个值
    a=pinv(D(:,indx(1:i)))*X;%indx(1:j)表示第一列前j个元素
    residual=X-D(:,indx(1:i))*a;
    res=norm(residual);
    if res< 1e-6
        break;
    end
end
A=zeros(size(D,2),1);
A(indx(indx~=0))=a;
end

　　

　　结果：

A =

    0.1450
    0.9391
         0
         0
    0.4210
    0.1049
         0
         0
   -0.5503
         0


res =

   3.1402e-16


epsilon =

   3.1402e-16

从系数中可以看出，从10个原子共选出了5个原子进行表示，最后的残差非常小，说明稀疏表示的结果和原数据非常接近。

二、在opencv2中实现稀疏表示 

void getData(Mat &data, Mat &signal);
int main(int argc, char* argv[])
{
    Mat dic, signal;
    getData(dic, signal);  //获取模拟数据
    Mat temp(1, dic.cols, CV_32F);  //用一个矩阵保存每个原子的模长
    for (int i = 0; i<dic.cols; i++)
    {
        temp.col(i) = norm(dic.col(i));  //每个原子的模长
    }
    divide(dic, repeat(temp, dic.rows, 1), dic); //字典原子单位化
    signal = signal / norm(signal);  //信号单位化
    Mat A=src.OMP(dic, signal, 8);  //调用OPM求解
    float res =(float)norm(signal - dic*A); //计算残差
    cout << "系数：" <<endl<< A << endl;
    cout<<endl<<"残差:"<< endl<<res << endl; //输出残差
    waitKey(0);
    return 0;
}
void getData(Mat &dic, Mat &signal)
{
    dic = (Mat_<float>(5, 10) <<
        6, 7, 9, 9, 7, 0, 6, 3, 6, 9,
        1, 8, 7, 8, 5, 3, 8, 1, 7, 3,
        3, 3, 5, 4, 8, 2, 6, 1, 2, 2,
        6, 1, 0, 7, 3, 5, 0, 6, 3, 3,
        7, 5, 0, 5, 3, 0, 2, 7, 1, 7);
    signal = (Mat_<float>(5, 1) << 9, 8, 8, 3, 9);
}

其中，OMP函数为：

Mat SRC::OMP(Mat& dic, Mat& signal,int sparsity)

{
    if (signal.cols>1)
    {
        cout << "wrong signal" << endl;
        exit(-1);
    }
    vector<int> selectedAtomOrder;   //保存所有选出的字典原子序号
    Mat coef(dic.cols, 1, CV_32F, Scalar::all(0)); //需要返回的系数    
    Mat residual = signal.clone();  //初始化残差
    Mat indx(0, 1, CV_32F);//初始化临时系数
    Mat phi;    //保存已选出的原子向量
    float max_coefficient;
    unsigned int atomOrder;  //每次所选择的原子的序号

    for (;;)
    {
        max_coefficient = 0;
        //取出内积最大列
        for (int i = 0; i <dic.cols; i++)
        {
            float coefficient = (float)dic.col(i).dot(residual);

            if (abs(coefficient) > abs(max_coefficient))
            {
                max_coefficient = coefficient;
                atomOrder = i;
            }
        }
        selectedAtomOrder.push_back(atomOrder); //添加选出的原子序号        
        Mat& temp_atom = dic.col(atomOrder); //取出该原子
        if (phi.cols == 0)
            phi = temp_atom;
        else
            hconcat(phi, temp_atom, phi); //将新原子合并到原子集合中（都是列向量）

        indx.push_back(0.0f);    //对系数矩阵新加一项
        solve(phi, signal, indx, DECOMP_SVD);    //求解最小二乘问题
        residual = signal - phi*indx;  //更新残差
        float res_norm = (float)norm(residual);
        if (indx.rows >= sparsity || res_norm <= 1e-6) //如果残差小于阈值或达到要求的稀疏度，就返回
        {
            for (int k = 0; k < selectedAtomOrder.size(); k++)
            {
                coef.row(selectedAtomOrder[k]).setTo(indx.row(k));  //得到最终的系数
            }
            return coef;
        }
    }
}

最终输出结果为：

系数：
[0.14503297;
 0.9391216;
 0;
 0;
 0.42096639;
 0.1048916;
 0;
 0;
 -0.55029994;
 0]

残差:
1.70999e-007

看以看出，opencv得到的系数和matlab得到的系数基本是一样，只是小数点后保留的位数区别。因为小数位数不相同，所以最后残差有点不同，但不影响最终结果，我们只需要系数相同即可。