Output of a System with a White Gaussian Process as an Input I want to present a question which seems to contradict a theorem I know.
Hence I guess I'm missing something and would be happy to understand what I'm missing.
Here's is the system:

The signals $ {X}_{1}, {X}_{2} $ are Independent White Gaussian Noise Process.
The filters $ {H}_{1}(f), {H}_{2}(f) $ are identical, namely, $ {H}_{1}(f) = {H}_{2} $.
The phase, $ \theta \sim U[0, 2\pi] $ and independent of any other parameter or signal.
The questions are

*

*Are $ {Y}_{1}(t), \  {Y}_{2}(t) $ Gaussian Processes?

*Is $ r(t) $ Gaussian Process?

The answers, according to a book are $ {Y}_{1}, \ {Y}_{2} $ aren't Gaussian Processes and $ r(t) $ is a Gaussian Process.
Yet $ {Y}_{1}, \ {Y}_{2} $ are obviously independent (Edit: They are dependent, see answers and comments below) and according to Cramer's Theorem:
http://en.wikipedia.org/wiki/Normal_distribution#Infinite_divisibility_and_Cram.C3.A9r.27s_theorem
If the sum of two independent R.V. is Gaussian then each of them must be Gaussian.
Yet it's pretty clear to me that Y1 / Y2 aren't Gaussian.
I wrote the following MATLAB Simulation:
numRealizations = 1e4;
numSamples      = 1e3;

h1Coef = ones(10, 1) / 10;
h2Coef = h1Coef;

harmonicFreq    = 10; %<! [Hz]
timeVector      = [0:(numSamples - 1)];
randomPhaseVec  = 2 * pi * rand(numRealizations, 1);

% White Noise is Ergodic
% Each Row is a realization (Over the Time Vector)
x1SignalMatrix  = randn(numRealizations, numSamples);
x2SignalMatrix  = randn(numRealizations, numSamples);

% Filtering along the columns
z1SignalMatrix = filter(h1Coef, 1, x1SignalMatrix, [], 2);
z2SignalMatrix = filter(h2Coef, 1, x2SignalMatrix, [], 2);

y1SignalMatrix = z1SignalMatrix .* cos(bsxfun(@plus, (2 * pi * harmonicFreq * timeVector), randomPhaseVec));
y2SignalMatrix = z2SignalMatrix .* sin(bsxfun(@plus, (2 * pi * harmonicFreq * timeVector), randomPhaseVec));

rSignalMatrix = y1SignalMatrix + y2SignalMatrix;

timeIndex           = randi([1, numSamples], [1, 1]);
realizationIndex    = randi([1, numRealizations], [1, 1]);

figure;
hist(x1SignalMatrix(:, timeIndex), 500);
figure;
normplot(x1SignalMatrix(:, timeIndex));


figure;
hist(z1SignalMatrix(:, timeIndex), 500);
figure;
normplot(z1SignalMatrix(:, timeIndex));


figure;
hist(y1SignalMatrix(:, timeIndex), 500);
figure;
normplot(y1SignalMatrix(:, timeIndex));

figure;
hist(y1SignalMatrix(realizationIndex, :), 500);
figure;
normplot(y1SignalMatrix(realizationIndex, :));


figure;
hist(rSignalMatrix(:, timeIndex), 500);
figure;
normplot(rSignalMatrix(:, timeIndex));

figure;
hist(rSignalMatrix(realizationIndex, :), 500);
figure;
normplot(rSignalMatrix(realizationIndex, :));

This is the distribution of Y1 over the Ensemble:

Which is clearly not Gaussian as I expected.
This is the distribution of Y2 over the time (One Realization):

Now, the results for $ r(t) $:
The Distribution of the Ensemble:

The Distribution of 1 Realization:

Well, for Y1, I can see why the ensemble won't Gaussian as for each 't' the sample is multiplied by a different factor hence each of them, though is Gaussian, is coming from a different distribution.
I even wrote a proof why it doesn't comply with Gaussian Distribution properties (The fourth Moment) - http://imgur.com/l017RrR.
Yet I can't get how $ r(t) $ is Gaussian.
What am I missing?
Thank You.
 A: I believe that, stripped to its essentials, this question is asking the following.

Let $X, Y, \Theta$ denote independent random variables and suppose
  that $X, Y \sim N(0,1)$. Then,
  
  
*
  
*Is $X\cos \Theta + Y\sin \Theta$ a normal random variable?
  
*Are $X\cos \Theta$ and $Y\sin\Theta$ normal random variables? Are they independent random variables? Are they uncorrelated random variables?
  
*Do any of the answers to the above questions change if $\Theta$ has a
  specific distribution, e.g. $\Theta \sim U[0,2\pi)$?


Let $Z = X\cos\Theta+Y\sin\Theta$. Since $X$ and $Y$ are independent of $\Theta$, their conditional
joint density given that $\Theta = \theta$ is the same as their
unconditional joint density.  Thus,
$$f_{X,Y\mid \Theta}(x,y\mid \Theta=\theta) 
= \frac{1}{2\pi}\exp\left(-\frac{x^2+y^2}{2}\right), 
~ -\infty < x, y < \infty.\tag{1}$$


Edit added in response to OP's question:
Proof of claim: If $(X,Y)$ is given to be independent
of $\Theta$, we have that
$$f_{X,Y,\Theta}(x,y,\theta) = f_{X,Y}(x,y)f_\Theta(\theta), \qquad
\text{definition of independence}$$
Consequently,
$$\begin{align}
f_{X,Y\mid \Theta}(x,y\mid \Theta=\theta) 
&= \frac{f_{X,Y,\Theta}(x,y,\theta)}{f_\Theta(\theta)}
&\text{definition of conditional density}\\
&=  \frac{f_{X,Y}(x,y)f_\Theta(\theta)}{f_\Theta(\theta)}
&\text{substitute from above}\\
&= f_{X,Y}(x,y)\end{align}$$
Thus, $f_{X,Y\mid \Theta}(x,y\mid \Theta=\theta)$,
the conditional joint density of $(X,Y)$ given $\Theta = \theta$
is the same as $f_{X,Y}(x,y)$, the unconditional joint density
of random variables $X$ and $Y$ that have never heard of $\Theta$
and don't know that they are independent of $\Theta$.
Note that it is not necessary to assume that $X$ and $Y$ are independent
of each other: just that they are independent of $\Theta$.
But, in this case, $X$ and $Y$ are assumed to be independent 
standard normal random variables with joint density $f_{X,Y}(x,y)$
given by the right side of $(1)$, and thus 
$f_{X,Y\mid \Theta = \theta}(x,y\mid\Theta = \theta)$
equals the right side of $(1)$.
For future reference, note that since
$$\begin{align}
f_{X,Y\mid \Theta = \theta}(x,y\mid \Theta = \theta)
&= \frac{1}{2\pi}\exp\left(-\frac{x^2+y^2}{2}\right),\\
&= \frac{1}{\sqrt{2\pi}}\exp\left(-\frac{x^2}{2}\right)
\times \frac{1}{\sqrt{2\pi}}\exp\left(-\frac{y^2}{2}\right),\\
&= f_X(x)f_Y(y)\\
&=f_{X\mid \Theta=\theta}(x\mid \Theta = \theta)
f_{Y\mid \Theta=\theta}(y\mid \Theta =\theta)
\end{align}$$
where the last equality follows from the observations that, by
independence,
$$\begin{align}
f_{X\mid \Theta=\theta}(x\mid \Theta = \theta) 
&= \frac{f_{X,\Theta}(x,\theta)}{f_\Theta(\theta)}
= \frac{f_X(x)f_\Theta(\theta)}{f_\Theta(\theta)}=f_X(x)\\
f_{Y\mid \Theta=\theta}(y\mid \Theta = \theta) 
&= \frac{f_{Y,\Theta}(y,\theta)}{f_\Theta(\theta)}
= \frac{f_Y(y)f_\Theta(\theta)}{f_\Theta(\theta)}=f_Y(y)
\end{align}$$
Thu, $X$ and $Y$ are conditionally independent of each other given $\theta$,
and of course they are unconditionally independent as well.
End of edit added in response to OP's question: 


Thus, conditioned on $\Theta = \theta$, 
$Z = X\cos\theta + Y\sin\theta$ is a weighted sum of two independent
standard normal random variables. So, the conditional distribution of
$Z$ is a normal distribution, and in fact, a standard normal distribution.
This follows readily from
$$\begin{align}
E[Z\mid \Theta=\theta] &= E[X\cos\theta+Y\sin\theta] = 0,\\
\operatorname{var}(Z\mid\Theta=\theta) &= \operatorname{var}(X)\cos^2\theta
+ \operatorname{var}(Y)\sin^2\theta = 1\cdot\cos^2\theta + 1\cdot\sin^2\theta
= 1.
\end{align}$$

Since the conditional density of $Z$ given $\Theta = \theta$ is the
  standard normal density regardless of the value of $\Theta$, it follows
  that 
$\quad$ the unconditional density of $Z = X\cos\Theta + Y\sin\Theta$
  is the standard normal density.
Note that the distribution of $\Theta$ is immaterial. Indeed, $\Theta$
  could be a degenerate random variable that takes on a fixed value
  $\theta$ with probability $1$, and the unconditional density
  of $Z$ would still be the standard normal density. In particular, 
  $X\cos\Theta+Y\sin\Theta$ is a standard normal random variable
  when $\Theta \sim U[0,2\pi)$.


Let $\hat{X} = X\cos\Theta$ and $\hat{Y} = Y\sin\Theta$.
These are uncorrelated random variables since independence
of $X,Y,\Theta$ gives us that 
$$\begin{align}
E[\hat{X}] &= E[X\cos\Theta] = E[X]E[\cos\Theta]=0\cdot E[\cos\Theta] = 0,\\
E[\hat{Y}] &= E[Y\sin\Theta] = E[Y]E[\sin\Theta]=0\cdot E[\sin\Theta] = 0,\\
E[\hat{X}\hat{Y}] &= E[XY\cos \Theta\sin \Theta] 
= E[X]E[Y]E[\cos \Theta \sin \Theta] = 0\cdot 0\cdot E[\cos \Theta \sin \Theta].
\end{align}$$ 
Note that, conditioned on $\Theta = \theta$, $\hat{X}$ and
$\hat{Y}$are conditionally independent zero-mean
normal random variables with variances $\cos^2\theta$ and
$\sin^2\theta$. They are thus also conditionally uncorrelated
random variables. The conditional joint density is
$$f_{\hat{X},\hat{Y}\mid \Theta}(\hat{x},\hat{y}\mid \Theta=\theta) 
= \frac{1}{2\pi\cos\theta\sin\theta}\exp\left(-\frac{\hat{x}^2}{2\cos^2\theta}
+\frac{\hat{y}^2}{2\sin^2\theta}\right), 
~ -\infty < x, y < \infty.$$
While $\hat{X}$ and $\hat{Y}$ are conditionally independent
given $\Theta = \theta$, their joint density is very much
dependent on $\theta$, and it is not immediately obvious that
their unconditional density
$$f_{\hat{X},\hat{Y}}(\hat{x},\hat{y})
= \int_{-\infty}^\infty f_{\hat{X},\hat{Y}\mid \Theta}(\hat{x},\hat{y}\mid \Theta=\theta)f_\Theta(\theta)\,\mathrm d\theta$$
factors into the product of the marginal densities
$$\begin{align}
f_{\hat{X}}(\hat{x})
&= \int_{-\infty}^\infty f_{\hat{X}\mid \Theta}(\hat{x}\mid \Theta=\theta)f_\Theta(\theta)\,\mathrm d\theta
= \int_{-\infty}^\infty \frac{1}{\sqrt{2\pi}\cos\theta}\exp\left(-\frac{\hat{x}^2}{2\cos^2\theta}
\right)f_\Theta(\theta)\,\mathrm d\theta\\
f_{\hat{Y}}(\hat{y})
&= \int_{-\infty}^\infty f_{\hat{Y}\mid \Theta}(\hat{y}\mid \Theta=\theta)f_\Theta(\theta)\,\mathrm d\theta
= \int_{-\infty}^\infty \frac{1}{\sqrt{2\pi}\sin\theta}\exp\left(-\frac{\hat{y}^2}{2\sin^2\theta}
\right)f_\Theta(\theta)\,\mathrm d\theta
\end{align}$$
as would be needed for $\hat{X}$ and $\hat{Y}$ to be
unconditionally independent random variables.  Nor is it immediately obvious
that the marginal distributions are normal distributions. However, as noted above,
$\hat{X}$ and $\hat{Y}$ are nonetheless unconditionally uncorrelated
random variable.  
Some special cases when $\hat{X}$ and $\hat{Y}$
are indeed unconditionally independent normal random variables
are as follows.


*

*$\Theta$ equals a constant $\theta$ with probability $1$. In this case,
the unconditional density is the same as the conditional density
and thus $\hat{X}$ and $\hat{Y}$ are independent zero-mean
normal random variables with variances $\cos^2\theta$ and $\sin^2\theta$
respectively.

*$\Theta$ is a discrete random variable taking on values
$\pi/4, 3\pi/4, 5\pi/4, 7\pi/4$ with equal probability. In this
case, for all $4$ values of $\Theta$, the conditional joint
distribution is that of independent $N(0,\frac{1}{2})$ random
variables and so the
unconditional distribution also enjoys the same properties.
Note however that the latter choice does not extend to
$\Theta$ being a uniform distribution over $N\geq 5$ phases,
and so it is doubtful that the  $U[0,2\pi)$) for $\Theta$
will yield either normal distributions for $\hat{X}$ and $\hat{Y}$
or that the joint distribution will be the product of the marginals
as required for independence.  The same comments apply to other distributions
that one might choose for $\Theta$. But, $\hat{X}$ and $\hat{Y}$ 
will nonetheless be conditionally as well as unconditionally
uncorrelated random variables regardless of
the distribution of $\Theta$.
Thus, the OP's claim that "$Y1$ and $Y2$ are obviously independent"
(emphasis added) should be taken with a considerably large grain of
salt. If they are indeed independent, there is nothing obvious about
it (to my poor brain, at least; YMMV).
A: Solution:  


*

*The processes $ {Y}_{1} $ and $ {Y}_{2} $ are dependent.
$ {Y}_{1} $ and $ {Y}_{2} $ are dependent intuitively by looking at the problem from a different angle.
Given a uniformly distributed random variable $ \theta $ the random variables $ sin(\theta) $ and $ cos(\theta) $ are clearly dependent. Multiplying them by a Gaussian Random Variable doesn't decorrelate them.

*The processes $ {Y}_{1} $ and $ {Y}_{2} $ aren't Gaussian (And not Ergodic).
Clearly over time (The Phase is constant) the process $ {Y}_{i} $ is Gaussian. Yet the ensemble isn't Gaussian since each of its realization is retrieved from Gaussian Distribution with different parameters.

*The process $ r(t) $ is indeed Gaussian.
This could be proved by using the Characteristic Function of the random process.
$ X \sim N(0, 1) $, $ Y \sim N(0, 1) $, $ \theta \sim U(0, 2\pi) $ all are independent.
Let $ Z = X sin(\theta) + Y cos(\theta) $.  Looking at its Characteristic Function and applying the Smoothing Theorem yields:
\begin{align}
{ \varphi }_{Z}(t) & = & \mathbb{E}[{e}^{itZ}] \\
& = & \mathbb{E} [ \mathbb{E}[{e}^{itZ} | \theta] ] \\
& = & \mathbb{E} [ \mathbb{E}[{e}^{it(X cos(\theta) + Y sin(\theta))} | \theta] ]
\end{align}
Looking at the last equation per realization of $ \theta $:
$$ \mathbb{E}[{e}^{it(X cos(\theta) + Y sin(\theta))} | \Theta = \theta] = \mathbb{E}[{e}^{itX sin(\theta)}] \mathbb{E}[{e}^{itY cos(\theta)}] $$
Each of the item is the Characteristic Function of a scaled Normally Distributed Random Process:
$$ \mathbb{E}[{e}^{itX sin(\theta)}] \mathbb{E}[{e}^{itY cos(\theta)}] = {e}^{-\frac{{t}^{2}}{2} ({sin}^{2}(\theta) + {cos}^{2}(\theta))} = {e}^{-\frac{{t}^{2}}{2}} $$
Hence we get:
$$ { \varphi }_{Z}(t) = \mathbb{E} [ \mathbb{E}[{e}^{itZ} | \Theta = \theta] ] = \mathbb{E} [{e}^{-\frac{{t}^{2}}{2}}] = {e}^{-\frac{{t}^{2}}{2}} $$
Namely, it has the Characteristic Function of a Normalized Gaussian Variable -> $ Z \sim N(0, 1) $.

