www.gusucode.com > 瑞利信道下噪声能量的估计源码程序 > 瑞利信道下噪声能量的估计源码程序/MSP_estimate_version1/ReferenceMethods/Brummer_Aja.m
% Brummer-Aja's method based on eq. 20, 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 % % 31/12/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 = Brummer_Aja(data, [7, 7], 1000) function [sigma] = Brummer_Aja(data, window, bins) N = window(1) * window(2); FUNCTION = @(x)( (1/N)*sum(x) ); % function calculates first-order moment EX = colfilt(data, window, 'sliding', FUNCTION); % initial noise level searching [histogram_p, histogram_x] = hist(EX(:), 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 [h, l] = hist(EX(EX <= fc), bins); % letters according to eq. 20 in 'Noise estimation in single-and multiple-coil magnetic resonance data based on statistical models', p. 1400 % according to eq. 20, p. 1400, x(1) - sigma, x(2) - K N_gamma_N = N^N ./ gamma(N); F_MIN = @(x) sum( (h - x(2) .* ((l.^(2*N-1)) ./ (2^(N-1) .* (pi/4 .* x(1).^2).^N)) .* N_gamma_N .* exp(-(N .* l.^2) ./ (2*pi/4 .* x(1).^2)) ).^2); sigma0 = histogram_x(index); K0 = 1; [sigma_K, fval, exitflag] = fminsearch(F_MIN, [sigma0, K0]); sigma = sigma_K(1); if(exitflag == 0) fprintf('<strong> Brummer-Aja error! </strong>\n') end