function [Z, V] = whiten(x)
epsilon = 1e-4;

% x is the input patch data of size
% z is the ZCA transformed data. The dimenison of z = x.

x_mean = mean(x, 1);
x = bsxfun(@minus, x, x_mean);
sigma = x * x' / size(x,2);
[U, S, ~] = svd(sigma);
V = U * diag(1 ./ sqrt(diag(S) + epsilon)) * U';
Z = V * x;
end