稀疏信号代码详解（MATLAB版）

MATLAB代码是压缩感知/稀疏信号恢复的经典仿真示例。以下机型详细解释其数学意义和应用场景：

代码详解

1. 参数设置

m = 128; n = 256;

n = 256: 原始信号 u 的维度（高维空间）
m = 128: 观测数据 b 的维度（低维空间）
关键: m < n，观测数量远小于信号维度，这是欠定系统

2. 测量矩阵

A = randn(m, n);

创建一个 128×256 的随机高斯矩阵
为什么用随机矩阵？
- 满足限制等距性质（RIP） 的概率很高
- 可以高效地压缩稀疏信号的信息
- 在压缩感知理论中，随机测量是最优的

3. 稀疏信号

u = sprandn(n, 1, 0.1);

生成 256×1 的稀疏向量
0.1 表示稀疏度：大约10%的元素非零
实际非零元素数 ≈ 256 × 0.1 = 25.6 个
稀疏性是能够恢复的关键前提

4. 观测数据

b = A * u;

通过线性测量获得压缩后的数据
b 是 128×1 的向量
无噪声的理想情况：b = A * u

数学问题描述

已知：A ∈ ℝ^(m×n), b ∈ ℝ^m, 且 m < n
求解：u ∈ ℝ^n 满足 A u = b

核心矛盾：

方程 A u = b 有无穷多解（因为A是"矮胖"矩阵）
但u具有稀疏性先验知识

解决方案：稀疏恢复

理论依据：压缩感知理论

如果满足：

信号u是稀疏的（或可压缩的）
测量矩阵A满足RIP性质
观测数量 m ≥ O(k log(n/k))，其中k是稀疏度

则可以通过ℓ₁最小化精确恢复原始信号。

恢复算法

通常使用以下方法求解：

1. 基追踪（Basis Pursuit）

% 无噪声情况下的精确恢复
min ||u||₁  s.t.  A*u = b

2. LASSO/基追踪去噪

% 有噪声情况下的稳健恢复
min ½||A*u - b||₂² + λ||u||₁

3. 贪婪算法

正交匹配追踪（OMP）
压缩采样匹配追踪（CoSaMP）

完整实验示例

% 生成数据 
m = 128; n = 256;
A = randn(m, n);
u = sprandn(n, 1, 0.1);
b = A * u;

% 方法1：使用LASSO恢复（需要优化工具箱）
lambda = 0.1;  % 正则化参数
u_recovered = lasso(A, b, 'Lambda', lambda);

% 方法2：使用OMP（需要SPAMS或其他工具箱）
% u_recovered = omp(A'*b, A'*A, 25);  % 25是稀疏度估计

% 方法3：使用CVX工具箱求解基追踪
% cvx_begin
%     variable x(n)
%     minimize(norm(x, 1))
%     subject to
%         A*x == b
% cvx_end
% u_recovered = x;

% 评估恢复效果
recovery_error = norm(u_recovered - u) / norm(u);
sparsity = nnz(abs(u_recovered) > 1e-3) / n;

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

% 可视化
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('幅值');

实际应用场景

1. 医学成像（MRI加速）

原始信号：高分辨率图像（稀疏在小波域）
测量矩阵：部分傅里叶采样
观测数据：部分k空间数据
恢复结果：完整图像

2. 单像素相机

原始信号：场景图像
测量矩阵：随机掩模模式
观测数据：单个传感器的测量序列
恢复结果：完整图像

3. 无线传感网络

原始信号：多个传感器的读数
测量矩阵：随机投影
观测数据：压缩后的聚合数据
恢复结果：所有传感器的原始数据

4. 故障诊断

原始信号：机器多个位置的振动信号
测量矩阵：随机采样
观测数据：少数几个传感器的读数
恢复结果：完整振动图谱，用于故障定位

关键理论要点

1. 稀疏表示

信号u本身可能不稀疏，但在某个变换域（如傅里叶、小波、DCT）是稀疏的：

u = Ψ * α  % α是稀疏系数，Ψ是稀疏基
b = A * Ψ * α = Θ * α  % Θ = AΨ是感知矩阵

2. 受限等距性质（RIP）

对于任意k-稀疏向量v，矩阵A满足：

(1 - δ_k) ||v||₂² ≤ ||A v||₂² ≤ (1 + δ_k) ||v||₂²

其中δ_k是RIP常数，0 < δ_k < 1

3. 恢复保证

当A满足RIP且δ_2k < √2 - 1时，通过ℓ₁最小化可以：

精确恢复无噪声情况下的信号
稳定恢复有噪声情况下的信号

总结

这段代码虽然简短，但包含了压缩感知的核心思想：

维数灾难的突破：从少量观测恢复高维信号
稀疏性的价值：稀疏性作为有效的先验信息
随机性的力量：随机测量是信息保留的高效方式

这个框架不仅在理论上优美，在实际应用中也产生了深远影响，特别是在数据采集、传输和存储效率方面带来了革命性的改进。

稀疏信号代码详解

一、代码整体架构

这段代码演示了压缩感知/稀疏信号恢复的核心思想：从一个低维观测向量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) > 1
    u_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) > 1
    u_recovered = B(:, 1);
else
    u_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));

七、理论背景总结

压缩感知核心思想：

欠定系统可解：当m < n时，A*u = b有无穷多解
稀疏性先验：真实信号u是稀疏的（大部分元素为0）
L1最小化：在所有满足A*u = b的解中，寻找L1范数最小的解
恢复保证：如果A满足RIP条件且u足够稀疏，L1最小化的解等于L0最小化的解

关键公式：

观测模型：b = A*u（无噪声）
LASSO：min ½‖A*u-b‖² + λ‖u‖₁
基追踪：min ‖u‖₁ s.t. A*u = b
恢复误差：‖u_rec - u‖₂ / ‖u‖₂

这段代码简短，完整展示了压缩感知的基本流程，是理解现代信号处理中稀疏恢复技术的优秀起点。

posted @ 2026-02-02 23:48 kkman2000 阅读(0) 评论(0) 收藏举报

刷新页面返回顶部

kkman2000