首先，我们需要确保数据的均值（近似）为零。对于自然图像，我们通过减去每个图像块的均值（近似地）来达到这一目标。为此，我们计算每个图像块的均值，并从每个图像块中减去它的均值。Matlab实现如下：
avg = mean(x, 1); % 分别为每个图像块计算像素强度的均值。
x = x - repmat(avg, size(x, 1), 1);
下面计算得出样本数据x的协方差矩阵Σ：
sigma = x * x' / size(x, 2); % x 是一个n × m的矩阵，每一列表示一个训练样本
接下来，计算 Σ 的特征向量。你可以使用Matlab的 eig 函数来计算。但是由于 Σ 是对称半正定的矩阵，用 svd 函数在数值计算上更加稳定。
[U,S,V] = svd(sigma);
矩阵 U 将包含 Sigma 的特征向量（一个特征向量一列，从主向量开始排序），矩阵S 对角线上的元素将包含对应的特征值（同样降序排列）。矩阵V等于U的转置，可以忽略。
（注意 svd 函数实际上计算的是一个矩阵的奇异值和奇异向量，就对称半正定矩阵的特殊情况来说，它们对应于特征值和特征向量 ）
最后，计算 xRot和xTilde ：
xRot = U' * x; % 数据旋转后的结果。
xTilde = U(:,1:k)' * x; % 数据降维后的结果，这里k希望保留的特征向量的数目。
计算PCA白化后的数据:
xPCAwhite = diag(1./sqrt(diag(S) + epsilon)) * U' * x; % S的对角线包括了特征值λ, epsilon(ε)为大约等于10^-5的常数，用于进行正则化
计算ZCA白化后的数据:
xZCAwhite = U * diag(1./sqrt(diag(S) + epsilon)) * U' * x;
实践中实现PCA白化或ZCA白化时，有时一些特征值λ在数值上接近于0，这样在缩放步骤时将导致除以一个接近0的值，因而在实践中，我们使用少量的正则化实现这个缩放过程，即在取平方根和倒数之前给特征值加上一个很小的常数ε ：
当x 在区间[-1,1] 上时, 一般取值为10^-5 。
对图像来说, 这里加上ε ，对输入图像也有一些平滑(或低通滤波)的作用。这样处理还能消除在图像的像素信息获取过程中产生的噪声，改善学习到的特征。
注：上述代码引自Standford大学Andrew Ng的文章