image thumbnail

SpinEchoMLE

by

 

11 Jul 2013 (Updated )

Maximum likelihood estimator for spin echo time decay constant.

SpinEchoMLE(Y,two_tau,Avg,x,opt)
function [T2] = SpinEchoMLE(Y,two_tau,Avg,x,opt)
% description:  
%                Maximum likelihood estimator for parameter T2 in the model
%                         Y = a * x * exp(-two_tau/T2)' + G
%                where G is zero-mean gaussian noise, and  
%           T2 - unknown time-decay constant (scalar)
%            x - unit norm vector for the echo waveform (M x 1)
%                with unknown amplitude, a; MLE is computed for either
%                of two cases: x known or unknown.            
%      two_tau - column vector of echo times (N x 1).
%                Thus, column k of matrix Y is the measured echo at 
%                time two_tau(k).
% inputs:
%            Y - M x N matrix of echo data, one sampled echo per column
%      two_tau - N x 1 column vector of echo times
%          Avg - N x 1 list of number of pulses averaged at each echo time
%                (default is equal number of pulses averaged at each echo)
%            x - M x 1 column vector of known echo shape (optional)
%          opt - Optimization parameters for line search (optional)
% output:
%           T2 - MLE for decay time.
%                time unit is merely consistent with usage in two_tau
% examples of usage:
%           for unknown echo shape and equal averaging per pulse
%               >> T2 = SpinEchoMLE(Y,two_tau)
%           for unknown echo shape and unequal averaging per pulse
%               >> T2 = SpinEchoMLE(Y,two_tau,Avg)
%           for known echo shape with default 1D search parameters
%               >> T2 = SpinEchoMLE(Y,two_tau,Avg,x)
%           for known echo shape and equal averaging, default search
%               >> T2 = SpinEchoMLE(Y,two_tau,[],x)
% authors:
%           lcp and jna, 11-July-2013
%           Copyright 2013
%           Revision 1.0
%           Reference: "Estimation of spin-echo relaxation time," Golub,
%                      Potter, Ash, Blank and Ahmad, 2013.


%% defaults on inputs
switch nargin
    case 1
        error('myApp:argChk',...
            'SpinEchoMLE requires data array and list of echo times');
    case 2
        Avg = ones(size(two_tau));
        x=[];
        opt = optimset('TolX',1e-6,'Display','off','MaxFunEvals',2000, ...
            'Algorithm','levenberg-marquardt');%default values
    case 3
        x=[];
        opt = optimset('TolX',1e-6,'Display','off','MaxFunEvals',2000, ...
            'Algorithm','levenberg-marquardt');%default values
        if((min(Avg) <= 0) || (norm(Avg - round(Avg)) ~= 0)  ...
                || (length(Avg) ~= length(two_tau)))
            error('myApp:argChk', ...
            'Invalid input: Avg should be %d positive integers',...
                length(two_tau));
        end
    case 4
        opt = optimset('TolX',1e-6,'Display','off','MaxFunEvals',2000, ...
            'Algorithm','levenberg-marquardt');%default values
end
if isempty(opt) || ~exist('opt','var')  %default values
    opt = optimset('TolX',1e-6,'Display','off','MaxFunEvals',2000, ...
            'Algorithm','levenberg-marquardt');
end
if isempty(Avg)
    Avg = ones(size(two_tau));
end
two_tau=two_tau(:); %force column


%% MLE
Winv = diag(Avg);
if isempty(x), %unknown echo shape
    Q=Winv*(Y'*Y)*Winv; 
else %known echo shape
    x = x(:)/norm(x); %force column and normalize length
    Q=Winv*Y'*(x*x')*Y*Winv; 
end
% Initialize T2 optimization at mean echo time
SEfun = @(T) exp(-two_tau/T)'*Q*exp(-two_tau/T) / ...
    ( - exp(-two_tau/T)'*Winv*exp(-two_tau/T) );
T2 = fminsearch(@(T) SEfun(T),mean(two_tau), opt);
end

Contact us