基于稀疏辅助信号平滑的心电信号降噪方法（Matlab R2021B）

基于形态成分分析理论（MCA）的稀疏辅助信号分解方法是由信号的形态多样性来分解信号中添加性的混合信号成分，它最早被应用在图像处理领域，后来被引入到一维信号的处理中。

在基于MCA稀疏辅助的信号分析模型中，总变差方法TV是其中一个原型，稀疏辅助平滑方法结合并统一了传统的LTI低通滤波和总变差算法，兼具LTI低通滤波和总变差算法的优势，稀疏辅助平滑降噪的适用性更广泛，降噪的表现更好。已有研究说明，稀疏辅助平滑降噪相比低通滤波器能够有效保留瞬态冲击的幅值。

鉴于此，提出一种基于稀疏辅助信号平滑的心电信号降噪方法，运行环境为Matlab R2021B。

function [A, B, B1, D, a, b, b1, H1norm, HTH1norm] = Lowpass(deg, fc, N, K)
% [A, B, B1] = ABfilt(d, fc, N, K)
%
% Banded matrices for zero-phase high-pass recursive filter.
% The filter is H = inv(A) * B where the matrices A and B
% are created as 'sparse' matrices.
%
% INPUT
%   d  : degree of filter is 2d
%   fc : cut-off frequency (normalized frequency, 0 < fc < 0.5)
%   N  : length of signal
%   K  : order of difference matrix D (need 1 <= K <= 2*d) (default K = 1)
%
% OUTPUT
%   A, B, B1 : banded filter matrices
%       with B = B1*D where D is the K-th order difference (up to sign)
%
% Use [A, B, B1, D, a, b, b1] = ABfilt(...) to return
% filter coefficient vectors a, b, b1.

if nargin < 4
    K = 1;
end

if K > 2*deg
    error('ABfilt: K > 2*d')
end

omc = 2*pi*fc;
t = ((1-cos(omc))/(1+cos(omc)))^deg;

% Define p such that P(z)P(1/z) = B(z), i.e., P'*P = B
p = 1;
for k = 1:deg
    p = conv(p, [-1 1]);
end
P = spdiags( p(ones(N-deg,1), :), 0:deg, N-deg, N);        % banded matrix

B = P' * P;
% pp = conv(p, p(end:-1:1));

q = sqrt(t);
for i = 1:deg
    q = conv(q, [1 1]);
end
Q = spdiags( q(ones(N-deg,1), :), 0:deg, N-deg, N);    % banded matrix

A = P'*P + Q'*Q;

if K <= deg
    d = 1;
    for i = 1:K
        d = conv(d, [-1 1]);
    end
    D = spdiags(d(ones(N,1), :), 0:K, N-K, N);               % D: banded matrix
    
    p1 = deconv(p, d);
    P1 = spdiags( p1(ones(N-deg,1), :), 0:deg-K, N-deg, N-K);  % banded matrix
    B1 = P' * P1;

    b1 = conv(p1, p(end:-1:1));
else
    % deg < K <= 2*deg
    K2 = 2*deg - K;
    d = 1;
    for i = 1:K2
        d = conv(d, [-1 1]);
    end
    B1 = spdiags(d(ones(N,1), :), 0:K2, N-K2, N)';               % D: banded matrix
    
    p1 = deconv(p, d);
    D1 = spdiags( p1(ones(N-deg,1), :), 0:deg-K2, N-deg, N-K2);  % banded matrix
    D = D1'*P;

    b1 = d;

end

a = conv(p, p(end:-1:1)) + conv(q, q(end:-1:1));
b = conv(p, p(end:-1:1));

% verify that B = B1*D
err = B - B1*D;
mae = max(abs(err(:)));
if mae > 1e-10
    disp('Error in ABfilt (B1*D not equal to B)')
end


% Calculate filter norms

imp = zeros(size(B1,2), 1);
imp(round(N/2)) = 1;                        % imp : impulse signal (located at center to avoid transients)

h1 = A \ (B1 * imp);
H1norm = sqrt( sum( abs( h1 ).^2 ) );       % norm of filter inv(A)*B1

hh = B' * ((A*A') \ (B1 * imp));
%完整代码：https://mbd.pub/o/bread/mbd-ZZaZmJtx

HTH1norm = sqrt( sum( abs( hh ).^2 ) );     % norm of filter B'*inv(A*A')*B1