www.gusucode.com > 瑞利信道下噪声能量的估计源码程序 > 瑞利信道下噪声能量的估计源码程序/MSP_estimate_version1/ReferenceMethods/Sijbers_Aja.m

    % Sijbers-Aja's method based on eq. 24, page 1400 in
% S. Aja-Fern醤dez, et al., Noise estimation in single-and multiple-coil magnetic resonance data based on statistical models, 
% Magnetic Resonance Imaging, vol. 27(10), pp. 1397-1409, 2009
%
% 05/01/2013
% 
% Tomasz Pieciak
% AGH university of Science and Technology, Krakow, Poland
% pieciak@agh.edu.pl, http://home.agh.edu.pl/pieciak/
%
% ARGUMENTS
%   data - single-coil MRI data
%   window - sliding window size
%   bins - histogram bins
%
% FUNCTION RETURNS
%   sigma - estimated noise level (sigma)
%
% USAGE
%    sigma_estimated = Sijbers_Aja(data, [7, 7], 1000)

function [sigma] = Sijbers_Aja(data, window, bins)


N = window(1) * window(2);
FUNCTION = @(x)( (1/(N-1))*sum(x.^2) ); % function calculates second-order moment
EX2 = colfilt(data, window, 'sliding', FUNCTION);
EX2 = sqrt(EX2);

% initial noise level searching
[histogram_p, histogram_x] = hist(EX2(:), bins);
histogram_p_filtered = filtfilt(ones(1,25), 1, histogram_p); % low-pass filter (window 1x25)
[value, index] = max(histogram_p_filtered);
fc = histogram_x(2*index);

% minimization procedure - least-squares fitting procedure
[n, l] = hist(EX2(EX2 <= fc), bins);  % letters according to eq. 2 in 'Noise estimation in single-and multiple-coil magnetic resonance data based on statistical models', p. 1400


Nk = sum(n);
K = bins;
F_MIN = @(x) (  Nk*log( gammainc(N, (l(1).^2) .* (N./(2*x.^2)), 'upper') - gammainc(N, (l(K).^2) .* (N./(2*x.^2)), 'upper') )   -   sum(n(2:K) .* log( gammainc(N, (l(1:K-1).^2) .* (N./(2*x.^2)), 'upper') - gammainc(N, (l(2:K).^2) .* (N./(2*x.^2)), 'upper') ))  );   
                        % according to eq. 24, p.1340
                        % x - sigma

sigma0 = histogram_x(index);
[sigma, fval, exitflag] = fminsearch(F_MIN, sigma0);

if(exitflag == 0)
    fprintf('<strong> Sijbers-Aja error! </strong>\n')
end