目录
- 前言
- 一、压缩感知基本原理
- 二、代码仿真
- 1. CVX工具箱求解L1范数
- 2. CVX学习视频
- 3. 仿真实现
- 三、 重点参考
前言
压缩感知（Compressive Sensing,CS）与传统的香农采样定理（奈奎斯特采样定理）有着明显区别，香农采样定理明确采样率应为信号带宽的2倍以上，一般工程上取2.56倍到5倍，采样率越高数据量越大。
一、压缩感知基本原理
压缩感知（Compressive Sensing，CS）。相对于传统的奈奎斯特采样定理——要求采样频率必须是信号最高频率的两倍或两倍以上（这就要求信号是带限信号，通常在采样前使用低通滤波器使信号带限），压缩感知则利用数据的冗余特性，只采集少量的样本还原原始数据
总结：压缩感知方法是利用重建信号在变换域上的稀疏作为先验信息，用测量矩阵观测被测信号，由观测值结合重建算法重建出完整的被测信号。
1、目标信号稀疏表示。将目标信号变换到尽可能稀疏的稀疏变换域，即为信号的稀疏表示。常用的稀疏基有离散余弦变换基（DCT）,快速傅里叶变换（FFT）.
2、构建测量矩阵。测量矩阵即为压缩感知采样的实现部分，通过测量矩阵控制采样过程，采样量下降同时保证目标信号所含有效信息不丢失，能够由压缩采样值重建出目标信号。常用的有高斯随机矩阵，部分哈达玛矩阵等。
3、设计重建算法。重建算法是从采样值求解最优化问题寻找到目标信号最优解。重建算法的准确性，高效性和稳定性是其设计的关键。常用的重建算法有L1范数和正交匹配追踪算法（OMP）等。
二、代码仿真
1. CVX工具箱求解L1范数
CVX是一个基于matlab的凸优化建模系统。CVX将Matlab转换为一种建模语言，允许使用标准的Matlab表达式语法指定约束和目标。例如，考虑以下凸优化模型:
在默认模式下，CVX支持一种特殊的凸优化方法，我们称之为有纪律的凸编程。在这种方法下，凸函数和凸集是从凸分析的一个小规则集出发，从凸函数和凸集的一个基本库出发，建立起来的。使用这些规则表示的约束和目标会自动转换为规范形式并得到解决。要了解更多关于约束凸编程的信息，请参阅这些参考资料;有关凸分析和凸优化的基础知识，请参阅《凸优化》一书。
注意CVX不是非线性优化的通用工具，也不是检查你的模型是否凸的工具。在开始使用CVX之前，务必确认您的模型可以表示为MIDCP或GP。如果两者都不是，那么CVX就不是这个任务的正确工具。
官方Introduction地址：
http://cvxr.com/cvx/doc/intro.html#what-is-cvx
下载地址
CVX Home Page：
http://cvxr.com/cvx/
下载地址：
http://cvxr.com/cvx/download/
我们需要的操作是：
1.从官网下载CVX的安装包；
2.解压到自己想要的文件夹中，此时会生成一个cvx的文件夹。
3.打开MATLAB，不要手动将CVX添加到路径中。
4.在MATLAB命令行输入下列命令：
cd C:\personal\cvx
cvx_setup
这里cvx_setup函数执行各种任务来验证您的安装是正确的，设置Matlab搜索路径，以便它可以找到所有的CVX程序文件，并运行一个简单的测试问题来验证安装。
2. CVX学习视频
CVX工具箱学习教程
b站有位老师的入门视频讲得很好：
https://www.bilibili.com/video/BV1UQ4y1K7Vf?spm_id_from=333.999.0.0
3. 仿真实现
本文分别以稀疏基有离散余弦变换基（DCT）和快速傅立叶变换基（FFT）做为稀疏基，高斯随机矩阵、部分哈达玛矩阵为测量矩阵，L1范数、正交匹配追踪算法（OMP）为重建算法进行压缩感知算法实现。
本文以f = cos(2pi/256t) + sin(2pi/128t)做为原信号，取原信号f的20%做为输入进行压缩感知重建。
main.m
% 该程序用于验证压缩感知理论（包含了L1最小范数求解和OMP求解）
%
%
%
clear all; close all;
%% 产生信号
choice_transform = 1; % 选择正交基，1为选择DCT变换，0为选择FFT变换
choice_Phi = 0; %选择测量矩阵，1为部分哈达玛矩阵，0为高斯随机矩阵
%-----------------------利用三角函数生成频域或DCT域离散信号--------------------------
n = 512;
t = [0: n-1];
f = cos(2*pi/256*t) + sin(2*pi/128*t); % 产生频域稀疏的信号
%-------------------------------信号降采样率-----------------------
n = length(f);
a = 0.2; % 取原信号的 a%
m = double(int32(a*n));
%--------------------------------------画信号图--------------------------------------
switch choice_transformcase 1ft = dct(f);disp('ft = dct(f)')case 0ft = fft(f);disp('ft = fft(f)')
end
disp(['信号稀疏度：',num2str(length(find((abs(ft))>0.1)))])
figure('name', 'A Tone Time and Frequency Plot');
subplot(2, 1, 1);
plot(f);
xlabel('Time (s)');
% ylabel('f(t)');
subplot(2, 1, 2);
switch choice_transformcase 1plot(ft)disp('plot(ft)')case 0plot(abs(ft));disp('plot(abs(ft))')
end
xlabel('Frequency (Hz)');
% ylabel('DCT(f(t))');
%% 产生感知矩阵和稀疏表示矩阵
%--------------------------利用感知矩阵生成测量值---------------------
switch choice_Phicase 1Phi = PartHadamardMtx(m,n); % 感知矩阵（测量矩阵） 部分哈达玛矩阵case 0Phi = sqrt(1/m) * randn(m,n); % 感知矩阵（测量矩阵） 高斯随机矩阵
end
% Phi = randn(m,n); %randn 生成标准正态分布的伪随机数（均值为0，方差为1）
% Phi = rand(m,n); % rand 生成均匀分布的伪随机数。分布在（0~1）之间
f2 = (Phi * f')'; % 通过感知矩阵获得测量值
% f2 = f(1:2:n);
switch choice_transformcase 1Psi = dct(eye(n,n)); %离散余弦变换正交基 代码亦可写为Psi = dctmtx(n);disp('Psi = dct(eye(n,n));')case 0Psi = inv(fft(eye(n,n))); % 傅里叶正变换，频域稀疏正交基（稀疏表示矩阵）disp('Psi = inv(fft(eye(n,n)));')
end
A = Phi * Psi; % 恢复矩阵 A = Phi * Psi
%% 重建信号
%---------------------使用CVX工具求解L1范数最小值-----------------
cvx_begin;variable x(n) complex;
% variable x(n) ;minimize( norm(x,1) );subject toA*x == f2' ;
cvx_end;
figure;subplot(2,1,2)
switch choice_transformcase 1plot(real(x));disp('plot(real(x))')case 0plot(abs(x));disp(' plot(abs(x))')
end
title('Using L1 Norm（Frequency Domain）');
% ylabel('DCT(f(t))'); xlabel('Frequency (Hz)');
switch choice_transformcase 1sig = dct(real(x));disp('sig = dct(real(x))')case 0sig = real(ifft(full(x)));disp(' sig = real(ifft(full(x)))')
end
subplot(2,1,1);
plot(f)
hold on;plot(sig);hold off
title('Using L1 Norm (Time Domain)');
% ylabel('f(t)'); xlabel('Time (s)');
legend('Original','Recovery')
%-----------------------------使用OMP算法重建-----------------------
% figure;plot(dct(theta));title(['K=',num2str(K)])switch choice_transformcase 1re(K) = norm(f'-(dct(theta)));case 0re(K) = norm(f'-real(ifft(full(theta))));end
end
theta = CS_OMP(f2,A,find(re==min(min(re))));
disp(['最佳稀疏度K=',num2str(find(re==min(min(re))))]);
% theta = CS_OMP(f2,A,10);
figure;subplot(2,1,2);
switch choice_transformcase 1plot(theta);disp('plot(theta)')case 0plot(abs(theta));disp('plot(abs(theta))')
end
title(['Using OMP(Frequence Domain) K=',num2str(find(re==min(min(re))))])
switch choice_transformcase 1sig2 = dct(theta);disp('sig2 = dct(theta)')case 0sig2 = real(ifft(full(theta)));disp('sig2 = real(ifft(full(theta)))')
end
subplot(2,1,1);plot(f);hold on;
plot(sig2)
hold off;
title(['Using OMP(Time Domain) K=',num2str(find(re==min(min(re))))]);
legend('Original','Recovery')
%%部分哈达玛矩阵：partHadamadMtx.m
function [ Phi ] = PartHadamardMtx( M,N )
%PartHadamardMtx Summary of this function goes here
% Generate part Hadamard matrix
% M -- RowNumber
% N -- ColumnNumber
% Phi -- The part Hadamard matrix
% 来源http://blog.csdn.net/jbb0523/article/details/44700735
%% parameter initialization
%Because the MATLAB function hadamard handles only the cases where n, n/12,
%or n/20 is a power of 2 L_t = max(M,N);%Maybe L_t does not meet requirement of function hadamard L_t1 = (12 - mod(L_t,12)) + L_t; L_t2 = (20 - mod(L_t,20)) + L_t; L_t3 = 2^ceil(log2(L_t)); L = min([L_t1,L_t2,L_t3]);%Get the minimum L
%% Generate part Hadamard matrix Phi = []; Phi_t = hadamard(L); RowIndex = randperm(L); Phi_t_r = Phi_t(RowIndex(1:M),:); ColIndex = randperm(L); Phi = Phi_t_r(:,ColIndex(1:N));
end 正交匹配追踪算法
function [ theta ] = CS_OMP( y,A,t )
% 实现压缩感知OMP算法
%CS_OMP Summary of this function goes here
%Version: 1.0 written by jbb0523 @2015-04-18
% Detailed explanation goes here
% y = Phi * x
% x = Psi * theta
% y = Phi*Psi * theta
% t 稀疏度
% 令 A = Phi*Psi, 则y=A*theta
% 现在已知y和A，求theta
% 来源：http://blog.csdn.net/jbb0523/article/details/45130793[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(列向量) At = zeros(M,t);%用来迭代过程中存储A被选择的列 Pos_theta = zeros(1,t);%用来迭代过程中存储A被选择的列序号 r_n = y;%初始化残差(residual)为y for ii=1:t%迭代t次，t为输入参数 product = A'*r_n;%传感矩阵A各列与残差的内积 [val,pos] = max(abs(product));%找到最大内积绝对值，即与残差最相关的列 At(:,ii) = A(:,pos);%存储这一列 Pos_theta(ii) = pos;%存储这一列的序号 A(:,pos) = zeros(M,1);%清零A的这一列，其实此行可以不要，因为它与残差正交 %y=At(:,1:ii)*theta，以下求theta的最小二乘解(Least Square) theta_ls = (At(:,1:ii)'*At(:,1:ii))^(-1)*At(:,1:ii)'*y;%最小二乘解 %At(:,1:ii)*theta_ls是y在At(:,1:ii)列空间上的正交投影 r_n = y - At(:,1:ii)*theta_ls;%更新残差 end theta(Pos_theta)=theta_ls;%恢复出的theta
end 三、 重点参考
https://blog.csdn.net/Di_Wong/article/details/81191211
https://blog.csdn.net/gongshouxiayin/article/details/122477776
![[声纹识别]基于MFCC的声纹识别算法](https://img-blog.csdnimg.cn/408da1548faa47ac8f1d922457bc268f.png)
