Why does phase randomizing the Fourier transform of a data set render it Gaussian? Let's say I have a data set $s_n$.  I take the Fourier transform of this data set to obtain $\tilde{s}_n$. I randomize the complex phases of $\tilde{s}_n$ and I take an inverse transform to obtain $s'_n$. This way $s_n$ and $s'_n$ have the same power spectrum.  On the other hand, the PDF of $s'_n$ is Gaussian (at least asymptotically) irrespective of the PDF of $s_n$.  Why is this so?
 A: This is a hand-waving answer that might provide some intuition.
Given a data set $(x_0, x_1, \ldots x_{N-1})$ which is a set of $N$ (independent) samples from an unknown distribution (which we assume
has zero mean and finite varince), we take its
Discrete Fourier Transform (DFT) to obtain a complex vector 
$(X_0, X_1, \ldots, X_{N-1})$ where
\begin{align}
\Re(X_k) &= \frac{1}{\sqrt{N}}\sum_{m=0}^{N-1} x_m \cos(2\pi km/N)\\
\Im(X_k) &= -\frac{1}{\sqrt{N}}\sum_{m=0}^{N-1} x_m \sin(2\pi km/N)
\end{align}
Vigorous hand-waving and the murmuring of shibboleths that sound
a lot like "Sea Elle Tea" allow us to suggest that maybe, just 
maybe, $\Re(X_k)$ and $\Im(X_k)$ can be treated as (independent)
zero-mean normal random variables (approximately). Thus, $|X_k|$ is
Rayleigh distributed (approximately) and when you replace
$\angle X_k$ in $X_k = |X_k|\exp({i\angle X_k})$ by $\Theta_k \sim U[0,2\pi)$, the alleged resemblance of $\Re(\hat{X}_k)$ and 
$\Im(\hat{X}_k)$ to zero-mean normal random variables only
increases. (For best results, choose $\Theta_{N-k} = 2\pi-\Theta_k$
in order to maintain the necessary conjugate symmetry of the DFT 
result.)  Finally, apply the inverse DFT transformation (giving you
linear combinations of normal random variables) and
shout "See Elle Tea" in a more strident voice as you insist
to all and sundry that $(\hat{x}_0, \hat{x}_1, \ldots, \hat{x}_{N-1})$
sure looks like $N$ samples from a normal distribution
