\documentclass[10pt, aspectratio=169]{beamer}
\usefonttheme{professionalfonts}
%\usetheme{CambridgeUS}
%
% Choose how your presentation looks.
%
% For more themes, color themes and font themes, see:
% http://deic.uab.es/~iblanes/beamer_gallery/index_by_theme.html
%
\mode<presentation>
{
  \usetheme{default}      % or try Darmstadt, Madrid, Warsaw, ...
  \usecolortheme{beaver} % or try albatross, beaver, crane, ...
  \usefonttheme{default}  % or try serif, structurebold, ...
  \setbeamertemplate{navigation symbols}{}
  \setbeamertemplate{caption}[numbered]
} 

\usepackage[english]{babel}
\usepackage[utf8x]{inputenc}
\usepackage{tikz}
\usepackage{pgfplots}
\usepackage{array}  % for table column M
\usepackage{makecell} % to break line within a cell
\usepackage{verbatim}
\usepackage{graphicx}
\usepackage{epstopdf}
\usepackage{amsfonts}
\usepackage{xcolor}
%\captionsetup{compatibility=false}
%\usepackage{dsfont}
\usepackage[absolute,overlay]{textpos}
\usetikzlibrary{calc}
\usetikzlibrary{pgfplots.fillbetween, backgrounds}
\usetikzlibrary{positioning}
\usetikzlibrary{arrows}
\usetikzlibrary{pgfplots.groupplots}
\usetikzlibrary{arrows.meta}
\usetikzlibrary{plotmarks}
\usetikzlibrary{calc}

\usepgfplotslibrary{groupplots}
\pgfplotsset{compat=newest} 
%\pgfplotsset{plot coordinates/math parser=false}

\usepackage{hyperref}
\hypersetup{
    colorlinks=true,
    linkcolor=blue,
    filecolor=magenta,      
    urlcolor=cyan,
}

%% 
\def\EXTERNALIZE{1} % for externalizing figures
\input{header.tex}

\title[EE 264]{Discrete-Time Random Signals}
\author{Jose Krause Perin}
\institute{Stanford University}
\date{June 27, 2018}

\begin{document}

\begin{frame}
  \titlepage
\end{frame}

\begin{frame}{Last lecture}

\begin{block}{Review of discrete-time signals and systems}
	\begin{itemize}
		\item Systems can be linear, time-invariant, memoryless, causal, and stable \\
		\item Linear and time-invariant (LTI) systems are completely characterized by their impulse response
		\item We can use the convolution sum to calculate the output of an LTI system to any signal \\
		\item The complex exponential $e^{j\omega n}$, and more generally $z^n$, are eigenfunctions of LTI systems
		\item Frequency-domain representation of signals and systems
		\begin{itemize}\normalsize
			\item Discrete-time Fourier transform (DTFT)
			\item $z$-transform and ROC. Without specifying the ROC the $z$-transform is ambiguous
		\end{itemize}
		\item The DTFT is equivalent to the $z$-transform evaluated on the unit circle: $H(e^{j\omega}) = H(z = e^{j\omega})$, provided that the unit circle is in the ROC of the $z$-transform
	\end{itemize}
\end{block}

\end{frame}


\section{Introduction}
\begin{frame}{Today's lecture} 

How to analyze linear and time-invariant (LTI) systems when the input is a random signal.

\begin{block}{Motivation}
	\begin{itemize}
		\item Many signals vary in complicated patterns that
		cannot be easily described analytically
		\item It is often convenient and useful to consider
		such signals as being created by some sort of
		random mechanism
	\end{itemize}   	
\end{block}
\end{frame}

%
\begin{frame}{Example: speech signals}

Speech signals are well described by the \textbf{Laplacian distribution}

\begin{figure}
	\centering
	\resizebox{0.9\linewidth}{!}{\input{figs/speech_and_distribution.tex}}
	\label{fig:speech_and_dist}
\end{figure} 

\end{frame}


%
\begin{frame}{Example: quantization}

Quantization error is well described by an \textbf{uniform distribution}, even though quantization is a deterministic operation 

\begin{center}
	\resizebox{0.7\linewidth}{!}{\input{figs/quantization_example.tex}}
\end{center}

\onslide<5>{More about quantization in lectures 5.}

\end{frame}

%
\begin{frame}{Example: noise and interference}

Noise and interfering signals are typically modeled as \textbf{random processes}

\begin{enumerate}
	\item What's the effect of the noise on the output? 
	\item How can we design the system to minimize the noise at the output?	
\end{enumerate}

\begin{figure}
	\centering
	\resizebox{0.7\linewidth}{!}{\input{figs/system_with_signal_and_noise_inputs.tex}}
\end{figure} 
\end{frame}

\begin{frame}{Outline}

\textbf{Our goal:} analyze LTI systems when the input signal is random

\begin{enumerate}
	\item Random processes
	\begin{itemize}
		\item Averages of a random variable
		\item Stationary random processes
		\item Time averages and ergodic random processes
	\end{itemize}
	\item LTI systems with a random input
	\item White noise
	\item Examples
\end{enumerate}
\end{frame}

%
\section{Random processes}
\begin{frame}{Random processes}

\begin{block}{Definition}
	A random process (or \textit{stochastic process}) is an indexed set of random variables $x_n$, which are distributed according to some probability distribution $p_{x_n}(x)$
\end{block}	

\begin{block}{Examples}
\begin{columns}
	\begin{column}{0.5\textwidth}
		Consecutive coin tosses
	\end{column}
	\begin{column}{0.5\textwidth}  %%<--- here
		Random bit stream
	\end{column}
\end{columns}

\begin{columns}
	\begin{column}{0.5\textwidth}
		\begin{tikzpicture}[draw=black!50, node distance=1cm]
		\tikzstyle{block}=[draw=none,rectangle,fill=none,minimum size=1.5cm, inner sep=0pt]
		\node[block] (C1) {\resizebox{1.5cm}{!}{\includegraphics{figs/US_One_Cent_Obv.png}}};
		\node[block,right of=C1] (C2) {\resizebox{1.5cm}{!}{\includegraphics{figs/US_One_Cent_Rev.png}}};
		\node[block,right of=C2] (C3) {\resizebox{1.5cm}{!}{\includegraphics{figs/US_One_Cent_Rev.png}}};
		\node[block,right of=C3] (C4) {\resizebox{1.5cm}{!}{\includegraphics{figs/US_One_Cent_Obv.png}}};
		\node[block,right of=C4] (C5) {\Large$~~~\cdots$};
		\end{tikzpicture}
	\end{column}
	\begin{column}{0.5\textwidth}  %%<--- here
		\resizebox{0.8\linewidth}{!}{\input{figs/random_bit_sequence.tex}}
	\end{column}
\end{columns}
\end{block}
\end{frame}

\begin{frame}{Random processes}
We can think of a random process as a function $X(n, \chi)$ of two variables, time $n$ and the outcome of the underlying random experiment $\chi$. 
\begin{itemize}
	\pause\item For fixed $n$, $X(n, \chi)$ is a random variable.
	
	In the example of the fair coin tossing, 
	\begin{equation*}
	X(n=1, \chi) = \begin{cases}
	\mathrm{H},~\text{with probability}~0.5 \\
	\mathrm{T},~\text{with probability}~0.5
	\end{cases}
	\end{equation*}
	\pause\item For fixed $\chi$, $X(n, \chi)$ is a deterministic function of $n$ called a \textbf{sample function} or \textbf{sample sequence}
	
	\vspace{3mm}
	\centering
	\begin{tikzpicture}[draw=black!50, node distance=1cm]
		\tikzstyle{block}=[draw=none,rectangle,fill=none,minimum size=1.5cm, inner sep=0pt]
		\node[block] (C1) {\resizebox{1.5cm}{!}{\includegraphics{figs/US_One_Cent_Obv.png}}};
		\node[block,right of=C1] (C2) {\resizebox{1.5cm}{!}{\includegraphics{figs/US_One_Cent_Rev.png}}};
		\node[block,right of=C2] (C3) {\resizebox{1.5cm}{!}{\includegraphics{figs/US_One_Cent_Rev.png}}};
		\node[block,right of=C3] (C4) {\resizebox{1.5cm}{!}{\includegraphics{figs/US_One_Cent_Obv.png}}};
		\node[block,right of=C4] (C5) {\Large$~~~\cdots$};
		\node[below of=C1, scale=0.7] (res1) {$X(1, \chi) = \mathrm{H}$};
		\node[below of=C4, scale=0.7] (res4) {$X(4, \chi) = \mathrm{H}$};
		\node[scale=0.7] at ($(res1.east)!0.5!(res4.west)$) {$\cdots$};
	\end{tikzpicture}
\end{itemize}
\end{frame}

\begin{frame}{Ensemble of sample sequences}
The ensemble of sample sequences is a collection af all sequences generated by a random process.
\vspace{-0.3cm}
\begin{center}
	\resizebox{0.6\linewidth}{!}{\input{figs/ensemble_transitions.tex}}
\end{center}
\vspace{-0.3cm}
\pause[4]
\textbf{Recall:}  
\begin{itemize}
	\item For fixed $n$, we have a random variable (must use probability theory)
	\item For a fixed $\chi$, we have a deterministic signal
\end{itemize}
\end{frame}

\begin{frame}<beamer:0|handout:1>

With this interpretation of random processes, we have two axes:

\begin{itemize}
	\item The probability axis (fixed $n$), whereby the samples are random variables. Hence, we must use probability theory to analyze them.
	\item And we have the sample function axis (fixed $\chi$), whereby we have a deterministic signal generated by a realization of the random process. 
\end{itemize}

~~We'll start by dealing with random variables. We'll see what we can learn about a random process by calculating probability averages such as mean, variance, and autocorrelation function.

~~Then, we'll move to dealing with sample functions, and we'll calculate their time averages.

~~Later we'll see how we can relate probability averages to time averages.

\end{frame}

\subsection{Probability averages}
\begin{frame}{Averages of a random variable}

\textbf{Mean or expected value}
\begin{equation*}
	\mu_{x_n} = \E(x_n) = \int_{-\infty}^{\infty}xp_{x_n}(x)dx
\end{equation*}

\pause
\textbf{Average power or second moment}
\begin{equation*} 
	\E(|x_n|^2) = \int_{-\infty}^{\infty}|x|^2p_{x_n}(x)dx
\end{equation*}

\pause
\textbf{Variance}
\begin{align*}
\sigma^2_{x_n} &= \E(|x_n-\mu_{x_n}|^2) = \int_{-\infty}^{\infty}|x-\mu_{x_n}|^2p_{x_n}(x)dx \\
&= \E(|x_n|^2) - \mu_{x_n}^2
\end{align*}

The integrals should be replaced by sums when the random variable is discrete.

\end{frame}

\begin{frame}{Joint averages of random variables}
Expected value of a function of two random variables
\begin{equation*}
\E(g(x_n,y_m)) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}g(x, y)p_{x_n, y_n}(x, y)dxdy
\end{equation*}

\pause
Two random variables are \textbf{uncorrelated} if
\begin{equation*}
\E(x_ny_m) = \E(x_n)\E(y_m)
\end{equation*}

\pause
Two random variables are \textbf{statistically independent} if
\begin{equation*}
p_{x_n, y_n}(x,y) = p_{x_n}(x)p_{y_n}(y)
\end{equation*}

Independent random variables are also uncorrelated, but not all uncorrelated random variables are independent.

\textbf{Exception:} uncorrelated Gaussian random variables are always independent.
\end{frame}

\begin{frame}{Correlation functions}

\textbf{Autocorrelation}

\begin{equation*}
\phi_{xx}[n, m] = \E(x_nx_m^*) 
\end{equation*}

\textbf{Autocovariance}

\begin{equation*}
\gamma_{xx}[n, m] = \E((x_n-\mu_{x_n})(x_m - \mu_{x_m})^*) 
\end{equation*}

\textbf{Cross-correlation}

\begin{equation*}
\phi_{xy}[n, m] = \E(x_ny_m^*) 
\end{equation*}

\textbf{Cross-covariance}

\begin{equation*}
\gamma_{xy}[n, m] = \E((x_n-\mu_{x_n})(y_m - \mu_{y_m})^*) 
\end{equation*}

\textbf{Note:} By knowing the mean and the autocorrelation function, we can determine average power, variance, and autocovariance.

\end{frame}

\begin{frame}{Example: Bernoulli random process}
	
	\begin{itemize}
		\item A Bernoulli random process is a sequence of binary random variables $\{x_n \sim \mathcal{B}(\rho)\}$. Canonically, 
		\begin{columns}
			\begin{column}{0.5\linewidth}
				\begin{equation*}
				x_n = \begin{cases}
				1,~\text{with probability}~\rho, \\
				0,~\text{with probability}~1-\rho,
				\end{cases}
				\end{equation*}
			\end{column}
			\begin{column}{0.5\linewidth}
				\begin{equation*}
				p_{x_n}(x) = \begin{cases}
				\rho, &x = 1 \\
				1-\rho, &x = 0 \\
				0, &\text{otherwise}
				\end{cases}
				\end{equation*}
			\end{column}
		\end{columns}
		\item A Bernoulli process is \textbf{independent and identically distributed (IID)}. That is, each $x_n$ is picked independently from the same distribution $\mathcal{B}(\rho)$.
	\end{itemize}
	From this we can conclude:
	\begin{align*} 
	\mu &= 1\cdot\rho + 0\cdot(1-\rho) = \rho \\
	\E(x_n^2) &= 1^2\cdot\rho + 0^2\cdot(1-\rho) = \rho \\
	\sigma^2 &= \E(x_n^2) - \mu^2 = \rho(1 - \rho) \\
	\phi_{xx}[m] &=\E(x_{n+m}x_n) = \rho\delta[m] \tag{since it is IID}
	\end{align*}
	
\end{frame}

\begin{frame}{Example: Uniform random process}
	
	\begin{itemize}
		\item An uniform random process is a sequence of uniform random variables $\{x_n \sim \mathcal{U}[a, b]\}$.
		\begin{equation*}
		p_{x_n}(x) = \begin{cases}
		\displaystyle\frac{1}{b-a} & a\leq x\leq b \\
		0, &\text{otherwise}
		\end{cases}
		\end{equation*}
	\end{itemize}
	From this we can conclude:
	\begin{align*} 
	\mu &= \int_{a}^{b} \frac{x}{b-a}dx = \frac{b+a}{2} \\
	\E(x_n^2) &= \int_{a}^{b} \frac{x^2}{b-a}dx = \frac{b^2+ab+a^2}{3} \\
	\sigma^2 &= \E(x_n^2) - \mu^2 = \frac{(b-a)^2}{12} \\
	\phi_{xx}[m] &=\E(x_{n+m}x_n) = \E(x_n^2)\delta[m] \tag{Assuming IID}
	\end{align*}
\end{frame}

\begin{frame}{Example: Gaussian random process}
	
	\begin{itemize}
		\item A Gaussian random process is a sequence of Gaussian random variables $\{x_n \sim \mathcal{N}(\mu, \sigma^2)\}$.
		\begin{equation*}
		p_{x_n}(x) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\bigg(-\frac{(x-\mu)^2}{2\sigma^2}\bigg)
		\end{equation*}

	\end{itemize}
	From this we can conclude:
	\begin{align*} 
	\E(x_n^2) &= \sigma^2 + \mu^2 \\
	\phi_{xx}[m] &=\E(x_{n+m}x_n) = \E(x_n^2)\delta[m] \tag{Assuming IID}
	\end{align*}	
\end{frame}

%%%%%%%%%%%%%%%
\subsection{Stationary random processes}
\begin{frame}<5>{Stationary random processes}
\begin{itemize}
	\item Stationarity refers to \textbf{time invariance} of some, or all, of the statistics of a random process, such as mean, autocorrelation, joint distributions, etc
	\item A random process is \textbf{strict-sense stationary (SSS)}, if its finite-order distributions do not change over time. For the first-order distributions, this means $p_{x_m}(x_m) = p_{x_n}(x_n),~\forall n, m$.
\end{itemize}
 
\centering
\resizebox{0.7\linewidth}{!}{\input{figs/ensemble_transitions.tex}}
\end{frame}

\begin{frame}{Stationary random processes}

All statistics of a SSS random process are time invariant. 

As a result, the mean, average power, and variance are constant with $n$:
\begin{align*}
\mu &= \E(x_n) \\
\sigma^2 &= \E(|x_n|^2) - \mu^2
\end{align*}

And the autocorrelation function only depends on the time difference $m$:
\begin{align*}
\phi_{xx}[m] &= \phi_{xx}[n+m, n] = \E(x_{n+m}x_n^*)
\end{align*}

\textbf{Question:} What is an example of SSS random process?

\end{frame}

\begin{frame}{Wide-sense stationary (WSS) random processes}

\begin{itemize}
	\item Strict sense stationarity is a strong condition that is hard to verify in practice.
	\pause\item A weaker (and more useful) condition is \textbf{wide-sense stationarity}
	\pause\item A random process is \textbf{wide-sense stationary (WSS)} if its mean and autocorrelation function are \textbf{time invariant}. 
	\begin{align*}
	\mu &= \E(x_n) \\
	\phi_{xx}[m] &= \phi_{xx}[n+m, n] = \E(x_{n+m}x_n^*)
	\end{align*}
	The mean is constant, and the autocorrelation function only depends on the time difference $m$.
	\pause\item SSS implies WSS, but WSS does not mean SSS. \\
	\textbf{Exception:} WSS Gaussian random processes are also SSS.	
\end{itemize}
\end{frame}

%
\begin{frame}{Autocorrelation function of WSS processes}

The autocorrelation function $\phi_{xx}[m]$ of a WSS process $x[n]$ has the following properties

\begin{enumerate}
	\item $\phi_{xx}[m]$ is \textbf{real valued}
	\item $\phi_{xx}[m]$ is \textbf{even symmetric}, i.e., $\phi_{xx}[m] = \phi_{xx}[-m]$
	\item The DTFT of $\phi_{xx}[m]$ must be \textbf{non-negative at all frequencies}
	\begin{equation}
	\mathcal{F}\{\phi_{xx}[m]\} \geq 0, \forall~\omega\in[-\pi, \pi]
	\end{equation}
	$\mathcal{F}\{\cdot\}$ denotes the DTFT.	
\end{enumerate}

	Properties 1 to 3 are \textbf{necessary and sufficient} for a function to be an autocorrelation function of a WSS process.
\end{frame}

%
\begin{frame}
More properties
\begin{enumerate}\setcounter{enumi}{3}
	\item $|\phi_{xx}[m]| \leq \phi_{xx}[0] = \E(|x[n]|^2) = \text{\textbf{average power} of}~x[n]$
	
\textit{Proof:}
\begin{align*}
\phi_{xx}^2[m] &= [\E(x[m+n]x[n])]^2 \\
&\leq \E(|x[m+n]|^2)\E(|x[n]|^2) \tag{by Schwarz inequality} \\
&=\phi_{xx}^2[0] \tag{by stationarity}
\end{align*}
	\item If $\phi_{xx}[T] = \phi_{xx}[0]$ for some $T\neq 0$, then $\phi_{xx}[m]$ is periodic with period $T$.
\end{enumerate}
\end{frame}

%
\begin{frame}{Which functions can be $\phi_{xx}[m]$ of a WSS process?}
\begin{enumerate}
	\begin{columns}
		\begin{column}{0.3\linewidth}
			\item 
			\begin{center}
				\resizebox{\linewidth}{!}{\input{figs/right_sided_exp_curve.tex}}
			\end{center}	
		\end{column}
		\begin{column}{0.3\linewidth}
			\item 
			\begin{center}
			\resizebox{\linewidth}{!}{\input{figs/two_sided_exp_curve.tex}}
			\end{center}
		\end{column}
		\begin{column}{0.3\linewidth}
			\item 
			\begin{center}
			\resizebox{\linewidth}{!}{\input{figs/sinc.tex}}
			\end{center}
		\end{column}
	
	\end{columns}
	\begin{columns}
	\begin{column}{0.3\linewidth}
		\item 
		\begin{center}
			\resizebox{\linewidth}{!}{\input{figs/not_autocorrelation_example.tex}}
		\end{center}		
	\end{column}
	\begin{column}{0.3\linewidth}
		\item 
		\begin{center}
		\resizebox{\linewidth}{!}{\input{figs/saw.tex}}
		\end{center}
	\end{column}
		\begin{column}{0.3\linewidth}
		\item 
		\begin{center}
		\resizebox{\linewidth}{!}{\input{figs/constant.tex}}
		\end{center}
	\end{column}
\end{columns}
\end{enumerate}
\end{frame}

\subsection{Time averages and ergodic random processes}

\begin{frame}{Ensemble of sample functions}
\begin{itemize}
	\item So far we have focused on random variables
	\item What can we learn from a sample function?
\end{itemize}
\begin{center}
	\resizebox{0.75\linewidth}{!}{\input{figs/ensemble_stationary.tex}}
\end{center}
\end{frame}

%
\begin{frame}{Time averages}

\begin{itemize}
	\item Until now we have focused on probability averages $\E(\cdot)$
	\item We can also define time averages $\langle\cdot\rangle$
\end{itemize}

\begin{equation*}
\langle x[n] \rangle = \lim_{L\to\infty}\frac{1}{2L + 1}\sum_{n=-L}^L x[n]
\end{equation*}

\begin{equation*}
\langle x[n+m]x^*[n] \rangle = \lim_{L\to\infty}\frac{1}{2L + 1}\sum_{n=-L}^L x[n+m]x^*[n]
\end{equation*}

We can calculate time averages of any deterministic signal

\begin{center}
	\resizebox{0.6\linewidth}{!}{\input{figs/random_sequence.tex}}
\end{center}

\textbf{Careful with notation:} We use $x_n$ to refer to the \textit{random variables} of some random process, whereas $x[n]$ denotes a \textit{sample function} of some random process.

\end{frame}

\begin{frame}{Ergodic random processes}

A random process is \textbf{ergodic} if its time averages are equal to its probability averages:

\begin{align*}
\text{Time averages} &= \text{Probability averages} \\
\langle x[n] \rangle &= \E(x_n) = \mu \tag{expected value} \\
\langle x[n+m]x^*[n] \rangle &= \E(x_{n+m}x_n^*) = \phi_{xx}[m] \tag{autocorrelation function}
\end{align*}

\begin{itemize}
	\pause\item In practice, we don't have an ensemble of sample functions that we can use to estimate the mean and autocorrelation function.
	\pause\item We generally have only one sample function.
	\pause\item With the \textbf{ergodic assumption}, we can estimate probability averages from a single sample function
\end{itemize}

\end{frame}

%
\section{LTI systems with random input}
\begin{frame}{LTI system with a random input}
\begin{center}
\resizebox{\linewidth}{!}{\input{figs/DT_response_to_random.tex}}
\end{center}

\pause
For a given \textbf{sample function} $x[n]$, we can apply the convolution sum as usual:

\begin{equation*}
y[n] = \sum_{m=-\infty}^{\infty} x[n-m]h[m] 
\end{equation*}

\pause\textbf{More importantly:} what is the effect of the system on the statistics (e.g., mean and autocorrelation function) of the random process?

\end{frame}

%
\begin{frame}{LTI system with a random input}

\begin{block}{Mean or expected value}
	\begin{align*}
	\mu_y = \E(y[n]) &= \E\bigg(\sum_{n=-\infty}^{\infty} x[m-n]h[m]\bigg) \\
	&= \sum_{n=-\infty}^{\infty} \E(x[m-n])h[m] \tag{Expectation is a linear operator and $h[n]$ is not random} \\
	&= \mu_x\sum_{n=-\infty}^{\infty} h[m] \tag{Assuming $x[n]$ is WSS} \\
	&= \mu_xH(e^{j0})
	\end{align*}
\end{block}

The mean is scaled by the gain of the LTI system at zero frequency.
\end{frame}

\begin{frame}{LTI system with a random input}

\begin{block}{Autocorrelation function}
	\begin{align} \nonumber
	\phi_{yy}[m] &= \E(y[n+m]y^*[n]) \tag{by definition}  \\ \nonumber
	&= \E\bigg\lbrace\bigg(\sum_{r=-\infty}^{\infty} x[n+m-r]h[r]\bigg)\cdot\bigg(\sum_{k=-\infty}^{\infty} x^*[n-r]h^*[k]\bigg)\bigg\rbrace  \\ \nonumber
	&= \sum_{r=-\infty}^{\infty} h[r] \sum_{k=-\infty}^{\infty} h^*[k]\E(x[n+m-r]x^*[n-k]) \\
	&= \sum_{l=-\infty}^{\infty} \bigg(\sum_{k=-\infty}^{\infty} h[l+k]h^*[k]\bigg)\phi_{xx}[m-l] \tag{change variables $r = l+k$}
	\end{align}
\end{block}
\end{frame}

\begin{frame}
Let's define the \textbf{autocorrelation function of deterministic signals}

\begin{equation*}
c_{hh}[l] \equiv \displaystyle\sum_{k=-\infty}^{\infty} h[l+k]h^*[k]
\end{equation*} 

\pause
Note that the autocorrelation function of deterministic signals and convolution are closely related:
\begin{equation*}
c_{hh}[l] = h[l]\ast h^*[-l]
\end{equation*} 

\pause
Now we can rewrite the autocorrelation function of the output of a LTI system to a random process more compactly:
\begin{align*}
\phi_{yy}[m] &= \sum_{l=-\infty}^{\infty} c_{hh}[l]\phi_{xx}[m-l] \\
&= c_{hh}[m]\ast \phi_{xx}[m] = h[l]\ast h^*[-l] \ast \phi_{xx}[m]
\end{align*} 

\pause
The autocorrelation function of the input random process is \textit{filtered} by the deterministic autocorrelation function of the impulse response.

\end{frame}

\begin{frame}{In the frequency domain}

From the previous derivation in the time domain:
\begin{align*}
\phi_{yy}[m] &= c_{hh}[m]\ast \phi_{xx}[m] \\
&= h[l]\ast h^*[-l]\ast \phi_{xx}[m]
\end{align*}

\pause
In the frequency domain:
\begin{align*}
\mathcal{F}\{\phi_{yy}[m]\} &= H(e^{j\omega})\cdot  H^*(e^{j\omega})\cdot \mathcal{F}\{\phi_{xx}[m]\} \\
&= |H(e^{j\omega})|^2\cdot\mathcal{F}\{\phi_{xx}[m]\}
\end{align*}

\begin{itemize}
	\pause\item The DTFT of the autocorrelation function of a random process is called \textbf{power spectrum density (PSD)}
	\pause\item The PSD has units of W/Hz or dBm/Hz
	\pause\item \textbf{Notation:} $\Phi_{xx}(e^{j\omega}) \equiv \mathcal{F}\{\phi_{xx}[m]\}$
\end{itemize}
\end{frame}

\begin{frame}{Properties of the power spectrum density}

\begin{enumerate} 
	\item The PSD is \textbf{real valued}
	\begin{equation*}
	\Phi_{xx}(e^{j\omega}) = \Phi^*_{xx}(e^{j\omega}),
	\end{equation*}
	since the autocorrelation function has\textbf{ even symmetry}: $\phi_{xx}[m] = \phi_{xx}[-m]$.
	
	\item The PSD is \textbf{even symmetric} 
	\begin{equation*}
	\Phi_{xx}(e^{j\omega}) = \Phi_{xx}(e^{-j\omega}),
	\end{equation*}
	since the autocorrelation function is always real. 
	
	\item The PSD is \textbf{non negative} 
	\begin{equation*}
	\Phi_{xx}(e^{j\omega}) \geq 0,
	\end{equation*}
	This is the same condition required by an autocorrelation function of a WSS random process.
	
	\item The area under $\Phi_{xx}(e^{j\omega})$ is the \textbf{average power}
	\begin{equation*}
	\frac{1}{2\pi}\int_{-\pi}^{\pi} \Phi_{xx}(e^{j\omega})d\omega = \phi_{xx}[0] = E(|x[n]|^2)
	\end{equation*}
	
\end{enumerate} 

\end{frame}

\begin{frame}{Effect of LTI system on the autocorrelation function}
\begin{center}
	\resizebox{0.6\linewidth}{!}{\input{figs/LTI_system_random_input.tex}}
\end{center}

\begin{itemize}
	\item The input autocorrelation function is \textit{filtered} by the LTI system defined by $h[n]\ast h^*[-n] \leftrightarrow |H(e^{j\omega})|^2$ 
	\item The random process remains WSS after an LTI system
\end{itemize}

\end{frame}

\section{White noise}
\begin{frame}{White noise}

\begin{itemize}
	\item White noise is a particularly important class of random process
	\item A \textbf{zero-mean} random process whose \textbf{PSD is constant over all frequencies} is commonly referred to as \textbf{white noise}. Constant PSD over all frequencies implies that the samples are \textbf{uncorrelated} in the time domain
\end{itemize}

\begin{equation*}
\Phi_{xx}(e^{j\omega}) = \sigma_x^2, |\omega|\leq\pi \Longleftrightarrow \phi_{xx}[m] = \sigma_x^2\delta[m]  
\end{equation*}

\begin{center}
	\resizebox{0.6\linewidth}{!}{\input{figs/white_noise.tex}}
\end{center}

\end{frame}

\begin{frame}{White noise into LTI system}

For a white noise input, 
\begin{equation*} \tag{output autocorrelation function}
\phi_{yy}[m] = c_{hh}[m]\ast \phi_{xx}[m] = \sigma_x^2c_{hh}[m]
\end{equation*}

\begin{align*} 
\Phi_{yy}(e^{j\omega}) &= |H(e^{j\omega})|^2\Phi_{xx}(e^{j\omega}) \\
&= \sigma_x^2|H(e^{j\omega})|^2 \tag{output PSD}
\end{align*}

\begin{itemize}
	\item Note that the output noise PSD may not be white. Hence, we say that the filter $H(e^{j\omega})$ \textbf{colored} the noise or \textbf{shaped} the noise.
	\item It is typically easier to analyze systems with white noise. As a result, it is common to employ a \textbf{noise whitening filters} to make the noise white. 
\end{itemize}


\end{frame}

%
\section{Examples}
\begin{frame}{Example: moving average system}

Recall the 4-point moving average system defined by the difference equation:
\begin{equation*}
y[n] = \frac{1}{4}(x[n] + x[n-1] + x[n-2] + x[n-3])
\end{equation*}

This system has impulse response:
\begin{equation*}
h[n] = \frac{1}{4}(\delta[n] + \delta[n-1] + \delta[n-2] + \delta[n-3])
\end{equation*}

And frequency response:
\begin{align*}
H(e^{j\omega}) = \frac{1}{4}(1 + e^{-j\omega} + z^{-j2\omega} + z^{-j3\omega})=\frac{\sin(2\omega)}{4\sin(\omega/2)}e^{-j2\omega}
\end{align*} 

\end{frame}

%
\begin{frame}{White noise into moving average filter}
In \textbf{time domain}:

\begin{center}	
	\resizebox{0.65\linewidth}{!}{
	\begin{tikzpicture}[draw=black!50, node distance=4cm, font=\sffamily]
		\tikzstyle{node}=[circle,fill=black,minimum width=2pt,inner sep=0pt]
		\tikzstyle{block}=[draw=black,rectangle,fill=none,minimum width=3cm, minimum height=1cm, inner sep=0pt]
		
		\node[node] (phix) {};
		\node[block, text width = 4cm, align= center, right of=phix] (DSP2) {$c_{hh}[m] = h[m]\ast h^*[-m]$};
		\coordinate[right of=DSP2] (phiy) {};
		
		\path[->, >=stealth, shorten >= 0pt] (phix) edge (DSP2);
		\path[->, >=stealth, shorten >= 0pt] (DSP2) edge (phiy);
		
		\node[below = 0.5mm of phix, text width=3cm, align=center] { $\phi_{xx}[m] = \sigma_x^2\delta[m]$};	
		\node[below = 0.5mm of phiy] {$\phi_{yy}[m]$};
	\end{tikzpicture}
	}
	
	\resizebox{0.55\linewidth}{!}{\input{figs/moving_average4.tex}}
\end{center}
\end{frame}

%
\begin{frame}{White noise into moving average filter}

In the \textbf{frequency domain}:

\begin{center}	
	\resizebox{0.9\linewidth}{!}{
		\begin{tikzpicture}[draw=black!50, node distance=4cm, font=\sffamily]
		\tikzstyle{node}=[circle,fill=black,minimum width=2pt,inner sep=0pt]
		\tikzstyle{block}=[draw=black,rectangle,fill=none,minimum width=3cm, minimum height=1cm, inner sep=0pt]
		
		\node[node] (phix) {};
		\node[block, text width = 3cm, align= center, right of=phix] (DSP2) {$|H(e^{j\omega})|^2$};
		\coordinate[right of=DSP2] (phiy) {};
		
		\path[->, >=stealth, shorten >= 0pt] (phix) edge (DSP2);
		\path[->, >=stealth, shorten >= 0pt] (DSP2) edge (phiy);
		
		\node[below = 0.5mm of phix, text width=4cm, align=center] { $\Phi_{xx}(e^{j\omega}) = \sigma_x^2$ \\ $\omega\in[-\pi, \pi]$};	
		\node[below = 0.5mm of phiy] {$\Phi_{yy}(e^{j\omega})$};
		\end{tikzpicture}
	}
\end{center}
\vspace{-0.5cm}
\begin{columns}
	\begin{column}{0.5\linewidth}
		\begin{align*}
		H(e^{j\omega}) = \frac{\sin(2\omega)}{4\sin(\omega/2)}e^{-j2\omega}
		\end{align*} 
		The output PSD is therefore
		\begin{align*}
			\Phi_{yy}(e^{j\omega}) &= \sigma_x^2|H(e^{j\omega})|^2 \\
			&= \sigma^2_x\bigg(\frac{\sin(2\omega)}{4\sin(\omega/2)}\bigg)^2
		\end{align*} 
	\end{column}
	\begin{column}{0.5\linewidth}
		\begin{center}
			\resizebox{0.8\linewidth}{!}{\input{figs/psd_moving_average4_to_white_noise.tex}}
		\end{center}
	\end{column}
\end{columns}
\end{frame}

\begin{frame}{Simulation example}

We'll generate an uniform random process and filter it by a 4-point moving average filter: $h[n] = \frac{1}{4}(\delta[n] + \delta[n-1] + \delta[n-2] + \delta[n-3])$

\begin{center}
	\resizebox{0.8\linewidth}{!}{
	\begin{tikzpicture}[draw=black!50, node distance=4cm, font=\sffamily]
		\tikzstyle{node}=[circle,fill=black,minimum size=2pt,inner sep=0pt]
		\tikzstyle{block}=[draw=black,rectangle,fill=none,minimum width=3cm, minimum height=1cm, inner sep=0pt]
		
		\node[node] (xc) {};
		\node[block, right of=xc, text width = 3cm, align= center] (DSP) { $h[n]\leftrightarrow H(e^{j\omega})$};
		\coordinate[right of=DSP] (yc) {};
		
		\path[->, >=stealth, shorten >= 0pt] (xc) edge (DSP);
		\path[->, >=stealth, shorten >= 0pt] (DSP) edge (yc);
		
		\node[block, draw=none] at ($(xc)+(0cm, 0.75cm)$) {\resizebox{7cm}{!}{\input{figs/random_sequence.tex}}};
		\node[below = 0.5mm of xc] {$x[n]$};	
		\node[below = 0.5mm of yc] {$y[n]$};

	\end{tikzpicture}
	}
\end{center}
In Matlab:\\
Generate $1\times 1000$ random vector uniformly distributed in $[-1, 1]$ \\
\texttt{>> x = 2*\textbf{rand}(1, 1000) - 1} \\
Calculate output:\\
\texttt{>> y = \textbf{filter}([1, 1, 1, 1]/4, 1, x)}\\
	
\end{frame}

\begin{frame}{Estimating the autocorrelation function}
\begin{block}{Theoretical autocorrelation} 
	\texttt{>> maxLag = 10 {\color{matlabcomment} \% maximum lag of our autocorrelation function}}\\
	\texttt{>> h = [1, 1, 1, 1]/4 {\color{matlabcomment}\% impulse response}} \\
	\texttt{>> chh = \textbf{conv}(h, conj(fliplr(h))) {\color{matlabcomment}\% deterministic autocorrelation}} \\
	\texttt{>> phi\_xx\_theory = zeros(1, 2*maxLag+1)} \\
	\texttt{>> phi\_xx\_theory(maxLag+1) = 1/3 {\color{matlabcomment}\% theoretical $\phi_{xx}[0]$ for $x_n\sim\mathcal{U}[-1, 1]$}} \\
	\texttt{>> phi\_yy\_theory = \textbf{conv}(phi\_xx\_theory, chh, 'same')} \\
	\texttt{>> \textbf{stem}(-maxLag:maxLag, phi\_xx\_theory) {\color{matlabcomment}\% Plot}} \\
	\texttt{>> \textbf{stem}(-maxLag:maxLag, phi\_yy\_theory)}
\end{block}
\end{frame}

\begin{frame}{Estimating the autocorrelation function}
\begin{block}{Empirical autocorrelation}
	\texttt{>> phi\_yy = \textbf{xcorr}(y, y, maxLag, 'unbiased') {\color{matlabcomment}\% estimate autocorrelation}} \\
	\texttt{>> \textbf{stem}(-maxLag:maxLag, phi\_yy)} \\
\end{block}
\begin{columns}
	\begin{column}{0.5\textwidth}
		\only<2|handout:1>{
		\begin{figure}
			\centering
			\includegraphics[width=\linewidth]{figs/lec2_random_experiment1_phi_xx.eps}
		\end{figure}
	}
	\only<3|handout:2>{
		\begin{figure}
			\centering
			\includegraphics[width=\linewidth]{figs/lec2_random_experiment2_phi_xx.eps}
		\end{figure}
	}
	\end{column}

	\begin{column}{0.5\textwidth}
		\only<2|handout:1>{
		\begin{figure}
			\centering
			\includegraphics[width=\linewidth]{figs/lec2_random_experiment1_phi_yy.eps}
		\end{figure}
		}
		\only<3|handout:2>{
			\begin{figure}
				\centering
				\includegraphics[width=\linewidth]{figs/lec2_random_experiment2_phi_yy.eps}
			\end{figure}
		}
	\end{column}
\end{columns}
\onslide<3|handout:2>{Now the random vector is $10000\times 1$}
\end{frame}

\begin{frame}{Estimating the probability density function}
\texttt{>> Nbins = 20} \\
\texttt{>> [counts, centers] = \textbf{hist}(x, Nbins)} \\
{\color{matlabcomment}\% normalize to make area under pdf equal to 1} \\
\texttt{>> x\_pdf = Nbins/(centers(end) - centers(1))*counts/sum(counts)} \\
\texttt{>> bar(centers, x\_pdf) {\color{matlabcomment}\% plot}}
\begin{columns}
	\begin{column}{0.5\textwidth}
			\begin{figure}
				\centering
				$p_x(x)$ \\
				\includegraphics[width=0.8\linewidth]{figs/lec2_random_experiment1_hist_x.eps}
			\end{figure}		
	\end{column}
	\begin{column}{0.5\textwidth}
	\begin{figure}
		\centering
		$p_y(y)$ \\
		\includegraphics[width=0.8\linewidth]{figs/lec2_random_experiment1_hist_y.eps}
	\end{figure}		
\end{column}
\end{columns}
\end{frame}

\begin{frame}{Central limit theorem}
The \textbf{central limit theorem} states that the probability density function of the sum of a large number of independent random variables approaches a Gaussian distribution.
\begin{equation*}
Z = X_1 + X_2 + \ldots + X_N \implies p_Z(z) \xrightarrow{N\to\infty} \mathcal{N}(\mu, \sigma^2)
\end{equation*}

\begin{itemize}
	\item In digital filters, we're essentially performing a weighted sum of samples of the input. As a result, the output of a digital filter for a random input is approximately Gaussian distributed.
	\item In the example of the 4-point moving average system, we were only adding 4 random variables and the resulting output pdf was close to Gaussian. 
	\item This effect is more visible in filters with larger memory, and in particular in filters with \textbf{infinite impulse response (IIR)} 
\end{itemize}

\end{frame}

\begin{frame}{Another probability theorem}

If we add two independent random variables $Z = X + Y$, the pdf of $Z$ is given by the convolution of the pdfs of $X$ and $Y$

\begin{equation*}
p_Z = p_X \ast p_Y
\end{equation*}

\pause
We can extend that to a sum of several independent random variables:
\begin{equation*}
Z = \sum_{k=1}^N X_k \implies p_Z = p_{X_1} \ast \ldots \ast p_{X_N}
\end{equation*}

\pause
From the \textbf{central limit theorem}, we know that as $N\to\infty$,  $p_Z\to\mathcal{N}(\mu, \sigma^2)$

\begin{itemize}
	\item This shows that the convolution of many signals $h_1(t) \ast h_2(t) \ast \ldots \ast h_N(t) \approx g(t)$, where $g(t)$ is the Gaussian function.
	\item Hence, if we cascade many LTI systems, the impulse response of the equivalent LTI system is approximately the Gaussian function.
\end{itemize}


\end{frame}

\begin{frame}{Summary}
\begin{itemize}
	%\item We use random processes to model signals that cannot be easily described by simple equations
	\item A random process is an indexed collection of random variables
	\item A random process is strict-sense stationary (SSS) if all its finite-order statistics are time invariant. That's hard to verify in practice.
	\item A random process is wide-sense stationary (WSS) if its mean is constant and if its autocorrelation function is only a function of the time difference. 
	\item A random process is ergodic if its time averages are equal to its probability averages
	\item The Fourier transform of the autocorrelation function is called the power spectrum density (PSD). The PSD has units of W/Hz or dBm/Hz.
	\item When a random signal is filtered by an LTI system defined by $h[n]\leftrightarrow H(e^{j\omega})$, its autocorrelation function is filtered by an LTI system defined by $h[n]\ast h^*[-n]$, and its PSD is shaped by $|H(e^{j\omega})|^2$
	\item Random processes that have PSD constant over all frequencies are called white noise
	\item By the central limit theorem, the output of an LTI system to a random input is approximately Gaussian distributed
	
\end{itemize}
\end{frame}

\end{document}
