function [m, sd, amp, phase] = mri_phasemap(images, freq, smooth);

% function [m, sd, amp, phase] = mri_phasemap(images, freq, smooth);
%
% Aug 2003, Farshad Moradi
% California Institute of Technology
%

m = 0;
x = 0;
sd = 0;
N = length(images);
if nargin<3, smooth = 0; end;

for i= 1: N,
    d = images(i).data;
    x = x + d * exp(-j*2*pi*freq*(i-1)/N);
    m = m + d/N;
    sd = sd + d.^2/N;
end;

sd = sd-m.^2;

if smooth,
	[xx,yy] = meshgrid(-3:1.5:3); b = exp(-(xx.^2+yy.^2)/2); b = b / sum(b(:));
	disp('smoothing images...');
    for i=1:size(x,3),
		x(:,:,i) = conv2(x(:,:,i),b,'same');
		sd(:,:,i) = conv2(sd(:,:,i),b,'same');
    end;
end;

sd = sqrt(sd);
amp = abs(x);
phase = angle(x);