非局部均值(NLM)信号去噪方法.rar

function [denoisedSig,debug] = NLM(signal,lambda,P,PatchHW) % Implements fast NLM method of Darbon et al, for a 1-D signal % Inputs: % signal: input signal (vector) % lambda: Gaussian scale factor % P: max search distance % PatchHW: patch half-width % Output: % denoisedSig: the NLM-denoised signal % debug: structure containing various quantitities that can help debug % Reference:Tracey Brian H, Miller Eric L. Nonlocal means denoising of ECG signals[J]. IEEE Transactions on Biomedical Engineering 2012,59(9):2383-2386 if length(P)==1 % scalar has been entered; expand into patch sample index vector Pvec = -P:P; else Pvec = P; % use the vector that has been input end debug=[]; N = length(signal); denoisedSig = NaN*ones(size(signal)); % to simpify, don't bother denoising edges iStart=1+PatchHW+1; iEnd = N-PatchHW; denoisedSig(iStart:iEnd) = 0; debug.iStart = iStart; debug.iEnd = iEnd; % initialize weight normalization Z = zeros(size(signal)); cnt = zeros(size(signal)); % convert lambda value to 'h', denominator, as in original Buades papers Npatch = 2*PatchHW+1; h = 2*Npatch*lambda^2; for idx = Pvec % loop over all possible differences: s-t % do summation over p - Eq. 3 in Darbon k=1:N; kplus = k+idx; igood = find(kplus>0 & kplus<=N); % ignore OOB data; we could also handle it SSD=zeros(size(k)); SSD(igood) = (signal(k(igood))-signal(kplus(igood))).^2; Sdx = cumsum(SSD); for ii=iStart:iEnd % loop over all points 's' distance = Sdx(ii+PatchHW) - Sdx(ii-PatchHW-1); % Eq 4; this is in place of point-by-point MSE % but note the -1; we want to icnlude the point ii-iPatchHW w = exp(-distance/h); %Eq 2 in Darbon t = ii+idx; % in the papers, this is not made explicit if t>1 && t<=N denoisedSig(ii) = denoisedSig(ii) + w*signal(t); Z(ii) = Z(ii) + w; cnt(ii) = cnt(ii)+1; end end end % loop over shifts % now apply normalization denoisedSig = denoisedSig./(Z+eps); debug.Z = Z; end