function [T2,X,K] = tsvdfy(tau2,M,N,B)
%tau2：信号采集时刻 M：CPMG信号幅值 N：迭代次数 B:布点数
%T2:布点数 X：反演结果 K：反演核心矩阵
T2 = logspace(-2,4,B);%T2布点
[T,TAU2] = meshgrid(T2,tau2);
K = exp(-TAU2./T); %生成反演核
SNR = (M(1)+M(2))/2/std(M(end-10:end));%计算信噪比
[a,~] = size(K);
K = [K ones(a,1)]; %将基底信号拟合进去
[U,S,V] = svd(K,'econ'); %奇异值分解
f = S(1,1)/SNR; %计算截断阈值
Thresh = sum(sum(S>f)); %计算截断位置
% Thresh = Thresh-2;
sc = diag(S(1:Thresh,1:Thresh)); %奇异值截断
X = V(:,1:Thresh)*diag(1./sc)*U(:,1:Thresh)'*M; %得到初解
%非负约束
k = 1;
while min(X) < 0 && k<= N
X(X<0) = 0;
AY = M-K*X; %计算误差
ef = norm(AY)/norm(M);
if ef<0.01
break;
end
AX = V(:,1:Thresh)*diag(1./sc)*U(:,1:Thresh)'*AY; %误差分摊到X
X = X+AX;
k = k+1;
end
X(X<0) = 0;
X = X(1:end-1); %去除基线
K = K(:,1:end-1);
评论0