压缩感知中OMP算法的C/C++实现

压缩感知中OMP算法的C/C++实现

• 背景介绍
• 算法实现部分
• 总结

---

阅读之前注意：

本文阅读建议用时：30min
本文阅读结构如下表：

项目	下属项目	测试用例数量
背景介绍	无	0
算法实现部分	无	1
总结	无	0

背景介绍

国内关于压缩感知的研究正在迅速发展，这里简要介绍一下压缩感知：所谓的压缩感知，就是把采样和压缩结合起来同时进行，这突破了传统的奈奎斯特采样定理（也叫香农定理），是信号采集理论的一次重大突破。压缩感知对于信号采集端的要求很低，并且只需要采集少量的信息，但采集信息后，在后台进行信号恢复时，需要进行大量的运算，因此信号重构算法和重构算法的GPU加速也是一大研究方向。
更多关于压缩感知，以及压缩感知的信号恢复重构算法，可以参考这两个博主：彬彬有礼的压缩感知专栏和Rachel Zhang的专栏，文章都写得很棒。所谓见贤思齐焉，我还得多多向他们学习！

算法实现部分

这一部分的算法实现（C/C++），改写自博客中的Matlab程序。
以下代码在配置Eigen库后，即可直接运行。

//本程序需要配置Eigen库
//OMP算法参考博客：https://blog.csdn.net/jbb0523/article/details/45130793
//程序改写自博客中的Matlab程序
//OMP_Version_1.0 2018.03.28_16:35 by Cooper Liu
//Questions ? Contact me : angelpoint@foxmail.com
#include"iostream"
#include <ctime>  
#include <Eigen/Eigenvalues>
using namespace std;
using namespace Eigen;

void swap(int *a, int i, int j)
{
	int tmp = a[i];
	a[i] = a[j];
	a[j] = tmp;
}

void randperm(int *a,int N)
{
	int i = 0;
	for (i = 0; i < N; i++)
		a[i] = i + 1;
	for (i = 0; i < N; i++)//交换a[i]和a[i]后面的随机序号
		swap(a, i, i + rand() % (N - i));
}

void main()
{
	int i = 0, j = 0;
	int M = 64, N = 256, K = 10;
	//int M = 64*5, N = 256*5, K = 10*5;
	int *Index_K = (int *)malloc(N*sizeof(int));
	randperm(Index_K, N);
	
	VectorXd X = VectorXd::Zero(N);//列向量
	for (i = 0; i < K; i++)
		X(Index_K[i]) = rand() / (double)(RAND_MAX);

	cout << "原始信号" << endl << X.transpose() << endl;
		
	MatrixXd Psi = MatrixXd::Identity(N, N);//单位阵
	MatrixXd Phi = (Eigen::MatrixXd::Random(M, N)).array();//随机高斯矩阵，保证任意2K列不线性相关
	//std::cout << Phi << std::endl;
	MatrixXd A(M, N);
	A = Phi*Psi;//A矩阵
	double miss = 0.;//检验A和Phi的差距
	for (i = 0; i < M; i++){
		for (j = 0; j < N; j++)
			miss = A(i, j) - Phi(i, j);
	}
	cout << "A和Phi miss is " << miss << endl;
	VectorXd y(M);
	y = Phi*X;//测量值
	//cout << y;

	//以下为OMP算法
	VectorXd theta = VectorXd::Zero(N);
	MatrixXd At = MatrixXd::Zero(M, K);
	VectorXd Pos_theta = VectorXd::Zero(K);
	VectorXd r_n = y;
	int pos = 0;
	VectorXd theta_ls;

	for (i = 0; i < K; i++)
	{
		VectorXd product = A.transpose()*r_n;
		VectorXd productTmp = product;
		productTmp = productTmp.cwiseAbs();
		productTmp.maxCoeff(&pos);
		for (j = 0; j < M; j++)
			At(j, i) = A(j, pos);
		Pos_theta(i) = pos;
		theta_ls = (At.leftCols(i+1).transpose()*At.leftCols(i+1)).inverse()*At.leftCols(i+1).transpose()*y;
		r_n = y - At.leftCols(i+1)*theta_ls;
	}
	for (i = 0; i < K; i++)
	{
		j = (int)Pos_theta(i);
		theta(j) = theta_ls(i);
	}
	
	MatrixXd x_r = Psi*theta;
	cout << "恢复的信号x_r is" << endl;
	cout << x_r.transpose() << endl;

	miss = 0.;
	for (i = 0; i < N; i++)
		miss = miss + abs(x_r(i) - X(i));
	cout << "info miss（信息损失） is " << miss << endl;
	system("pause");
}

总结

本篇博客的核心是：在理论上掌握OMP算法的流程，同时还应注意Eigen库中各种函数的调用，前一项要求线性代数的相关理论基础，后一项网上多搜索搜索就能掌握。
如果本文有帮助到你，不如请我一杯可乐吧 🍀