简介：这份Matlab代码资源面向从事高维数据分析、信号处理与机器学习的研究者和工程师，聚焦于用凸优化方法实现张量分解，从而完成数据恢复与去噪任务。张量作为矩阵向高维的推广，能有效刻画图像、视频、多模态数据的内在结构，而凸优化则保证低秩近似求解的稳定性与全局最优性，适合具备一定Matlab基础、希望深入张量处理与优化技术的读者。资源包共48个文件，以45个m脚本为核心，辅以README.md说明文档及git配置项，整体约55KB，涵盖张量分解、ADM求解、软阈值、矩阵补全、去噪与补全实验及绘图对比等模块，目录组织清晰。目前已有738人学习下载。读者可借助其中的算法函数、示例脚本与测试用例，快速理解CP、Tucker等分解思路，并依据自身数据类型调整迭代次数与阈值，完成从数据预处理到去噪效果可视化的完整流程。
1. 从一段被噪声淹没的信号说起：这套张量分解去噪代码能帮你做什么
如果你手头有一段多通道采集的信号——比如阵列传感器输出、高光谱波段数据、或者某个实验平台同时记录的多个通道——大概率会遇到同一个问题：单通道看还行，多通道拼在一起就全是噪声，传统的一维滤波按通道逐个处理，处理完通道之间的相关性也被切断了。这套基于凸优化的张量分解 Matlab 代码，解决的正是这个场景：它把多通道数据组织成高阶张量，用张量分解的方式把低秩的信号子空间和稀疏的噪声分离开，而不是逐通道去卷滤波核。代码是完整的 Matlab 工程，包含分解主流程、凸优化求解和去噪演示脚本，适合做信号处理、阵列数据处理、图像去噪方向的研究生和工程师直接拿来跑，也适合在此基础上做二次开发。它不依赖深度学习工具箱，纯数值计算路线，对没有 GPU 的环境很友好。
2. 技术选型与整体架构：为什么用张量分解而不是逐通道滤波
2.1 张量分解去噪的基本思路
先把问题说清楚。假设你采集到的数据是一个三阶张量 $\mathcal{X} \in \mathbb{R}^{I_1 \times I_2 \times I_3}$，其中两个维度是空间或通道维度，第三个维度是时间或样本维度。真实的干净信号 $\mathcal{L}$ 通常具有低秩结构——因为多通道之间的信号是高度相关的，而噪声 $\mathcal{E}$ 往往是稀疏的或者满秩的。于是去噪问题可以写成：
$$\min_{\mathcal{L}, \mathcal{E}} \ \text{rank}(\mathcal{L}) + \lambda |\mathcal{E}|_0 \quad \text{s.t.} \quad \mathcal{X} = \mathcal{L} + \mathcal{E}$$
这个问题的难点在于 rank 和 L0 范数都是非凸的，直接求解是 NP 难的。工程上的常见做法是松弛成凸问题：用张量的核范数（tensor nuclear norm）替代 rank，用 L1 范数替代 L0，得到：
$$\min_{\mathcal{L}} \ |\mathcal{L}|_* + \lambda |\mathcal{X} - \mathcal{L}|_1$$
这就是这套代码的核心优化目标。它属于凸优化范畴，可以用 ADMM（交替方向乘子法）或者近端梯度法稳定求解。相比逐通道小波阈值去噪，这种做法的优势在于：分解过程天然保留了通道间的耦合关系，去噪后的张量在结构上更一致，不会出现某个通道被过度平滑、另一个通道还残留噪声的情况。
2.2 代码的模块划分
拿到源码包后，目录结构大致是这样的（不同版本可能略有差异，但核心模块是固定的）：
tensor-denoise/ ├── main_denoise.m % 主入口，加载数据、调参、输出去噪结果 ├── core/ │ ├── tensor_nuclear_norm.m % 张量核范数计算（基于 t-SVD） │ ├── tsvd.m % 张量奇异值分解实现 │ ├── soft_threshold.m % 软阈值算子，L1 近端算子 │ └── admm_solver.m % ADMM 主循环 ├── utils/ │ ├── tensor2mat.m % 张量展开/折叠工具 │ ├── psnr_ssim.m % 去噪评价指标 │ └── gen_test_data.m % 生成带噪声的测试张量 └── data/ └── sample.mat % 示例数据这个划分很清晰：core/放算法核心，utils/放工具函数，main_denoise.m是唯一需要你手动改的入口。我一般拿到这种结构的项目，第一件事就是打开main_denoise.m看它默认加载什么数据、参数怎么设，然后直接跑一遍看输出。
2.3 凸优化求解器的选择
代码里用的是 ADMM，而不是 CVX 这类建模工具箱。原因很实际：CVX 虽然写起来简洁，但每次求解都要重新解析模型，对于需要反复迭代调参的去噪任务来说太慢，而且 CVX 对三阶以上张量的支持并不好。ADMM 把原问题拆成两个子问题交替求解，每个子问题都有闭式解，迭代几十次就能收敛，速度上完全可接受。
ADMM 的增广拉格朗日形式为：
$$\mathcal{L}\rho(\mathcal{L}, \mathcal{E}, \mathcal{Y}) = |\mathcal{L}|* + \lambda |\mathcal{E}|_1 + \langle \mathcal{Y}, \mathcal{X} - \mathcal{L} - \mathcal{E} \rangle + \frac{\rho}{2} |\mathcal{X} - \mathcal{L} - \mathcal{E}|_F^2$$
其中 $\mathcal{Y}$ 是拉格朗日乘子张量，$\rho$ 是惩罚参数。迭代更新规则在admm_solver.m里实现，下面会逐段拆。
3. 核心模块拆解：从 t-SVD 到 ADMM 迭代的完整链路
3.1 张量奇异值分解（t-SVD）的实现
张量核范数的计算依赖 t-SVD。和矩阵 SVD 不同，t-SVD 是在傅里叶域沿第三维做矩阵 SVD，然后逆变换回来。tsvd.m的核心逻辑如下：
function [U, S, V] = tsvd(X) % 输入 X: n1 x n2 x n3 张量 % 输出 U, S, V: t-SVD 分解结果 [n1, n2, n3] = size(X); % 沿第三维做 FFT，转到傅里叶域 Xf = fft(X, [], 3); U = zeros(n1, n1, n3); S = zeros(n1, n2, n3); V = zeros(n2, n2, n3); for i = 1:n3 % 对每个频率切片做矩阵 SVD [Uf, Sf, Vf] = svd(Xf(:,:,i), 'econ'); U(:,:,i) = Uf; S(:,:,i) = Sf; V(:,:,i) = Vf; end % 逆 FFT 回到原域 U = ifft(U, [], 3); S = ifft(S, [], 3); V = ifft(V, [], 3); end这里有几个关键点。第一，fft(X, [], 3)是沿第三维做变换，这个维度通常对应时间或样本，变换后每个频率切片是一个独立的矩阵。第二，svd(..., 'econ')用经济型分解，因为张量切片往往不是方阵，经济型能避免产生大量零奇异值。第三，逆变换后得到的 U、S、V 是复数张量，但实际使用时通常只取实部或者模值，因为原始数据是实数的。
注意：如果你的 Matlab 版本较老（R2016 之前），
fft对三维数组的支持可能有问题，需要手动 reshape 成二维再处理。常见做法是用reshape把n1 x n2 x n3变成n1 x (n2*n3)，做完 FFT 再变回来。
3.2 软阈值算子与核范数近端算子
ADMM 的每个子问题都涉及近端算子。对于 L1 范数，近端算子就是软阈值：
function Y = soft_threshold(X, tau) % 软阈值算子: Y = sign(X) * max(|X| - tau, 0) % X: 输入张量/矩阵 % tau: 阈值参数，对应 ADMM 中的 lambda/rho Y = sign(X) .* max(abs(X) - tau, 0); end这个函数看起来简单，但tau的取值直接决定去噪强度。tau太大，信号细节被一起抹掉；tau太小，噪声残留明显。代码里默认tau = lambda / rho，其中lambda是用户在主脚本里设的稀疏惩罚系数，rho是 ADMM 惩罚参数。我一般会先固定rho = 1，然后在一个对数尺度上扫lambda（比如 0.001 到 1），看哪个值对应的 PSNR 最高。
对于核范数，近端算子是奇异值软阈值（Singular Value Thresholding, SVT）：
function L = svt_tensor(X, tau) % 张量奇异值软阈值: 对奇异值做软阈值后重构 [U, S, V] = tsvd(X); S_thresh = soft_threshold(S, tau); L = tprod(U, tprod(S_thresh, ttranspose(V))); end这里用到了tprod（张量乘积）和ttranspose（张量转置），它们是 t-SVD 体系下的基本运算。tprod的实现本质是在傅里叶域做矩阵乘法再逆变换，和 t-SVD 的推导保持一致。
3.3 ADMM 主循环的逐段解析
admm_solver.m是整个项目的心脏。它的输入是观测张量X、参数lambda、rho、最大迭代次数max_iter和收敛阈值tol，输出是去噪后的张量L和噪声张量E。
function [L, E, obj_history] = admm_solver(X, lambda, rho, max_iter, tol) % ADMM 求解: min ||L||_* + lambda*||E||_1 s.t. X = L + E [n1, n2, n3] = size(X); % 初始化 L = zeros(n1, n2, n3); E = zeros(n1, n2, n3); Y = zeros(n1, n2, n3); % 拉格朗日乘子 obj_history = zeros(max_iter, 1); for k = 1:max_iter % --- L 子问题: 核范数近端算子 --- % L = argmin ||L||_* + (rho/2)||X - L - E + Y/rho||_F^2 temp = X - E + Y / rho; L = svt_tensor(temp, 1 / rho); % --- E 子问题: L1 近端算子 --- % E = argmin lambda*||E||_1 + (rho/2)||X - L - E + Y/rho||_F^2 temp = X - L + Y / rho; E = soft_threshold(temp, lambda / rho); % --- 乘子更新 --- Y = Y + rho * (X - L - E); % --- 收敛判断 --- primal_res = norm(X(:) - L(:) - E(:)); dual_res = rho * norm(L(:) - L_old(:)); % 需要保存上一步的 L obj_history(k) = tensor_nuclear_norm(L) + lambda * sum(abs(E(:))); if primal_res < tol && dual_res < tol fprintf('ADMM 在第 %d 次迭代收敛\n', k); obj_history = obj_history(1:k); break; end L_old = L; end end逐段说明。L 子问题里，temp = X - E + Y/rho是增广拉格朗日项对 L 的梯度方向，然后调用svt_tensor做奇异值软阈值，阈值取1/rho。E 子问题同理，temp = X - L + Y/rho，然后做元素级软阈值，阈值取lambda/rho。乘子更新Y = Y + rho*(X - L - E)是标准的对偶上升步骤。
收敛判断用了原始残差和对偶残差两个指标。原始残差衡量约束满足程度，对偶残差衡量 L 的变化量。两个都小于tol才认为收敛。实际跑的时候，tol设1e-4到1e-5比较合适，太小会迭代很多次但精度提升有限。
提示：
rho的取值对收敛速度影响很大。rho太大，L 子问题的阈值1/rho太小，去噪不充分；rho太小，E 子问题的阈值lambda/rho太大，噪声被过度剔除。常见做法是设rho = 1，然后调lambda。如果收敛慢，可以试试自适应rho，每 10 次迭代根据残差比调整。
3.4 主脚本的调用流程
main_denoise.m把上面这些模块串起来：
% 加载数据 load('data/sample.mat'); % 假设里面有变量 X_clean 和 X_noisy % 参数设置 lambda = 0.05; % 稀疏惩罚系数 rho = 1; % ADMM 惩罚参数 max_iter = 200; % 最大迭代次数 tol = 1e-5; % 收敛阈值 % 调用 ADMM 求解 [L_denoised, E_noise, obj_hist] = admm_solver(X_noisy, lambda, rho, max_iter, tol); % 评价指标 psnr_val = psnr_ssim(X_clean, L_denoised); fprintf('去噪后 PSNR: %.2f dB\n', psnr_val); % 可视化 figure; subplot(1,3,1); imshow(X_noisy(:,:,1), []); title('含噪'); subplot(1,3,2); imshow(L_denoised(:,:,1), []); title('去噪'); subplot(1,3,3); imshow(X_clean(:,:,1), []); title('干净');这段脚本可以直接跑。如果你的数据不是图像格式，把imshow换成plot或者imagesc就行。psnr_ssim是工具函数，输入干净张量和去噪张量，输出 PSNR 和 SSIM 两个指标。
4. 跑起来：环境配置、数据准备与参数调优
4.1 环境要求与依赖检查
这套代码对 Matlab 版本的要求不高，R2015b 以上基本都能跑。不需要任何额外工具箱——这是它比很多深度学习去噪代码方便的地方。如果你用的是 Octave，大部分函数也能兼容，但fft对三维数组的处理和 Matlab 略有差异，可能需要手动调整tsvd.m里的维度参数。
跑之前先检查两件事：一是core/和utils/是否都在 Matlab 的搜索路径里，二是data/sample.mat是否存在。如果路径不对，main_denoise.m会在load那一步报错。常见做法是在脚本开头加一行addpath(genpath(pwd))，把当前目录及子目录都加进路径。
4.2 用自己的数据替换示例数据
示例数据sample.mat里通常包含两个变量：X_clean（干净张量）和X_noisy（加噪张量）。你要用自己的数据，只需要保证变量名一致，或者改main_denoise.m里的加载逻辑。
数据格式要求是三维数组n1 x n2 x n3。如果你手头是二维矩阵（比如单张图像），需要先扩展成三维。常见做法是用repmat复制几份，或者用reshape把多个样本拼成第三维。比如你有 100 张 256x256 的图像，可以拼成256 x 256 x 100的张量。
% 把多个二维样本拼成三维张量 data = zeros(256, 256, 100); for i = 1:100 data(:,:,i) = load(sprintf('sample_%03d.mat', i)); end X_noisy = data;如果你的数据维度更高（比如四阶张量），这套代码需要改tsvd.m和tprod.m的维度处理逻辑。三阶是这套代码的原生支持范围，四阶以上需要额外做 matricization。
4.3 参数调优的实操方法
lambda是最关键的参数。它控制稀疏项和低秩项之间的权衡。lambda越大，越多能量被归到噪声项 E 里，去噪越激进；lambda越小，越多噪声被保留在 L 里。
我一般用网格搜索：先设lambda从 0.001 到 0.5 取 10 个对数等距值，每个值跑一遍 ADMM，记录 PSNR。然后在这个最优值附近再细化搜索。整个过程大概几分钟，比盲目试参数靠谱得多。
lambda_list = logspace(-3, log10(0.5), 10); psnr_list = zeros(size(lambda_list)); for i = 1:length(lambda_list) [L, ~, ~] = admm_solver(X_noisy, lambda_list(i), 1, 200, 1e-5); psnr_list(i) = psnr_ssim(X_clean, L); end [~, best_idx] = max(psnr_list); fprintf('最优 lambda: %.4f, PSNR: %.2f dB\n', lambda_list(best_idx), psnr_list(best_idx));rho一般不用调，固定 1 就行。但如果你的数据量很大（比如n1*n2*n3 > 1e6），ADMM 收敛会变慢，可以把rho设大一点（比如 5 或 10），同时把max_iter加到 500。
4.4 运行时间与内存估算
三阶张量的大小直接决定内存占用。tsvd需要对n3个切片做矩阵 SVD，每个切片是n1 x n2，所以内存峰值大约是O(n1*n2*n3)的几倍。对于256 x 256 x 100的数据，单精度下大约 25 MB，双精度 50 MB，普通笔记本完全能扛住。但如果你的数据是512 x 512 x 500，内存会到 GB 级别，需要改成分块处理或者用single精度。
运行时间方面，一次 ADMM 迭代的主要开销在tsvd和tprod。256 x 256 x 100的数据，一次迭代大约 0.5 到 1 秒，收敛通常需要 50 到 150 次迭代，总时间 1 到 3 分钟。如果开了parfor并行化tsvd里的频率切片循环，能快 2 到 4 倍。
5. 避坑与排查：跑这套代码时最容易翻车的几个地方
5.1 现象：去噪结果全黑或全白，PSNR 极低
原因通常是lambda设得太大或太小。lambda太大时，E 子问题把几乎所有能量都当成噪声剔除，L 接近零张量；lambda太小时，E 接近零，L 几乎等于原始含噪数据，去噪没效果。
解决方法是先跑一遍lambda扫描，看 PSNR 曲线的峰值位置。如果曲线单调下降，说明lambda的搜索范围整体偏大，把上限调小；如果单调上升，说明范围偏小，把上限调大。
5.2 现象：ADMM 不收敛，迭代到 max_iter 还在震荡
常见原因是rho和lambda的比例不合适。ADMM 的收敛性对rho敏感，rho太小会导致对偶变量更新步长过大，产生震荡。
解决办法是增大rho（比如从 1 调到 5），或者实现自适应rho：每 20 次迭代检查原始残差和对偶残差的比值，如果比值大于 10 就增大rho，小于 0.1 就减小rho。代码里可以加一个简单的自适应逻辑：
if primal_res > 10 * dual_res rho = rho * 2; Y = Y / 2; % 乘子需要同步缩放 elseif dual_res > 10 * primal_res rho = rho / 2; Y = Y * 2; end5.3 现象：tsvd报错「矩阵维度不一致」
这通常是因为输入张量的第三维大小n3和tsvd内部假设的不一致。有些版本的tsvd.m假设n3是奇数，因为 FFT 后频率切片是对称的，只处理一半。如果你的n3是偶数，需要检查代码里是否有floor(n3/2)之类的硬编码。
解决办法是打开tsvd.m，看for i = 1:n3这个循环的范围。如果是1:ceil(n3/2)，说明它只处理了前半部分频率，后半部分需要共轭对称补全。这种情况下，要么把数据第三维补零到奇数，要么改代码支持偶数。
5.4 现象：去噪后的张量出现块状伪影
这是张量分解类方法的典型问题，根源在于 t-SVD 沿第三维做 FFT 时，默认数据在第三维是周期性的。如果你的数据第三维不是周期信号（比如时间序列有趋势项），FFT 会在边界产生频谱泄漏，导致重构后出现块状伪影。
解决办法有两个：一是去趋势，在跑 ADMM 之前先对第三维做差分或减去均值；二是加窗，在 FFT 之前乘一个 Hann 窗，但这样会损失一部分能量，需要做重叠保留处理。我一般用第一种，简单有效。
5.5 现象：换用自己的数据后 PSNR 计算报错
psnr_ssim.m通常假设输入是uint8类型且范围在 0 到 255 之间。如果你的数据是double类型且范围在 0 到 1，直接传进去会得到错误的 PSNR 值。
解决办法是在调用前做类型转换和范围归一化：
if isa(X_clean, 'double') && max(X_clean(:)) <= 1 X_clean = uint8(X_clean * 255); L_denoised = uint8(L_denoised * 255); end psnr_val = psnr_ssim(X_clean, L_denoised);6. 二次开发：把 ADMM 换成近端梯度法，以及一个我常用的验证习惯
这套代码的架构留了很好的扩展点。admm_solver.m是一个独立的求解器，输入输出接口清晰，你完全可以写一个pgd_solver.m用近端梯度法替代它，然后在main_denoise.m里换一行调用就行。
近端梯度法的迭代格式是：
$$\mathcal{L}{k+1} = \text{prox}{\eta |\cdot|_*}\left(\mathcal{L}_k - \eta \nabla g(\mathcal{L}_k)\right)$$
其中 $g(\mathcal{L}) = \lambda |\mathcal{X} - \mathcal{L}|_1$，梯度 $\nabla g$ 在 L1 范数下是次梯度，实际实现时用sign函数近似。步长 $\eta$ 需要满足 $\eta < 2 / L_g$，其中 $L_g$ 是 $g$ 的 Lipschitz 常数。对于 L1 项，$L_g = \lambda$，所以 $\eta < 2/\lambda$。
function [L, obj_history] = pgd_solver(X, lambda, eta, max_iter, tol) % 近端梯度法求解: min ||L||_* + lambda*||X - L||_1 L = X; % 初始化 obj_history = zeros(max_iter, 1); for k = 1:max_iter % 梯度步: g(L) = lambda * ||X - L||_1 grad = lambda * sign(L - X); L_half = L - eta * grad; % 近端步: 核范数近端算子 L_new = svt_tensor(L_half, eta); % 收敛判断 obj_history(k) = tensor_nuclear_norm(L_new) + lambda * sum(abs(X(:) - L_new(:))); if norm(L_new(:) - L(:)) / norm(L(:)) < tol obj_history = obj_history(1:k); break; end L = L_new; end end近端梯度法的优势是实现简单，不需要维护对偶变量；劣势是收敛速度通常比 ADMM 慢，尤其是当问题条件数不好时。我一般会先用 ADMM 跑出结果作为基准，再用 PGD 验证，如果两者 PSNR 差距在 0.5 dB 以内，说明 ADMM 的结果可信。
另一个常用的改造方向是换用不同的张量分解。t-SVD 只是其中一种，还有 CP 分解、Tucker 分解。CP 分解的秩定义更直观，但计算 NP 难，通常用 ALS 交替最小二乘近似。Tucker 分解更灵活，但核范数的定义更复杂。如果你想试 CP 分解，可以把svt_tensor换成 CP-ALS 的软阈值版本，但要注意 CP 分解的秩选择比 t-SVD 更敏感。
验证改造是否成功，我有一个固定习惯：不管换什么求解器或分解方式，都先在gen_test_data.m生成的合成数据上跑一遍。合成数据的好处是干净张量已知，PSNR 可以直接算，而且噪声类型可控（高斯、拉普拉斯、椒盐）。如果合成数据上 PSNR 达标，再换真实数据。真实数据没有 ground truth，只能靠视觉检查和下游任务指标（比如分类准确率）来判断。从那以后我每次改求解器，都强制走一遍「合成数据验证 → 真实数据试跑 → 参数微调」这个流程，省了很多来回折腾的时间。希望帮到你。
本文还有配套的精品资源，点击获取
