当前位置：首页 > news >正文

稀疏信号代码详解

news 2026/3/26 21:39:14

稀疏信号代码详解

一、代码整体架构

这段代码演示了压缩感知/稀疏信号恢复的核心思想：从一个低维观测向量b中恢复出一个高维稀疏向量u。代码包含四个主要部分：

数据生成：创建仿真数据
恢复算法：三种恢复方法（注释中展示了两种）
性能评估：计算恢复误差和稀疏度
结果可视化：绘制信号对比图

二、数据生成部分详解

m = 128; n = 256;
A = randn(m, n);
u = sprandn(n, 1, 0.1);
b = A * u;

2.1 维度设置 (`m = 128; n = 256;`)

n = 256：原始信号维度
- 我们想要恢复的信号的维度
- 代表"完整信息"的空间
m = 128：观测数据维度
- 实际测量到的数据的维度
- 代表"压缩后信息"的空间

关键数学关系：m < n（128 < 256）

这表示方程组 A*u = b 是欠定系统（未知数多于方程）
在传统线性代数中，这样的方程组有无穷多解
但我们有额外的先验知识：u是稀疏的

2.2 测量矩阵 (`A = randn(m, n);`)

A = randn(m, n);  % 生成m×n的随机高斯矩阵

维度：128行×256列
元素分布：每个元素独立服从标准正态分布 N(0,1)
为什么用随机高斯矩阵？
1. 限制等距性质：以高概率满足RIP条件
2. 通用性：对任何稀疏信号都能很好地工作
3. 信息保留：能最大程度保留稀疏信号的信息

数学意义：

A 是一个线性映射：ℝ²⁵⁶ → ℝ¹²⁸
它将256维空间压缩到128维空间
测量方程：b = A * u

2.3 稀疏信号生成 (`u = sprandn(n, 1, 0.1);`)

u = sprandn(n, 1, 0.1);  % 生成256×1的稀疏向量

sprandn：生成稀疏随机矩阵
(n, 1)：维度为256×1
0.1：稀疏度参数
- 表示大约10%的元素是非零的
- 实际非零元素数 ≈ 256 × 0.1 = 25.6 ≈ 26个

稀疏信号示例：

u = [0, 0, 1.23, 0, 0, 0, -0.56, 0, ..., 0, 2.11, 0, 0]^T↑     ↑               ↑                 ↑全部256个元素中，大约只有26个非零

2.4 观测数据生成 (`b = A * u;`)

b = A * u;  % 通过线性测量获得观测数据

维度：b 是128×1的向量
物理意义：我们无法直接观测256维的u，只能通过128个"传感器"（对应A的128行）获得压缩测量b
无噪声假设：这是理想情况，实际中通常有噪声：b = A*u + noise

数学模型：

b = A * u
[128×1] = [128×256] × [256×1]

三、恢复算法部分详解

3.1 方法1：LASSO回归（代码中使用的方法）

lambda = 0.1;  % 正则化参数
u_recovered = lasso(A, b, 'Lambda', lambda);

LASSO数学模型：

minimize ½‖A*u - b‖²₂ + λ‖u‖₁

第一项：数据拟合项（最小二乘）
- 使恢复的信号u的观测A*u尽可能接近实际观测b
第二项：L1正则化项
- ‖u‖₁ = Σ|uᵢ|，即所有元素绝对值之和
- 促使解变得稀疏（很多元素变为0）

参数`lambda = 0.1`的作用：

λ太大：过度惩罚，所有系数都趋向0，模型过于简单（欠拟合）
λ太小：惩罚不足，系数不稀疏，可能过拟合
λ=0.1：经验值，需要根据实际情况调整

`lasso`函数输出：

在MATLAB中，lasso函数可能返回：

单个向量：对应指定λ的解
矩阵：如果λ是向量，返回多个λ对应的解
结构体：包含系数、拟合信息等

潜在问题：如果u_recovered是矩阵，需要取第一列：

if size(u_recovered, 2) > 1u_recovered = u_recovered(:, 1);
end

3.2 方法2：正交匹配追踪（注释中）

% u_recovered = omp(A'*b, A'*A, 25);

OMP算法思想：

初始化：残差 r = b，支撑集 S = ∅
迭代：
- 找到与残差最相关的原子（A的列）
- 将该原子索引加入支撑集
- 在支撑集上做最小二乘
- 更新残差
停止：达到预设稀疏度（25）或残差足够小

参数：

25：估计的稀疏度（已知u大约有26个非零元素）
A'*b：相关性向量
A'*A：Gram矩阵

3.3 方法3：基追踪（注释中）

% cvx_begin
%     variable x(n)
%     minimize(norm(x, 1))
%     subject to
%         A*x == b
% cvx_end
% u_recovered = x;

数学模型（无噪声情况）：

minimize ‖u‖₁
subject to A*u = b

这是LASSO在无噪声时的特例（λ→0）
使用CVX凸优化工具箱求解

四、性能评估部分详解

recovery_error = norm(u_recovered - u) / norm(u);
sparsity = nnz(abs(u_recovered) > 1e-3) / n;

4.1 恢复误差计算

recovery_error = norm(u_recovered - u) / norm(u);

分子：norm(u_recovered - u)
- 计算恢复信号与原始信号的差异
- 使用2-范数（欧几里得距离）
分母：norm(u)
- 原始信号的2-范数
相对误差：无量纲，便于比较不同信号

数学公式：

相对误差 = ‖u_recovered - u‖₂ / ‖u‖₂

4.2 稀疏度计算

sparsity = nnz(abs(u_recovered) > 1e-3) / n;

nnz()：计算非零元素个数
abs(u_recovered) > 1e-3：判断元素是否"显著非零"
- 阈值1e-3忽略数值极小的元素（可能是数值误差）
除以n：得到稀疏度比例

示例计算：

如果u_recovered有30个绝对值大于0.001的元素
稀疏度 = 30 / 256 ≈ 0.117 = 11.7%

4.3 结果输出

fprintf('恢复相对误差: %.4e\n', recovery_error);
fprintf('恢复信号的稀疏度: %.2f%%\n', sparsity*100);

%.4e：科学计数法，保留4位小数
%.2f%%：百分比形式，保留2位小数
理想情况：误差接近0，稀疏度接近10%

五、可视化部分详解

figure;
subplot(2,1,1);
stem(u, 'b.', 'MarkerSize', 10);
title('原始稀疏信号 u');
xlabel('索引'); ylabel('幅值');subplot(2,1,2);
stem(u_recovered, 'r.', 'MarkerSize', 10);
title('恢复的信号 u_{recovered}');
xlabel('索引'); ylabel('幅值');

5.1 图形窗口

figure;  % 创建新的图形窗口

5.2 第一个子图：原始信号

subplot(2,1,1);  % 2行1列，第1个位置
stem(u, 'b.', 'MarkerSize', 10);  % 蓝色点状杆图
title('原始稀疏信号 u');
xlabel('索引'); ylabel('幅值');

stem图：适合显示离散信号
蓝色('b')：表示原始信号
点标记('.')：小点表示每个采样点
MarkerSize：标记大小

预期效果：

大部分位置高度为0（稀疏性）
少数位置有显著非零值
非零值随机分布

5.3 第二个子图：恢复信号

subplot(2,1,2);  % 2行1列，第2个位置
stem(u_recovered, 'r.', 'MarkerSize', 10);  % 红色点状杆图
title('恢复的信号 u_{recovered}');
xlabel('索引'); ylabel('幅值');

红色('r')：表示恢复的信号
理想情况应与第一个子图基本一致

六、潜在问题与改进建议

6.1 代码中的问题

lasso输出格式：可能需要处理矩阵输出
工具箱依赖：需要Statistics and Machine Learning Toolbox
缺少网格和标签：可视化可读性可提升

6.2 改进版本

% 修正lasso输出
[B, FitInfo] = lasso(A, b, 'Lambda', lambda);
if size(B, 2) > 1u_recovered = B(:, 1);
elseu_recovered = B;
end% 改进可视化
figure('Position', [100, 100, 800, 600]);subplot(2,1,1);
stem(u, 'b.', 'MarkerSize', 10, 'LineWidth', 1.5);
title('原始稀疏信号 u (10%稀疏度)');
xlabel('索引'); ylabel('幅值');
grid on;
xlim([1, n]);subplot(2,1,2);
stem(u_recovered, 'r.', 'MarkerSize', 10, 'LineWidth', 1.5);
title(sprintf('恢复信号 (误差: %.2e, 稀疏度: %.1f%%)', ...recovery_error, sparsity*100));
xlabel('索引'); ylabel('幅值');
grid on;
xlim([1, n]);% 添加总标题
sgtitle(sprintf('稀疏信号恢复: m=%d, n=%d, λ=%.2f', m, n, lambda));

6.3 扩展功能

% 计算更多指标
support_original = find(abs(u) > 1e-3);
support_recovered = find(abs(u_recovered) > 1e-3);
intersection_rate = length(intersect(support_original, support_recovered)) / ...length(support_original) * 100;fprintf('支撑集重合率: %.1f%%\n', intersection_rate);
fprintf('原始信号非零数: %d\n', length(support_original));
fprintf('恢复信号非零数: %d\n', length(support_recovered));