在本文中，稀疏表示的原理不再具體講解，有需要的同學請自行百度。
本文采用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得到的系數基本是一樣，只是小數點后保留的位數區別。因為小數位數不相同，所以最后殘差有點不同，但不影響最終結果，我們只需要系數相同即可。