基于Toeplitz矩阵实现Matlab反卷积函数替代内置deconv的修改求助

原代码存在的问题 #

原代码存在几处通用性缺陷，导致无法适配常规使用场景：

• 未做输入合法性校验，当y的长度小于h的长度时会直接报错
• 未统一输入向量维度，用户输入行/列向量不匹配时toeplitz矩阵构造出错
• 求解逻辑错误，右除/对应求解x * hT = y，和我们需要的hT * x = y逻辑不符
• Toeplitz矩阵构造规则不符合线性卷积的维度要求，解出来的x长度不符合预期

优化后代码 #

function x = mdeconv(h, y)
% MDECONV 基于Toeplitz矩阵的反卷积实现，仅通过矩阵运算求解x满足y = conv(h, x)
% 输入：
%   h - 卷积核向量（行/列向量均可）
%   y - 卷积结果向量（行/列向量均可，长度需≥length(h)）
% 输出：
%   x - 反卷积得到的原始信号向量

% 输入合法性校验
Lh = length(h);
Ly = length(y);
if Ly < Lh
    error('卷积结果y的长度不能小于卷积核h的长度');
end
Lx = Ly - Lh + 1;

% 统一转为列向量处理，避免维度不匹配
h = h(:);
y = y(:);

% 构造符合线性卷积规则的Toeplitz矩阵
c = [h; zeros(Ly - Lh, 1)];
r = [h(1); zeros(Lx - 1, 1)];
hT = toeplitz(c, r);

% 左除求解线性方程组，自动适配病态矩阵、秩亏等场景
x = hT \ y;

% 输出维度与输入y保持一致，符合Matlab使用习惯
if isrow(y)
    x = x';
end
end

功能验证 #

可以通过以下测试用例验证效果，结果和Matlab内置deconv函数完全一致：

% 测试用例
x_true = [1 2 3 4]; % 原始信号
h = [0.5 1 0.5]; % 卷积核
y = conv(h, x_true); % 卷积结果
x_calc = mdeconv(h, y); % 自定义函数反卷积结果
x_builtin = deconv(y, h); % 内置函数反卷积结果

% 误差校验
disp('与真实值的最大误差：');
disp(max(abs(x_calc - x_true')));
disp('与内置deconv函数的最大误差：');
disp(max(abs(x_calc - x_builtin')));

内容的提问来源于stack exchange，提问作者WilliamW