\documentclass[12pt]{article}
\usepackage{graphicx}
\usepackage{float}
\usepackage{hyperref}
\usepackage{multirow}
\usepackage{amsmath}
\usepackage{enumitem}
\usepackage{subcaption}
\usepackage{url}
\usepackage{array}

\usepackage{listings}
\lstset{%
%backgroundcolor=\color{cyan!10},
basicstyle=\ttfamily,
numbers=left,numberstyle=\scriptsize
}

% adjust margins
\addtolength{\oddsidemargin}{-.5in}
\addtolength{\evensidemargin}{-.5in}
\addtolength{\textwidth}{1in}

\addtolength{\topmargin}{-.5in}
\addtolength{\textheight}{1in}

\renewcommand{\arraystretch}{2}


\title{Fourier Series}
\author{EE 224: Signals and Systems I}
\date{}

\begin{document}

\maketitle

\section{Overview}

In this lab you will experience the power of Fourier series analysis and synthesis.

\section{Learning Objectives}

By the end of this lab, students will be able to:
\begin{enumerate}
\item Use MATLAB to approximate Fourier series coefficients of a sampled signal.
\item Use MATLAB to synthesize a signal from its Fourier series coefficients.
\item Describe the difference in sound made by two musical instruments in terms of their Fourier series coefficients.
\item Create your own synthesized instrument using Fourier series (bonus).
\end{enumerate}

\section{Pre-Lab Reading}

\subsection{Fourier Analysis}
Begin by reviewing the Fourier analysis formula. Given a periodic signal $x(t)$ with period $T$, $x(t)$ can be
represented by
\begin{equation}
x(t)=\sum_{k=-\infty}^\infty a_k e^{j k \omega_0 t}
\end{equation}
where $\omega_0=\frac{2\pi}T$ and
\begin{equation}
a_k=\frac{1}{T}\int_{0}^T x(t) e^{-jk \omega_0 t} dt\qquad \mbox{ for all } k
\end{equation}
Although an infinite number of harmonics may be required for a general
signal, in most situations, a finite number of them provides a practically
good approximation.

If the signal $x(t)$ is not given as a mathematical function but we
have a sampled recorded version of it, the integrals to compute the
coefficients ($a_k$) cannot be evaluated precisely. Although there are more
efficient methods to perform the Fourier analysis directly on the sampled
signals (e.g., the Fast Fourier Transform (FFT)), we will use a simple
method to approximately evaluate those integrals: the Riemann
approximation. Assuming that the signal $x(t)$ is reasonably smooth (Riemann
integrable) and that the sampling time $T_s$ is an integer fraction of the
period $T$, i.e., $T = N T_S$ of samples in a period
(this is not really necessary but simplifies our code), then:
\begin{equation}\label{eq1.eq}
a_k=\frac{1}{T}\int_{0}^T x(t) e^{-jk \omega_0 t} dt\simeq
	\frac{T_s}{T}\sum_{n=0}^{T/T_s} x[T_s n]e^{-j k \omega_0 T_s n}.
\end{equation}
The approximation gets better and better as the number of samples goes to
infinity (meaning $T_s$ goes to zero).  In other words, fixing the sampling time limits the quality of the
approximation, especially for larger values of $k$. This is because $T_s$ becomes
too large with respect to $k\omega_0$, the frequency of the $k$th harmonic.

Find the Fourier series coefficients $a_0$, $a_1$, $a_2$, and $a_3$ for the square wave given below (assume the signal repeats with period $T=0.04$).  Use the integral in \eqref{eq1.eq} and include your answer in your report. %enter your answer into the Canvas prelab.
\begin{align*}
x(t) = \begin{cases} 1 & 0 \leq t < 0.02 \\
0 & 0.02 \leq t < 0.04
\end{cases}
\end{align*}

\subsection{Fourier Synthesis}
The synthesis formula allows one to generate periodic signals from the linear
combination of harmonic complex exponentials. Here we assume that the
number of non-zero coefficients $a_k$ is finite. Then $$
\tilde{x}(t)=\sum_{k}a_ke^{j k \omega_0 t} $$ In particular, $$
\tilde{x}{}[T_s n]=\sum_{k}a_ke^{j k \omega_0 T_s n} $$ 


\section{Lab Exercises}

\subsection{Fourier Series Analysis Function}
\label{sec:fanal}

Write a Matlab function ``findFS" which takes as inputs $x_T$, a vector of
samples of $x(t)$ covering one period of $x(t)$; $k$, the index of the desired
coefficient; the period $T$; and the sampling time $T_s$; and produces
as output the approximate $a_k$ coefficient according to (\ref{eq1.eq}).
Note that you can easily take advantage of vectorization by considering the
product of the row vector $x_T$ and the column vector $e^{-j k \omega_0
T_s n}$. A template for the function is provided.

%\vspace{\baselineskip}
%\noindent
%\emph{Get verification from a TA here.}


\subsection{Fourier Series Coefficients of a Square Wave}
\label{sec:square}

 Let's construct a test signal.
Define in Matlab \\
\centerline{\texttt{x\_T=[ones(1000,1);zeros(1000,1)];} }

This is one period of a square wave signal. Let the fundamental period be
$T=0.04$ sec.
\begin{itemize}
\item[a)] Calculate the sampling time $T_s$ and explain your answer.  (Hint: $T_s=T/2000$)
\item[b)] Use the function you wrote in the previous exercise to compute
$a_0$, $a_1$, $a_2$, and $a_3$.
\item[c)] Find $a_{-1}$, $a_{-2}$, and $a_{-3}$.  How do these relate to $a_1$, $a_2$, and $a_3$?
\item[d)] Verify the coefficients you have computed approximate reasonably
well the actual coefficients you obtained for the prelab.
\end{itemize}

%\vspace{\baselineskip}
%\noindent
%\emph{Get verification from a TA here.}

\subsection{Fourier Series Synthesis}
\label{sec:fsynt}

Write a Matlab function ``synthFS'' which has the following inputs: a column
vector $C=[a_0; a_1; \ldots; a_m]$ of coefficients where $m$ is the largest
integer corresponding to a non-zero coefficient, the sampling time $T_s$,
and the fundamental period $T$.  The output should be $\tilde{x}_T(T_s n)$, the vector of
$N = T/T_s$ ($N$ an integer) samples corresponding to one period of $\tilde{x}$. Note
that your program should verify that  $T/T_s$ is an integer.
You will also need to extend $C$ to include the complex conjugate
coefficients corresponding to the negative $k$s. Finally, you may want to
generate an array \texttt{F} whose columns are the vectors $e^{j k \omega_0
T_s n}$. $\tilde{x}_T$ is then computed as \texttt{F*CC}, where \texttt{CC} is a
vector containing all the coefficients from $-m$ to $m$. A template for this function is provided.

%\vspace{\baselineskip}
%\noindent
%\emph{Get verification from a TA here.}


\subsection{Spectra of a Trumpet and Whistle}
\label{sec:spectra}

Load the data file \texttt{trumpet\_whistle.mat}, which can be found on Canvas. It
contains vectors \texttt{trumpet} and \texttt{whistle}. Both signals are sampled at $f_0=44100$
Hz, and they represent one period of synthetic trumpet and whistle tones
generated by a toy electric keyboard.

\begin{itemize}
\item[a)] Compute the period of the two signals.
\item[b)] For each signal, compute and report the first 9 harmonic
coefficients, i.e., $a_k$s with $k=-9,\ldots,9$ using the function
``findFS'' you have developed.
\item[c)] For each signal, plot the magnitude and phase spectra of the frequency components from part b.
Briefly comment on their main differences. You may use the Matlab
function \texttt{stem} to plot the spectrum. 
\item[d)] For each signal, synthesize an approximation by using the periods
in part a, the coefficients in part b, the sampling time $T_s=1/f_0$, and
your function ``synthFS.''
\item[e)] For each signal, plot on the same plot the given signal and its
synthesized approximation. (You may use \texttt{spectraplot.m} found on Canvas. Make sure the axis labels are accurate.) Comment on the quality of the approximation in
both cases.
\item[f)] For each signal and its approximation, generate a new signal by
repeating them for 1000 periods (e.g., using the built in function \texttt{repmat}). Use the function \texttt{sound} in Matlab and
the sampling frequency $f_0$ to hear the sound you have generated. Can you
hear the difference between the originals and the approximations? Comment.
(Note you may need to scale the signals to have magnitude 1, see \texttt{help
sound}.)
\end{itemize}
Note that the trumpet sound has a much richer spectrum and corresponds to a
more complex sound. While the whistle sound essentially does not have any
harmonic components above the ninth and does not have any even harmonic
components, the trumpet sound has important harmonics above the ninth;
this can be argued from the error in the approximation.

%\vspace{\baselineskip}
%\noindent
%\emph{Get verification from a TA here.}

\subsection{Bonus}
\label{sec:bonus}

Synthesize your own ``instrument'' by repeating steps d-f of the last subsection.  Instead of using the $a_k$ that you calculated for the trumpet and whistle, come up with your own coefficients any way you want.  (You don't have to limit yourself to nine coefficients.) Upload your sound and describe it in your report.

%\vspace{\baselineskip}
%\noindent
%\emph{Demo your sound to a TA.}

\subsection{Report Checklist}

Be sure the following are included in your report.
\begin{enumerate}
\item Section~\ref{sec:fanal}: include your function in the appendix
\item Section~\ref{sec:square}: provide answers for a-d
\item Section~\ref{sec:fsynt}: include your function in the appendix
\item Section~\ref{sec:spectra}: provide answers for a-c, e, and f
\item Section~\ref{sec:bonus}: (bonus) describe the sound you created and explain how you created it
\end{enumerate}

%\section{References}

%If necessary

%\begin{thebibliography}{9}


%\end{thebibliography}





\end{document} 