Fourier bases for a stationary signal & relation to PCA for natural images Why does PCA of a translation-invariant signal give a Fourier basis?  
I've found proofs for this, but I'd love some intuitions.  
Any help is greatly appreciated!

EDIT:
Sorry that this question was unclear.
I study natural image statistics and here's an example of the kind of paragraph that inspired the question (from here):

When applying PCA to signals such as images, it is commonly assumed that the statistical properties of the image are translation invariant (also known as stationary). Specifically, one assumes that the correlation of the intensity at two locations in the image depends only on the displacement between the locations, and not on their absolute locations. In this case, the sinusoidal basis functions of the Fourier transform are guaranteed to be a valid set of principal component axes (although, as before, this set need not be unique). The variance along each of these axes is simply the Fourier power spectrum. Whitening may be achieved by computing the Fourier transform, dividing each frequency component by the square root of its variance, and (optionally) computing the inverse Fourier transform. This is further discussed below. [emphasis mine]

I don't need a formal proof, I just want an intuition and I feel like I'm getting near one, but I'd love any further insight.
Here's what I've got so far:
First: You draw a bunch of natural images from the world -- just take for granted that they are translation invariant, which means for the purposes here that the correlation between two pixels will only depend on distance between pixels, not on absolute position.  Also, for the sake of simplicity, let's just imagine we're only getting 1d images [i.e., slices through an image].
Sidenote: In the 1d case, the fact that correlation just depends on distance between pixels, not on absolute position means that the covariance matrix of these samples will have a diagonal structure -- the covariance along the diagonal will obviously be the same, but also the covariance along each off diagonal will be the same [i.e., cov(x1,x2) = cov(x2,x3) = cov(x3,x4), etc.].  In the case of natural images where correlation decreases monotonically with distance, the covariation will decrease monotonically as you move away from the diagonal, but you can imagine a kind of signal where there is more periodic structure -- the rest of the intuitions don't depend on which regime we're in.
Second: The Fourier basis is a basis that only considers relative, not absolute, position.  More precisely, the eigenfunctions of the translation operator are the Fourier basis functions: $e_k(t) = \frac{1}{T} e^{i\frac{2\pi kt}{T}}$.  Let's say you have a translation operator, $\mathcal{T}$, which just translates functions by a.  Thus: 
$$\mathcal{T} e_k(t) = e_k (t + a) = \frac{1}{T} e^{i\frac{2\pi k (t + a)}{T}}$$
$$ = \frac{1}{T} e^{i\frac{2\pi k a}{T}}  e^{i\frac{2\pi k t}{T}}$$
$$ = e^{i\frac{2\pi k a}{T}}  e_k(t)$$
Thus, by translating the Fourier function, we get out the Fourier function multiplied by some constant [its eigenvalue].  
Working through this gave me a flavor for why the Fourier basis is useful for translation-invariant signals, but as you can see this is a little rickety -- if anyone has further insight, I'd really appreciate it.
And one final specific question -- I get why Fourier bases are useful for a stationary signal, but there still is a question about PCA -- why is the orthonormal basis that captures maximal variance Fourier-like?
 A: I am afraid this will not fully answer your question, but I would still like to write it to give you some keywords and mainly hoping to sparkle further discussion. I will be happy if anybody provides a better answer.
Having said that, I disagree with what @whuber wrote in the comments above ("in no circumstances does a Fourier basis act like PCA" etc.). I think you are right in observing that PCA on smooth more-or-less-translation-invariant signals (or images) results in Fourier-like components. This is a real effect and it must have an explanation.
Mathematically, the covariance matrix with "diagonal structure" that you describe in your question is called Toeplitz matrix. Toeplitz matrices have many nice properties, one of them being that in the limit of infinite dimension (Toeplitz operator) the eigenvectors are sines and cosines with increasing frequencies. I think what you are asking to is to provide an intuitive explanation of this fact. I can't.
Of course sines and cosines are also eigenvectors of the Fourier operator, as you rightly observe. Which means that Toeplitz and Fourier operators are closely related. I would be very happy if somebody more savvy in math would explain how it works (I strongly suspect that with the correct perspective it becomes almost trivial).
Finally, if you are interested in a non-technical description of this phenomenon, I can recommend this paper Interpreting principal component analyses of spatial
population genetic variation from Nature Genetics. The authors remark that "sinusoidal mathematical artifacts ... arise generally when PCA is applied to spatial data" and give some nice examples.
