Previously, I showed how to whiten a matrix in Matlab. This involves finding the inverse square root of the covariance matrix of a set of observations, which is prohibitively expensive when the observations are high-dimensional – for instance, high-resolution natural images.

Thankfully, it’s possible to whiten a set of natural images approximately by multiplying the data by the inverse of the square root of the mean power spectrum in the Fourier domain. Assume that X is a matrix of size (nimages,height,width), then:

mX = bsxfun(@minus,X,mean(X)); %remove mean
fX = fft(fft(mX,[],2),[],3); %fourier transform of the images
spectr = sqrt(mean(abs(fX).^2)); %Mean spectrum
wX = ifft(ifft(bsxfun(@times,fX,1./spectr),[],2),[],3); %whitened X

Why this works

It’s not obvious at all that this should give similar results to the other method, which involved finding the inverse square root of the covariance matrix M. Let’s assume that the images are stationary – that is, the second order statistics are similar at any point in the image. That means that the elements of the covariance matrix should only depend on the relative position between two pixels. And of course, M is a symmetric matrix.

That means that M is a matrix such that the product Mx performs the convolution between x and the kernel associated with M – call this kernel m. That kernel is actually the sum of the autocorrelation of the images. If , then by the convolution theorem, . The Fourier transform F(x) is a linear mapping from an n-dimensional complex vector to another n-dimensional complex vector. Therefore, , where is an n-by-n dimensional matrix.

Then, we have that . If we take images X and transform them via , then it’s easy to verify that the covariance of Y will now be a scalar times the identity matrix. Therefore, modulo edge effects, which break the stationarity requirement, images can be whitened by dividing them in the Fourier domain by the inverse of the root of the mean power spectrum.