压缩感知（Compressive Sensing，CS）理论。搬运自《阵列信号处理与MATLAB实现》
信号的稀疏性是压缩感知的重要前提和理论基础，信号的稀疏性由如下定义：
- 定义 1：信号的稀疏性为信号中的非零元素数目较少。
- 定义 2：信号在某个变换域下近似稀疏，即为可压缩信号。从理论上讲对于任何信号都具有可压缩性，只要能找到其相应的稀疏表示空间，就可以有效地进行压缩采样。
- 定义 3：矩阵奇异值的稀疏性：矩阵奇异值中非零元素的个数（矩阵的秩）相对较少；也称为矩阵的低秩性，即矩阵的秩相对于矩阵的行数或列数而言很小。
1. 压缩感知的理论框架
1.1 信号的稀疏表示
$$ \newcommand{\vect}[1]{\boldsymbol{#1}} $$定义矢量$\vect{x} = [x_1, x_2,\dots,x_N]^T$的$L_p$范数。
$$ ||x||_p = (\sum_{i=1}^N|x_i|^p)^{\frac{1}{p}} $$$L_0$范数：向量中非零元素的个数
$$ |x|_0 = \sharp\{i \mid x_i \neq 0\} $$其中 $\sharp$ 表示集合中元素的个数（计数）
===
对于信号$x \in R^N$，其在标准正交基$\Psi$小的表示系数矢量为：
$$ \vect{a} = \Psi^T x $$根据$L_p$范数的定义，若$\vect{a}$满足
$$ ||\vect{a}||_p \le K $$对于实数$0< p < 2 $和$K>2$同时成立，则称$x$在变换域$\Psi$下是稀疏的。特别的，当$p=0$时，称信号$x$在变换域$\Psi$下是$K$稀疏的（不多于$K$个非零元素）。
首先考虑一般的信号重构问题，信号$x$在时域本身就是稀疏或可压缩的，即上述的基变换$\Psi$为$\delta(t)$函数（不做任何变换），给定一个投影测量矩阵$\Phi \in \mathbb{R}^{M\times N} \quad (M \ll N )$，则信号$x$在该测量矩阵$\Phi$下的线性投影测量值为：
$$ \begin{equation} y = \Phi x \tag{1} \end{equation} $$考虑由投影测量数据$y$来重构信号$x$。 由于$y$的维数$M$远小于$x$的维数$N$，方程(1)是欠定方程，有无穷多个解，直接通过解方程的方法无法重构原始信号。
压缩感知理论证明：如果原始信号$x$本身在时域是$K-$稀疏的或可压缩的，并且$y$与$\Phi$满足一定的条件， 信号$x$可以由投影测量值$y$通过求解以下最小$L_0$范数问题以极高概率得到其精确的重构。
自然信号在时域内几乎都是不稀疏的，但可能在某个变换基域稀疏。设自然信号$x$在基变换$\Psi$下具有稀疏性或可压缩性，即
$$ \vect{x}=\Psi\vect{a} $$$\vect{a}$为信号$x$在变换域$\Psi$下的$K$稀疏表示系数。于是信号$x$在投影测量矩阵$\Phi$下线性测量可以表示为：
$$ \begin{equation} y = \Phi x = \Phi\Psi a = \tilde{\Phi} a \tag{2} \end{equation} $$其中，$\tilde{\Phi} = \Phi\Psi$为$M\times N$维矩阵，称为感知矩阵。
$y$可以看作稀疏信号$a$关于$\Phi$的线性测量。由于变换域$\Psi$固定，要使$\tilde{\Phi} = \Phi\Psi$满足RIP条件。
受限等距属性（RIP）
什么样的测量矩阵 $\Phi$ 能保证这种恢复是唯一的、稳定的？RIP 条件就是用来刻画测量矩阵“好坏”的标准。
定义：设 $\Phi$ 是一个 $M \times N$ 的矩阵。如果对于所有的 $K$-稀疏向量 $x \in \mathbb{R}^N$（即 $|x|_0 \leq K$），都存在一个常数 $\delta_K \in (0, 1)$，使得以下不等式成立：
$$ (1 - \delta_K) |x|_2^2 \leq ||\Phi x||_2^2 \leq (1 + \delta_K) |x|_2^2 $$则称矩阵 $\Phi$ 满足$K$ 阶受限等距属性（RIP of order $K$）。其中，满足上述条件的最小常数 $\delta_K$ 被称为RIP 常数（RIP Constant）。
定义解析：
- 等距性（Isometry）： 如果 $\delta_k = 0$，则 $||\Phi x||_2 = ||x||_2$。这意味着变换 $\Phi$ 保持了向量的欧几里得长度（能量）不变。这通常发生在正交矩阵中。
- 限制性（Restricted）： 这个性质不需要对所有向量 $x$ 成立，只需要对稀疏向量成立。这是压缩感知的关键，因为我们要处理的信号本身就是稀疏的。
- 常数 $\delta_k$： $\delta_k$ 越小，说明矩阵 $\Phi$ 对稀疏向量的长度保持得越好。通常要求 $\delta_k$ 足够小（例如 $\delta_{2K} < \sqrt{2}-1 \approx 0.414$ 或更宽松的条件），才能保证恢复算法的成功。
1.2 测量矩阵的设计
在压缩感知理论中，得到信号的稀疏表示后，设计一个测量矩阵$\Phi$，使得在该测量矩阵上压缩投影得到的M个测量值能够保留原始信号的绝大部分信息，使原始信号的信息损失最小，从而保证从这些少量测量值中能够精确重构出长度为$N (M \ll N)$的原始信号。
基本准则：
- 非相干性
- 等距约束性
1.3 信号重构算法
信号重构：从$M$个测量值中重构出长度为$N$ （$M \ll N$）的稀疏信号。
若信号是稀疏的或可压缩的，且感知矩阵$\tilde{\Phi}$满足RIP等稀疏重构条件，则该信号可以很高的概率被稀疏重构出来。
信号重构问题为以下的最小$L_0$范数问题求解：
$$ \hat{\vect{a}} = \arg \min||\vect{a}||_0 \quad \text{s.t. } \quad \tilde{\Phi}\vect{a} = \vect{y} \tag{3} $$含义： 在所有能满足测量方程 $y = \Psi \alpha$ 的解中，寻找非零元素个数最少的那个解
求解上述问题需要列出$a$中所有非零项位置的 $\binom{N}{K}$ 种可能的组合才能得到最优解。 在$N$很大时无法有效实现。
近似等价的信号重建算法：松弛算法、贪婪算法、非凸方法
1.3.1 基追踪
采用$L_1$范数代替$L_0$范数可以凸优化求解
$$ \hat{\vect{a}} = \arg\min||\vect{a}||_1 \quad\text{s.t. }\quad \tilde{\Phi}\vect{a} = \vect{y} \tag{4} $$问题$(3)$和问题$(3)$在满足一定条件时是等价的，因此信号重构问题可以转化为一个线性规划问题加以求解，这种方法也称为基追踪（Basis Pursuit, BP）。
考虑噪声的存在，上述问题可以转换为如下最小$L_1$范数问题（基追踪降噪，Basis Pursuit De-Noising, BPDN）：
$$ \hat{\vect{a}} = \arg\min||\vect{a}||_1 \quad\text{s.t. }\quad ||\tilde{\Phi}\vect{a} - \vect{y}||_2 \le \sigma \tag{5} $$$\sigma$代表噪声的一个可能的标准差。
含义：在所有能够把观测误差控制在 $\epsilon$ 范围内的解（满足数据保真度）中，寻找一个 $L_1$ 范数最小的解（最稀疏的解）。
2. 基于压缩感知理论的DOA估计算法
假设$K$个远场窄带信号入射到$M$个天线的均匀线阵上，第$k$个信号入射角度为$\theta_k$。$t$时刻阵列接收的单快拍数据矢量可以表示为：
$$ \vect{x}(t) = \vect{A}(\theta)\vect{s}(t) + \vect{n}(t) $$其中$\vect{s}(t)= [ s_1 \quad s_2 \quad \dots \quad s_K]^T $表示$K$个信源矢量。 $\vect{A}(\theta) = [\vect{a}(\theta_1) \quad \vect{a}(\theta_2) \quad \dots \quad \vect{a}(\theta_K)]$为方向矩阵。 $\vect{n}(t)$为阵列接收噪声。
下面，对阵列流型$\vect{A}$进行扩展，形成完备冗余字典$\vect{G}$，使其包含所有可能的方位角：
$$ \vect{G}(\theta) = [\vect{a}(\theta_1) \quad \vect{a}(\theta_2) \quad \dots \quad \vect{a}(\theta_Q)] $$其中$\vect{a}(\theta_i)$为原子，Q为方位角$\theta_k$的等网格划分数目。利用完备冗余字典$\vect{G}$将上式转化为稀疏表示问题：
$$ \vect{x}= \vect{G}\vect{\delta} + \vect{n} $$稀疏信号$\delta$为Q维度稀疏向量，非零元素个数为$K$，即$\vect{x}$基于字典$\vect{G}$是$K$-稀疏的。非零元素的位置的对应角度代表了入射角$\theta$的值，非零元素幅值即信号在采样时刻的幅度。
对$\vect{\delta}$的求解：
$$ \hat{\vect{\delta}} = \arg\min||\delta||_{0} \quad\text{s.t. }\quad \vect{x} = \vect{G}\vect{\delta} + \vect{n} $$考虑噪声的存在，上述问题可以转换为如下最小$L_1$范数问题：
$$ \hat{\vect{\delta}} = \arg\min||\vect{\delta}||_1 \quad\text{s.t. }\quad ||\vect{x} - \vect{G}\vect{\delta}||_2^2 \le \gamma^2 $$基于压缩感知DOA估计算法MATLAB
classdef CSEstimator < handle
% CSEstimator 基于压缩感知的 DOA 估计器
% 使用稀疏信号恢复方法实现角度估计
%
% 使用示例:
% estimator = CSEstimator('NumAntennas', 8, 'AngleRange', -90:0.5:90);
% est_aoa = step(estimator, signal_data);
% plotSpectrum(estimator);
properties
% 系统参数
NumAntennas = 8; % 天线数量 (M)
OperatingFrequency = 5.8e9; % 中心频率 (Hz)
ElementSpacing = 0.0258; % 天线间距 (米, 默认为半波长)
SensorArray; % 天线阵列对象 (phased.ULA)
% 算法参数
NumSignals = 1; % 估计的路径数量
NumSignalsSource = 'Auto'; % 信号数估计方式: 'Auto' 或 'Property'
NumSignalsMethod = 'MDL'; % 自动估计方法: 'MDL' 或 'AIC'
% 搜索网格参数
AngleRange = -90:1:90; % 角度搜索范围 (度)
% 压缩感知参数
RegularizationParam = 0.1; % 正则化参数 lambda
MaxIterations = 1000; % 最大迭代次数
Tolerance = 1e-6; % 收敛阈值
Algorithm = 'OMP'; % 算法选择: 'OMP', 'LASSO', 'ISTA'
end
properties (Access = private)
% 内部缓存
DictionaryMatrix % 过完备字典矩阵
LastSpectrum % 最近一次计算的谱
LastSparseCoeff % 最近一次稀疏系数
end
methods
function obj = CSEstimator(varargin)
% 构造函数: 解析输入参数 (Name-Value pairs)
if nargin > 0
for i = 1:2:nargin
if isprop(obj, varargin{i})
obj.(varargin{i}) = varargin{i+1};
else
error('Property %s does not exist', varargin{i});
end
end
end
% 自动计算半波长间距
c = physconst('LightSpeed');
lambda = c / obj.OperatingFrequency;
if isempty(obj.SensorArray)
obj.ElementSpacing = lambda / 2;
else
obj.ElementSpacing = obj.SensorArray.ElementSpacing;
obj.NumAntennas = obj.SensorArray.NumElements;
end
% 构建字典矩阵
obj.buildDictionary();
end
function D = getNumSignals(obj, eigenvals, K, fb)
% 获取信号子空间维度
if strcmp(obj.NumSignalsSource, 'Auto')
if strcmp(obj.NumSignalsMethod, 'AIC')
D = phased.internal.aictest(eigenvals, K, fb);
else
D = phased.internal.mdltest(eigenvals, K, fb);
end
else
D = obj.NumSignals;
end
end
function [est_aoa, spectrum] = step(obj, signal_data)
% STEP 执行压缩感知 DOA 估计
% 输入 signal_data: [M x N] (天线数, 快拍数)
% 输出 est_aoa: 估计的角度 (度)
% 输出 spectrum: 稀疏恢复系数 (用于绘图)
[M, N] = size(signal_data);
% 确保天线数匹配
if M ~= obj.NumAntennas
error('输入数据的天线数 (%d) 与设置的天线数 (%d) 不匹配', M, obj.NumAntennas);
end
% 获取字典矩阵
A = obj.DictionaryMatrix;
% 对每个快拍进行稀疏恢复，然后取平均
n_angles = length(obj.AngleRange);
sparse_coeff = zeros(n_angles, 1);
switch obj.Algorithm
case 'OMP'
% 估计信号数量用于OMP
R = signal_data * signal_data' / N;
[~, eig_vals_mat] = eig(R);
eig_vals = sort(real(diag(eig_vals_mat)), 'descend');
D = obj.getNumSignals(eig_vals, N, false);
D = max(D, 1); % 至少估计一个信号
for n = 1:N
y = signal_data(:, n);
x = obj.omp(A, y, D);
sparse_coeff = sparse_coeff + abs(x).^2;
end
case 'LASSO'
for n = 1:N
y = signal_data(:, n);
x = obj.lasso_ista(A, y);
sparse_coeff = sparse_coeff + abs(x).^2;
end
case 'ISTA'
for n = 1:N
y = signal_data(:, n);
x = obj.lasso_ista(A, y);
sparse_coeff = sparse_coeff + abs(x).^2;
end
otherwise
error('未知算法: %s', obj.Algorithm);
end
sparse_coeff = sparse_coeff / N;
obj.LastSparseCoeff = sparse_coeff;
% 转换为谱 (归一化)
spectrum = sparse_coeff / max(sparse_coeff);
obj.LastSpectrum = spectrum;
% 寻找峰值
[peaks, locs] = findpeaks(spectrum, 'SortStr', 'descend');
if isempty(locs)
% 如果没有找到峰值，取最大值位置
[~, locs] = max(spectrum);
end
% 返回估计的角度
est_aoa = -obj.AngleRange(locs(1:min(obj.NumSignals, length(locs))));
spectrum = fliplr(spectrum.');
end
function plotSpectrum(obj)
% PLOTSPECTRUM 绘制空间谱
if isempty(obj.LastSpectrum)
error('请先调用 step() 方法计算谱。');
end
figure;
P_log = 10 * log10(obj.LastSpectrum + eps);
plot(obj.AngleRange, P_log, 'LineWidth', 1.5);
xlabel('角度 (度)');
ylabel('空间谱 (dB)');
title('压缩感知 DOA 估计谱');
grid on;
xlim([min(obj.AngleRange), max(obj.AngleRange)]);
end
end
methods (Access = private)
function buildDictionary(obj)
% 构建过完备导向矢量字典
n_angles = length(obj.AngleRange);
M = obj.NumAntennas;
d = obj.ElementSpacing;
lambda = physconst('LightSpeed') / obj.OperatingFrequency;
% 天线索引
idx_m = (0:M-1)';
% 构建字典矩阵 [M x n_angles]
obj.DictionaryMatrix = zeros(M, n_angles);
for i = 1:n_angles
theta = obj.AngleRange(i);
phi = -2 * pi * d / lambda * sind(theta);
obj.DictionaryMatrix(:, i) = exp(1j * idx_m * phi);
end
end
function x = omp(obj, A, y, K)
% OMP (Orthogonal Matching Pursuit) 算法
% A: 字典矩阵 [M x N]
% y: 观测向量 [M x 1]
% K: 稀疏度 (信号数)
[~, N] = size(A);
x = zeros(N, 1);
residual = y;
support = [];
for k = 1:K
% 计算相关性
correlations = abs(A' * residual);
% 排除已选择的原子
correlations(support) = 0;
% 选择最大相关的原子
[~, idx] = max(correlations);
support = [support; idx];
% 最小二乘求解
A_s = A(:, support);
x_s = A_s \ y;
% 更新残差
residual = y - A_s * x_s;
% 检查收敛
if norm(residual) < obj.Tolerance
break;
end
end
% 填充稀疏系数
x(support) = x_s;
end
function x = lasso_ista(obj, A, y)
% ISTA (Iterative Shrinkage-Thresholding Algorithm) 用于 LASSO
% min_x ||y - Ax||_2^2 + lambda * ||x||_1
[~, N] = size(A);
lambda = obj.RegularizationParam;
% 计算步长 (Lipschitz 常数的倒数)
L = max(eig(A' * A));
step_size = 1 / L;
% 初始化
x = zeros(N, 1);
for iter = 1:obj.MaxIterations
x_old = x;
% 梯度步
gradient = A' * (A * x - y);
z = x - step_size * gradient;
% 软阈值 (复数版本)
x = obj.soft_threshold_complex(z, lambda * step_size);
% 检查收敛
if norm(x - x_old) / (norm(x_old) + eps) < obj.Tolerance
break;
end
end
end
function y = soft_threshold_complex(~, x, thresh)
% 复数软阈值函数
magnitude = abs(x);
y = max(magnitude - thresh, 0) .* exp(1j * angle(x));
end
end
end
