Page MenuHomec4science
Files F63356418

fourier.tex
No OneTemporary
Actions

File Metadata

Created
Sun, May 19, 13:11

fourier.tex
View Options

\documentclass{article}
% Change "article" to "report" to get rid of page number on title page
\usepackage{amsmath,amsfonts,amsthm,amssymb,url}
\usepackage{setspace}
\usepackage{fancyhdr}
\usepackage{lastpage}
\usepackage{extramarks}
\usepackage{chngpage}
\usepackage{soul}
\usepackage[usenames,dvipsnames]{color}
\usepackage{graphicx,float,wrapfig}
\usepackage{ifthen}
\usepackage{listings}
\usepackage{courier}
\usepackage{mathtools}
\usepackage{gensymb}
\usepackage{subcaption}
% units of mesure
\usepackage{siunitx}
% comments
\usepackage{comment}
% pgfplots environment
%package setup for graphs
\usepackage{tikz}
\usepackage{pgfplots}
\usepackage{pgfplotstable}
\usetikzlibrary{patterns}
\usepgfplotslibrary{external}
\pgfplotsset{compat=newest}
%pgfplot setup
\makeatletter
\pgfplotsset{
/pgfplots/flexible xticklabels from table/.code n args={3}{%
\pgfplotstableread[#3]{#1}\coordinate@table
\pgfplotstablegetcolumn{#2}\of{\coordinate@table}\to\pgfplots@xticklabels
\let\pgfplots@xticklabel=\pgfplots@user@ticklabel@list@x
},
% layer definition
layers/my layer set/.define layer set={
background,
main,
foreground
}{
% you could state styles here which should be moved to
% corresponding layers, but that is not necessary here.
% That is why we don't state anything here
},
% activate the newly created layer set
set layers=my layer set
}
\makeatother
\IfFileExists{./tikzext.perm}{\tikzexternalize[prefix=tikzext/]}{\tikzexternalize}
% definition of the absolute value
\DeclarePairedDelimiter\abs{\lvert}{\rvert}
\usepackage{scalerel,stackengine}
\stackMath
\newcommand\reallywidehat[1]{%
\savestack{\tmpbox}{\stretchto{%
\scaleto{%
\scalerel*[\widthof{\ensuremath{#1}}]{\kern-.6pt\bigwedge\kern-.6pt}%
{\rule[-\textheight/2]{1ex}{\textheight}}%WIDTH-LIMITED BIG WEDGE
}{\textheight}%
}{0.5ex}}%
\stackon[1pt]{#1}{\tmpbox}%
}
% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%
% Here put your info (name, due date, title etc).
% the rest should be left unchanged.
%
%
%
% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% Homework Specific Information
\newcommand{\hmwkTitle}{Report 1}
\newcommand{\hmwkSubTitle}{Fourier transforms and analysis}
\newcommand{\hmwkDueDate}{March 27, 2014}
\newcommand{\hmwkClass}{Computational Physics III}
\newcommand{\hmwkClassTime}{\today}
%\newcommand{\hmwkClassInstructor}{Prof. Oleg Yazyev}
\newcommand{\hmwkAuthorName}{Raffaele Ancarola}
%
%
% In case you need to adjust margins:
\topmargin=-0.45in %
\evensidemargin=0in %
\oddsidemargin=0in %
\textwidth=6.5in %
\textheight=9.5in %
\headsep=0.25in %
% This is the color used for comments below
\definecolor{MyDarkGreen}{rgb}{0.0,0.4,0.0}
% For faster processing, load Matlab syntax for listings
\lstloadlanguages{Matlab}%
\lstset{language=Matlab, % Use MATLAB
frame=single, % Single frame around code
basicstyle=\small\ttfamily, % Use small true type font
keywordstyle=[1]\color{Blue}\bf, % MATLAB functions bold and blue
keywordstyle=[2]\color{Purple}, % MATLAB function arguments purple
keywordstyle=[3]\color{Blue}\underbar, % User functions underlined and blue
identifierstyle=, % Nothing special about identifiers
% Comments small dark green courier
commentstyle=\usefont{T1}{pcr}{m}{sl}\color{MyDarkGreen}\small,
stringstyle=\color{Purple}, % Strings are purple
showstringspaces=false, % Don't put marks in string spaces
tabsize=3, % 5 spaces per tab
breaklines=True,
linewidth=\textwidth,
%
%%% Put standard MATLAB functions not included in the default
%%% language here
morekeywords={xlim,ylim,var,alpha,factorial,poissrnd,normpdf,normcdf},
%
%%% Put MATLAB function parameters here
morekeywords=[2]{on, off, interp},
%
%%% Put user defined functions here
morekeywords=[3]{FindESS, homework_example},
%
morecomment=[l][\color{Blue}]{...}, % Line continuation (...) like blue comment
numbers=left, % Line numbers on left
firstnumber=1, % Line numbers start with line 1
numberstyle=\tiny\color{Blue}, % Line numbers are blue
stepnumber=1 % Line numbers go in steps of 5
}
% Setup the header and footer
\pagestyle{fancy} %
\lhead{\hmwkAuthorName} %
%\chead{\hmwkClass\ (\hmwkClassInstructor\ \hmwkClassTime): \hmwkTitle} %
\rhead{\hmwkClass\ : \hmwkTitle} %
%\rhead{\firstxmark} %
\lfoot{\lastxmark} %
\cfoot{} %
\rfoot{Page\ \thepage\ of\ \protect\pageref{LastPage}} %
\renewcommand\headrulewidth{0.4pt} %
\renewcommand\footrulewidth{0.4pt} %
% This is used to trace down (pin point) problems
% in latexing a document:
%\tracingall
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some tools
\newcommand{\enterProblemHeader}[1]{\nobreak\extramarks{#1}{#1 continued on next page\ldots}\nobreak%
\nobreak\extramarks{#1 (continued)}{#1 continued on next page\ldots}\nobreak}%
\newcommand{\exitProblemHeader}[1]{\nobreak\extramarks{#1 (continued)}{#1 continued on next page\ldots}\nobreak%
\nobreak\extramarks{#1}{}\nobreak}%
\newlength{\labelLength}
\newcommand{\labelAnswer}[2]
{\settowidth{\labelLength}{#1}%
\addtolength{\labelLength}{0.25in}%
\changetext{}{-\labelLength}{}{}{}%
\noindent\fbox{\begin{minipage}[c]{\columnwidth}#2\end{minipage}}%
\marginpar{\fbox{#1}}%
% We put the blank space above in order to make sure this
% \marginpar gets correctly placed.
\changetext{}{+\labelLength}{}{}{}}%
\setcounter{secnumdepth}{0}
\newcommand{\homeworkProblemName}{}%
\newcounter{homeworkProblemCounter}%
\newenvironment{homeworkProblem}[1][Problem \arabic{homeworkProblemCounter}]%
{\stepcounter{homeworkProblemCounter}%
\renewcommand{\homeworkProblemName}{#1}%
\section{\homeworkProblemName}%
\enterProblemHeader{\homeworkProblemName}}%
{\exitProblemHeader{\homeworkProblemName}}%
\newcommand{\problemAnswer}[1]
{\noindent\fbox{\begin{minipage}[c]{\columnwidth}#1\end{minipage}}}%
\newcommand{\problemLAnswer}[1]
{\labelAnswer{\homeworkProblemName}{#1}}
\newcommand{\homeworkSectionName}{}%
\newlength{\homeworkSectionLabelLength}{}%
\newenvironment{homeworkSection}[1]%
{% We put this space here to make sure we're not connected to the above.
% Otherwise the changetext can do funny things to the other margin
\renewcommand{\homeworkSectionName}{#1}%
%\settowidth{\homeworkSectionLabelLength}{\homeworkSectionName}%
\addtolength{\homeworkSectionLabelLength}{0.25in}%
\changetext{}{-\homeworkSectionLabelLength}{}{}{}%
\subsection{\homeworkSectionName}%
\enterProblemHeader{\homeworkProblemName\ [\homeworkSectionName]}}%
{\enterProblemHeader{\homeworkProblemName}%
% We put the blank space above in order to make sure this margin
% change doesn't happen too soon (otherwise \sectionAnswer's can
% get ugly about their \marginpar placement.
\changetext{}{+\homeworkSectionLabelLength}{}{}{}}%
\newcommand{\sectionAnswer}[1]
{% We put this space here to make sure we're disconnected from the previous
% passage
\noindent\fbox{\begin{minipage}[c]{\columnwidth}#1\end{minipage}}%
\enterProblemHeader{\homeworkProblemName}\exitProblemHeader{\homeworkProblemName}%
\marginpar{\fbox{\homeworkSectionName}}%
% We put the blank space above in order to make sure this
% \marginpar gets correctly placed.
}%
\newcommand{\unlimwrap}[3][0.45\linewidth]
{
\begin{minipage}{\linewidth}
\hspace{-0.5cm}
\begin{minipage}{#1}
#2
\end{minipage}
\hspace{1cm}
\begin{minipage}{0.7\linewidth}
#3
\end{minipage}
\end{minipage}
}
%%% I think \captionwidth (commented out below) can go away
%%%
%% Edits the caption width
%\newcommand{\captionwidth}[1]{%
% \dimen0=\columnwidth \advance\dimen0 by-#1\relax
% \divide\dimen0 by2
% \advance\leftskip by\dimen0
% \advance\rightskip by\dimen0
%}
% Includes a figure
% The first parameter is the label, which is also the name of the figure
% with or without the extension (e.g., .eps, .fig, .png, .gif, etc.)
% IF NO EXTENSION IS GIVEN, LaTeX will look for the most appropriate one.
% This means that if a DVI (or PS) is being produced, it will look for
% an eps. If a PDF is being produced, it will look for nearly anything
% else (gif, jpg, png, et cetera). Because of this, when I generate figures
% I typically generate an eps and a png to allow me the most flexibility
% when rendering my document.
% The second parameter is the width of the figure normalized to column width
% (e.g. 0.5 for half a column, 0.75 for 75% of the column)
% The third parameter is the caption.
\newcommand{\scalefig}[3]{
\begin{figure}[ht!]
% Requires \usepackage{graphicx}
\centering
\includegraphics[width=#2\columnwidth]{#1}
%%% I think \captionwidth (see above) can go away as long as
%%% \centering is above
%\captionwidth{#2\columnwidth}%
\caption{#3}
\label{#1}
\end{figure}}
% Includes a MATLAB script.
% The first parameter is the label, which also is the name of the script
% without the .m.
% The second parameter is the optional caption.
%\newcommand{\matlabscript}[2]
% {\begin{itemize}\item[]\lstinputlisting[caption=#2,label=#1]{#1.m}\end{itemize}}
\newcommand{\matlabscript}[2]
{\lstinputlisting[caption=#2,label=#1]{#1.m}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make title
%\title{\vspace{2in}\textmd{\textbf{\hmwkClass:\ \hmwkTitle\ifthenelse{\equal{\hmwkSubTitle}{}}{}{\\\hmwkSubTitle}}}\\\normalsize\vspace{0.1in}\small{Due\ on\ \hmwkDueDate}\\\vspace{0.1in}\large{\textit{\hmwkClassInstructor\ \hmwkClassTime}}\vspace{3in}}
\title{\vspace{2in}\textmd{\textbf{\hmwkClass:\ \hmwkTitle\ifthenelse{\equal{\hmwkSubTitle}{}}{}{\\\hmwkSubTitle}}}\\\normalsize\vspace{0.1in}\small{Due\ on\ \hmwkDueDate}\\\vspace{0.1in}\large{\textit{ \hmwkClassTime}}\vspace{3in}}
\date{}
\author{\textbf{\hmwkAuthorName}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\begin{spacing}{1.1}
\maketitle
% Uncomment the \tableofcontents and \newpage lines to get a Contents page
% Uncomment the \setcounter line as well if you do NOT want subsections
% listed in Contents
%\setcounter{tocdepth}{1}
\newpage
\tableofcontents
\newpage
\section{Introduction: theoretical basis}
\subsection{The fourier transform}
\paragraph{Definition in $L^2$}
Let $f \in L^2(\mathbb{R})$, then it's \textit{Fourier transform} $\hat{f}$ is defined as
\begin{align} \label{defnt}
\hat{f}(k) := \mathcal{F}\{f\}(k) := \int_{-\infty}^{\infty} f(x) \cdot e^{-2 \pi k x} dx & &
f(x) := \mathcal{F}^{-1}\{\hat{f}\}(x) := \int_{-\infty}^{\infty} \hat{f}(k) \cdot e^{2 \pi x k} dk
\end{align}
where the second definition is called the \textit{inverse Fourier transform} and
the bijection is a consequence of $L^2(\mathbb{R})$ [\cite{stubbe}].
\paragraph{Properties}
The following properties can be straigh-forward demonstrated:
\begin{itemize}
\item Linearity: $\forall f,g \in L^2, a \in \mathbb{R}: \reallywidehat{f + a \cdot g} = \hat{f} + a \cdot \hat{g}$
\item Scaling: $\forall a \in \mathbb{R}: \reallywidehat{f(a \cdot x)} = \frac{1}{\abs*{a}} \hat{f}(\frac{k}{a})$
\item Translation: $\forall x_0 \in \mathbb{R}: \reallywidehat{f(x - x_0)} = e^{-2 \pi i x_0 k} \hat{f}(k)$
\item Modulation: $\forall k_0 \in \mathbb{R}: \reallywidehat{e^{2 \pi i x k_0} f(x)} = \hat{f}(k - k_0)$
\item Parseval identity: $\forall f,g \in L^2: \langle f, g \rangle = \langle \hat{f}, \hat{g} \rangle$ where $\langle \cdot , \cdot \rangle$ is the canonical $L^2$ scalar product (see [\cite{stubbe}])
\end{itemize}
\paragraph{Convolution definition and theorem}
Let $f, g \in L^2(\mathbb{R})$, then the \textit{convolution} of $f$ with respect to $g$ is defined as follow:
\begin{equation} \label{convolution}
(f \circledast g)(x) := \int_{-\infty}^{\infty} f(y)g(x - y) dy
\end{equation}
Having listed enough properties of the \textit{Fourier transform}, it's now possible to relate these two concepts by stating the \textit{convolution} theorem:
\newline
\begin{equation} \label{conv_thm}
\reallywidehat{f \circledast g} = \hat{f} \cdot \hat{g}
\end{equation}
This result is very useful in its practical applications, because it allows to convert
a convolution into a product operation.
\newpage
\section{Discrete Fourier transform}
The passage from the continous to a discrete expression of the Fourier transform is
necessary in order to process digital signals. The interest on doing such an elaboration
is that an analogic circuit which computes the continous
Fourier transform on a generic signal doesn't forcely exist.
\subsection{Discretization}
In order to pass to a discrete representation of the Fourier transform, some conventions need to be set:
\begin{itemize}
\item The signal is evaluated on a limited set of time arguments.
\item The first time of the sequence is $0$
\item The function is supposed to be constant between one sample and another.
\item The function is zero in all other points.
\end{itemize}
Let $\{f_n\}_{0 < n < N}$ be the signal samples derived from a function $f(x)$, where
$N$ is the total number of samples and $T$ the sampling period.
Because $\forall x \in \mathbb{R} \setminus [0, (N-1)T]: f(x) = 0$ the fourier integration domain
restricts to $[0, (N-1)T]$, then applying a discretisation on $x$, given by $x_n = nT, 0 < n < N$,
the integral reduces itself into a sum:
\begin{equation} \label{dft_deriv}
\hat{f}(k) = \sum_{n=0}^{N-1} \int_{nT}^{(n+1)T} f_n \cdot e^{-2\pi i nT k} dx = T \sum_{n=0}^{N-1} f_n \cdot e^{-2\pi i nT k}
\end{equation}
Then taking the indicisation of the $k$-space as $k_m = \frac{m}{NT}, 0 < m < N$, the
output results a $N$-sized set of values $\{\hat{f}_m\}_{0 < m < N}$:
\begin{equation} \label{dft_fake}
\hat{f}_m = T \sum_{n=0}^{N-1} f_n \cdot e^{-2\pi i \frac{nm}{N}}
\end{equation}
Notice that this expression of the discrete Fourier transform (or DFT) still depends on $T$ and since it's a scaling term,
it's presence is arbitrary and it could be removed.
In order to guaratee the application bijectivity, it's important evaluate the
inverse discrete fourier transform by taking the adjoint of the sum given in equation (\ref{dft_fake})
on the $\{\hat{f}_m\}_{0 < m < N}$ set and setting it to $\{f_n\}_{0 < n < N}$.
\begin{equation} \label{inverse_dft}
f_n = \mathcal{F}^{-1}\{\hat{f}\}_n = \frac{T}{NT} \sum_{m=0}^{N-1} \hat{f}_m \cdot e^{2\pi i \frac{nm}{N}} = \frac{1}{N} \sum_{m=0}^{N-1} \hat{f}_m \cdot e^{2\pi i \frac{nm}{N}}
\end{equation}
The previous expression shows that the inverse DFT doesn't depend on the sampling period $T$, then
there's no reason to keep that dependency for the DFT. This consideration allows to redefine the DFT as follow:
\begin{equation} \label{dft}
\hat{f}_m := \mathcal{F}\{f\}_m = \sum_{n=0}^{N-1} f_n \cdot e^{-2\pi i \frac{nm}{N}}
\end{equation}
Notice then that the DFT is periodic in $N$:
\begin{equation} \label{dft_period}
\hat{f}_{m+N} = \sum_{n=0}^{N-1} f_n \cdot e^{-2\pi i \frac{n(m+N)}{N}} = \sum_{n=0}^{N-1} f_n \cdot e^{-2\pi i \frac{nm}{N}} \cdot e^{-2\pi i n} = \hat{f}_m
\end{equation}
\subsection{Cepstrum and echo} \label{cepstrum_math}
The \textit{cepstrum} is a numerical instrument that allows to detect echo in a signal.
In principle, a echoed signal is supposed to be composed of various convolutions of functions:
\begin{equation} \label{echoed_signal}
f(x) = (f_1 \circledast f_2 \circledast ... f_N)(x)
\end{equation}
where $N$ in this case is the total number of functions composing the signal.
Using the convolution theorem (equation (\ref{conv_thm})), by induction it's easy to demonstrate that
the fourier transformed signal $\hat{f}(k)$ becomes the product of all fourier transformed components $\hat{f_n}$.
\begin{equation} \label{fftechoed_signal}
\hat{f}(k) = \prod_{i=1}^{N} \hat{f_i}(k)
\end{equation}
This step is foundamental in order to understand the concept of \textit{cepstrum},
for which the purpose is to separate the signal by convolution components.
So, the main idea is to take the logarithm of the product in the identity (\ref{fftechoed_signal}),
then taking the inverse fourier transform by linearity the single components $f_i$ will show up
as isolated peaks. Practically this task is done by the \textit{real cepstrum} and the \textit{power cepstrum}:
\begin{align} \label{cepstrum}
\mathcal{RC}(f) = \mathcal{F}^{-1}\{\ln(\abs*{\mathcal{F}\{f\}})\} & &
\mathcal{C}(f) = \abs*{\mathcal{F}^{-1}\{\ln(\abs*{\mathcal{F}\{f\}}^2)\}}^2 & &
\mathcal{C} = 4 \cdot \abs*{\mathcal{RC}}^2
\end{align}
\newpage
\section{\textit{Matlab} implementation}
\begin{homeworkProblem}
% When problems are long, it may be desirable to put a \newpage or a
% \clearpage before each homeworkProblem environment
\begin{homeworkSection}{(1) Straigh forward implementation}
Using the formulas given in the equations (\ref{inverse_dft}) and (\ref{dft}),
the resulting \textit{Matlab} code for the DFT is straight forward as it
shows to listing (\ref{mydft}).
\textit{Matlab} storage system of vectors and matrices is very flexible but
it's not clear how a matrix is treated when the initialisation is handled externally.
So, in order to determine the size of the output, the best approach is to
adapt it to the input size. That's what the lines (3) and (5) do: they take the input size format
and expand it into the output initialization.
\end{homeworkSection}
\matlabscript{mydft}{A \textit{matlab} implementation of the DFT}
\begin{homeworkSection}{(2) Testing and complexity analysis}
\begin{minipage}{\textwidth}
\hspace{-0.5cm}
\begin{minipage}{0.5\textwidth}
There are two essential aspects to check whether the algorithm works:
\begin{enumerate}
\item No bugs, all behaviours are expected
\item Lowest possible time complexity
\end{enumerate}
Both tasks can be executed using a script which compares the code execution of (\ref{mydft}) with
the \textit{Matlab} built-in \texttt{fft} (see file \texttt{test.m}).
At first glance, the complexity is supposed to be $\mathcal{O}(N^2)$ where $N$ is the size of the input array,
which is justified by the fact that the code loops $N$ times (index $j$) nested into another
$N$ sized loop (index $i$).
\end{minipage}
\hspace{1cm}
\begin{minipage}{0.55\textwidth}
\resizebox{\textwidth}{!}{
\input{graphs/testmydft.tex}
}
\end{minipage}
\vspace{1cm}
\end{minipage}
On the other hand the \textit{Matlab} built-in function \texttt{assert} helps in order to verify
that the output of \texttt{mydft} and \texttt{fft} correspond.
The graph in figure (\ref{complexgraph}) shows that the elapsed time increases
polynomially at the second order at the augmentation of the size of the input buffer.
\end{homeworkSection}
\end{homeworkProblem}
\clearpage
\begin{homeworkProblem}
\begin{homeworkSection}{(1.a) Straigh forward implementation of the inverse DFT}
Analogly to the previous point, the \textit{Matlab} implementation of the inverse
fourier transform is coded directly from the equation (\ref{inverse_dft}).
Knowing that the form of the algorithm is exactly the same, then the time complexity
is always $\mathcal{O}(N^2)$.
\matlabscript{mydftinverse}{A \textit{Matlab} implementation of the inverse DFT}
\end{homeworkSection}
\begin{homeworkSection}{(1.b) Inverse DFT test}
The best method to prove that the inverse DFT is working fine, is evaluating
it to a set previously tranformed values with the \texttt{mydft} function.
Then it's enought to take a random set using the \textit{Matlab} function \texttt{rand},
transform it with \texttt{mydft} and retransform it back with \texttt{mydftinverse}, as
shown in the script listing (\ref{testinverse}).
\matlabscript{testinverse}{Inverse DFT testing script}
\end{homeworkSection}
\begin{homeworkSection}{(1.c) Power cepstrum implementation}
Taking the previously given definition (equation (\ref{cepstrum})),
the following code is produced:
\matlabscript{mypowercepstrum}{Cepstrum implementation}
\end{homeworkSection}
\newpage
\begin{homeworkSection}{(1.d) Real cepstrum relation}
Following the same principle of the previous point,
using a random generated sequence it's easy to prove that
the built-in function \texttt{rceps} follows the third relation
given in the equation (\ref{cepstrum}).
\matlabscript{testpowercepstrum}{Cepstrum testing script}
\end{homeworkSection}
\begin{homeworkSection}{(2.a, 2.b and 2.c) Loading and playing audio file}
The built-in function \texttt{audioread} takes the file name as parameter and
returns the sample array \texttt{S} and the sampling frequency \texttt{sf} as
it's shown in the script below.
Then the \texttt{sound} function allows to send the samples to the audio output sink,
which consequently transforms them into a physical sound by the speakers.
\matlabscript{playsound}{Loading and playing an audio file in \textit{Matlab}}
Graphically visualizing the output of a loaded file,
the plots in figure (\ref{samp_plots}) show the output
of four audio files choosen in this regard: \texttt{audio1.wav}, \texttt{audio2.wav}, \texttt{audio3.wav}, \texttt{audio4.wav}.
As they're supposed to be, all the four signals are periodic and "sinelike",
in the sense that they are clearly composed by a finite series of sine and cosine.
\end{homeworkSection}
\begin{figure}[p]
\hspace{-2cm}
\begin{subfigure}{0.6\textwidth}
\resizebox{\textwidth}{!}{
\input{graphs/audio1/sample.tex}
}
\caption{Audio 1}
\label{samp_plots1}
\end{subfigure}
\hspace{0.5cm}
\begin{subfigure}{0.6\textwidth}
\resizebox{\textwidth}{!}{
\input{graphs/audio2/sample.tex}
}
\caption{Audio 2}
\label{samp_plots2}
\end{subfigure}
% TODO, this funcking spacing
\hspace{-2cm}
\begin{subfigure}{0.6\textwidth}
\resizebox{\textwidth}{!}{
\input{graphs/audio3/sample.tex}
}
\caption{Audio 3}
\label{samp_plots3}
\end{subfigure}
\hspace{0.5cm}
\begin{subfigure}{0.6\textwidth}
\resizebox{\textwidth}{!}{
\input{graphs/audio4/sample.tex}
}
\caption{Audio 4}
\label{samp_plots4}
\end{subfigure}
\caption{Audio samples plots}
\label{samp_plots}
\end{figure}
\begin{homeworkSection}{(2.d and 2.e and 2.f) Fourier transform of the signals}
Taking the fourier transform of all the signals the foundamental frequencies
are given by the distances between two consequent enough high peaks.
This definition isn't formal at all but essentially the other "micro peaks"
showing up in the figures (\ref{fourier_plots}) are negligible and they are interpreted as noise.
\begin{table}[h]
\centering
\begin{tabular}{|c|c|c|c|}
\hline
Audio 1 & Audio 2 & Audio 3 & Audio 4 \\
\hline
$522.8$ & $259.9$ & $253.1$ & $515.0$ \\
\hline
\end{tabular}
\caption{Foundamental frequencies in expressed in \si{\hertz}}
\label{foundfreq}
\end{table}
In the table (\ref{foundfreq}) the foundamental frequences are presented for
each audio. The plot referenced as "$\sin$ check" is the fourier transform
of a pure sine function holding the founded frequency for each audio.
Precisely the form is the following:
\begin{equation} \label{sincheck}
f_{check}(t) = 0.25 * \sin(2 \pi f_0 t)
\end{equation}
where $f_0$ is the foundamental frequency.
What's found by this analysis is that the audio 1 and 4 share the same
foundamental frequency and the same could be said comparing the audio 2 and 3.
Furthermore audio 1 and 4 are doubled in frequency than 2 and 3, which means
that the played musical note is the same for all files but the first pair has a difference
of an octave (12 half tones) with respect to the second one.
I recall that the transformation between musical notes and the frequencies is logarithmic and
adding 12 half tone is equivalent to double the frequency.
The real difference between audio 1 and audio 4 is on the echo, which is detected via
the \textit{Cepstral} analysis.
\end{homeworkSection}
\begin{figure}[p]
\hspace{-2cm}
\begin{subfigure}{0.6\textwidth}
\resizebox{\textwidth}{!}{
\input{graphs/audio1/fourier.tex}
}
\caption{Audio 1}
\label{fourier_plots1}
\end{subfigure}
\hspace{0.5cm}
\begin{subfigure}{0.6\textwidth}
\resizebox{\textwidth}{!}{
\input{graphs/audio2/fourier.tex}
}
\caption{Audio 2}
\label{fourier_plots2}
\end{subfigure}
% TODO, this funcking spacing
\hspace{-2cm}
\begin{subfigure}{0.6\textwidth}
\resizebox{\textwidth}{!}{
\input{graphs/audio3/fourier.tex}
}
\caption{Audio 3}
\label{fourier_plots3}
\end{subfigure}
\hspace{0.5cm}
\begin{subfigure}{0.6\textwidth}
\resizebox{\textwidth}{!}{
\input{graphs/audio4/fourier.tex}
}
\caption{Audio 4}
\label{fourier_plots4}
\end{subfigure}
\caption{Audio fourier transform plots}
\label{fourier_plots}
\end{figure}
\begin{homeworkSection}{(3) Cepstral analysis}
\noindent
\begin{minipage}{\textwidth}
\begin{minipage}{0.6\textwidth}
As mentioned in the section \textit{Cepstrum and echo} the \textit{cepstrum}
is the correct tool in order to determine whether echo is present or not in a signal.
Particularly the \textit{cepstral} analysis is an anagram shifting of \textit{spectral}
analysis, but instead of detecting the frequencies each pure sine
composing the signal could have a phase shift and that's what the \textit{cepstral}
analysis shows up. More concretely a redondant signal is always essentially composed by
the same frequency and the characteristic which is creating the echo is the superposition
of phase shifted sines.
The code in listing (\ref{mypowercepstrum}) corresponds to a \textit{Matlab} implementation
of the power cepstrum. The domain of the resulting function of a cepstrum is exactly the same
as the one of the original function, then its unit of measure is cleary the one of a time,
the \texttt{second} for the \textit{SI} unit system.
\end{minipage}
\hspace{1cm}
\begin{minipage}{0.4\textwidth}
\includegraphics[width=\textwidth,keepaspectratio]{redaction/cepss.png}
\captionof{figure}{Visualisation of echo in signal \cite{echocepss}}
\end{minipage}
\vspace{0.5cm}
\end{minipage}
Taking a look to the fourier transforms in figure (\ref{fourier_log_plots}),
in audio 1 the main peaks are much higher than the ones in audio 4.
This fact suggests that the audio 4 is more complex as signal, in sense
that holds a major echo.
Then analysing the cepstrum in figure (\ref{cepstral_graph}),
it's possible to see that the area generated integrating the signal is much more
distributed along all the domain in audio 4. Given that the curve
integrals evaluating using the \textit{trapezium rule} (file \texttt{discr\_integral.m}),
are supposed to be $1$ because of the $L^1$ normalization, then it's justified to say
that the major peaks of audio 1 are higher than the ones of audio 4.
This conclusion confirms the preceding intuition.
In another point of view, in the middle of the domain the audio 4
presents a huge quantity of micro peaks highlighted by the cepstrum plot, when the
the audio 1 is approximatively flat in that region. This is clear sign of temporal
disorder in the information collection of the audio 4 because each micro peak corresponds
to a small contribution of a shifted sine, then an echo.
\begin{figure}[p]
\hspace{-2cm}
\begin{subfigure}{0.6\textwidth}
\resizebox{\textwidth}{!}{
\input{graphs/audio1/fourierlog.tex}
}
\caption{Audio 1}
\label{fourier_log_plots1}
\end{subfigure}
\hspace{0.5cm}
\begin{subfigure}{0.6\textwidth}
\resizebox{\textwidth}{!}{
\input{graphs/audio4/fourierlog.tex}
}
\caption{Audio 4}
\label{fourier_log_plots4}
\end{subfigure}
\caption{Audio fourier transform plots}
\label{fourier_log_plots}
\end{figure}
\begin{figure}[p]
\centering
\resizebox{\textwidth}{!}{
\input{graphs/cepstral.tex} % TODO, adjust number format
}
\caption{Power cepstrum comparison evaluated audio 1 and 4}
\label{cepstral_graph}
\end{figure}
\end{homeworkSection}
\end{homeworkProblem}
\newpage
\section{Fast fourier transform (FFT) - Radix-2 algorithm}
\subsection{Derivation}
Let $\{x_n\}_{0 \le n < N}$ where $N = 2^r$ and $r \in \mathbb{N}$.
Taking the DFT definition (equation (\ref{dft})) let's notate the
output as $\hat{X}_k$, precisely in the following way:
\begin{equation} \label{fft_def}
\hat{X}_k(\{x_n\}) = \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N}nk}, \quad 0 \le k < N
\end{equation}
Then by separing even indexes from odd ones, the previous expression takes another form:
\begin{equation} \label{fft_even_odd}
\hat{X}_k(\{x_n\}) = \hat{X}_k(\{x_{2n}\}) + e^{-\frac{2\pi i}{N}k} \cdot \hat{X}_k(\{x_{2n + 1}\})
\end{equation}
where $\{x_{2n}\}$ and $\{x_{2n + 1}\}$ denote $\{x_n\}$ evaluated on even and respectively odd indices.
Now taking the second part of the output, it's straight forward to demostrate the following identity:
\begin{equation} \label{fft_even_odd_second}
\hat{X}_{k + \frac{N}{2}}(\{x_n\}) = \hat{X}_{k+\frac{N}{2}}(\{x_{2n}\}) - e^{-\frac{2\pi i}{N}k} \cdot \hat{X}_{k+\frac{N}{2}}(\{x_{2n + 1}\})
\end{equation}
Notice those sublists have size $\frac{N}{2}$, then, using the DFT periodicity (equation (\ref{dft_period})), the right hand expressions
are periodic in their input size which is exactly $\frac{N}{2}$. This statement can be used to reduce the equation (\ref{fft_even_odd_second})
to:
\begin{equation} \label{fft_even_odd_second}
\hat{X}_{k + \frac{N}{2}}(\{x_n\}) = \hat{X}_{k}(\{x_{2n}\}) - e^{-\frac{2\pi i}{N}k} \cdot \hat{X}_{k}(\{x_{2n + 1}\})
\end{equation}
In summary the following relations are deduced:
\begin{equation} \label{radixtwo_sum}
\forall 0 \le k < \frac{N}{2}:
\begin{cases}
& \hat{X}_k(\{x_n\}) = \textcolor{red}{\hat{X}_k(\{x_{2n}\})} + e^{-\frac{2\pi i}{N}k} \cdot \textcolor{blue}{\hat{X}_k(\{x_{2n + 1}\})} \\
& \hat{X}_{k + \frac{N}{2}}(\{x_n\}) = \textcolor{red}{\hat{X}_{k}(\{x_{2n}\})} - e^{-\frac{2\pi i}{N}k} \cdot \textcolor{blue}{\hat{X}_{k}(\{x_{2n + 1}\})}
\end{cases}
\end{equation}
\paragraph{Time complexity evaluation}
This final expression allows to compute the first and the second part of the output in only one step.
The computation depth of the recursion is supposed to be $r$ and the output assignment takes exactly the size
of the input in the recursion step, which is initially $N$. Then the time complexity is $\mathcal{O}(N \cdot r)$,
but by definition $r = \log_2(N)$, which means that the resulting complexity becomes $\mathcal{O}(N \cdot \log_2(N))$.
\newpage
\subsection{Implementation}
% third section
\begin{homeworkProblem}
\begin{homeworkSection}{Better than a straigh forward implementation}
\noindent
\begin{minipage}{\textwidth}
\begin{minipage}{0.5\textwidth}
Using the relationships found in the previous section (equation \ref{radixtwo_sum}), the radix-2 algorithm
written in \textit{Matlab} takes the form shown in script listing (\ref{myfft}).
The first process consist in a control in the input size: using the bitwise operators' logic
it's possible to determine whether $N$ is a pure power of 2 in $\log_2(N)$ passages.
This task is performed in the \texttt{input\_legit} function.
Once verified the conditions, it's possible to enter into the \texttt{fft\_step} function,
which computes by \textit{radix-2} the discrete fourier transform of the input.
At the end another check on the output size is done in order to verify whether it has to be
transposed or not.
The graph in figure (\ref{complexFFT}) shows that the time complexity of the fast fourier transform
computed with the \textit{radix-2} algorithm follows the same polynomial degree of the built-in
\texttt{fft}.
\end{minipage}
\hspace{1cm}
\begin{minipage}{0.55\textwidth}
\resizebox{\textwidth}{!}{
\input{graphs/testmyfft.tex}
}
\captionof{figure}{Comparison in time complexity with \textit{Matlab} built-in \texttt{fft}}
\label{complexFFT}
\end{minipage}
\vspace{0.5cm}
\end{minipage}
The reason of the offset difference between the two curves is that the time that an elementary
operation is computed is a bit lower for the \texttt{fft} function. This is due to the fact that \texttt{fft} is precompiled, then
the sub-implementation calls directly the machine language, when \texttt{myfft} is using the \textit{Matlab} interpreter.
\matlabscript{myfft}{\textit{Matlab} implementation of Radix-2 FFT}
\end{homeworkSection}
\end{homeworkProblem}
\section{Image filtering and dislocation index analysis}
\begin{homeworkProblem}
\begin{homeworkSection}{Setup (1): image loading and constrast}
In order to load an image in \textit{Matlab} the built-it functions \texttt{imread} and \texttt{rgb2gray}
are what's needed. The first one loads the image in the \texttt{rgb} (red, green, blue) format, which means that the output is
a bidimensional discretized map $N_i \times N_j \times 3$. Each point of the picture, in fact, refers to a 3 sized array
of which the values are bounded between $0$ and $255$.
Being the purpose of this exercice working on a simple environment, the direct simplification is to convert the image into
a gray scaled format, where each point is described by a simple intensity value always bounded between $0$ and $255$.
This conversion is applied directly by the \texttt{rgb2gray} function and the result will be a $N_i \times N_j$ grid/matrix.
The constrast is applied manually using a linear scaling transformation such that shifts all minimal values
to zero and scales by the maximal gap between the two extrema, as it's shown in the listing (\ref{contrast}).
Sometimes, when a linear scaling isn't enough, a logarithmic transformation is taken on data before being amplified
by the \texttt{constrast} function.
\matlabscript{contrast}{Linear constrast function}
\end{homeworkSection}
\begin{homeworkSection}{Setup (2): filtering and developped tools}
Applying a band pass filter on a fourier transformed image means that
some regions are left unchanged, when some others are obscured (set to the minimal value or zero).
In order to fulfill this purpose in the most simple way, some important classes have been developped:
\begin{itemize}
\item \texttt{FDomain}: a general class handling a geometrical 2D domain
\item \texttt{FDiscMap}: a wrapper for the image handling transformation between the discretized grid and a continous space
\item \texttt{FEllipse}: a specialized \texttt{FDomain}, it describes an ellipse
\item \texttt{FArcEllipse}: a specialized \texttt{FEllipse}, it describes an angular portion of an ellipse
\item \texttt{FFilter}: a class which handles the filtering task, collecting domains and obscuring where it's specified
\item \texttt{FBandPassFilter}: a specialized \texttt{FFilter}, it obscures as it's shown in the figure (\ref{filtered})
\end{itemize}
\end{homeworkSection}
\begin{homeworkSection}{Resulting fourier space}
\noindent
\begin{figure}[h!]
\centering
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth,keepaspectratio]{redaction/grayscaled}
\caption{Grey-scaled image file \texttt{stm.png}}
\label{grayscaledimg}
\end{subfigure}
\hspace{0.5cm}
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth,keepaspectratio]{redaction/fouriertransf}
\caption{Fourier transformed image file \texttt{stm.png}}
\label{transformedimg}
\end{subfigure}
\end{figure}
The procedure used to obtain the fourier transformed image showed in figure (\ref{transformedimg})
is using the \texttt{fft2} function followed by a \texttt{fftshift} centering.
\texttt{fft2} is \textit{Matlab} built-in function and it performs a bidimensional discrete fourier transform
as the continous transform can be generalised on an integral in multiple dimensions.
In this case the peaks of the original periodic pattern are visualized as an intensity singularity (maximum white),
then just like the one dimensional expression, the properties and the behaviour remains unchanged.
This means also that the risulting matrix of a bidimensional FFT must be swapped and that's what the \texttt{fftshift}
performs in two dimensions.
\noindent
\begin{minipage}{\textwidth}
\begin{minipage}{0.6\textwidth}
\paragraph{About the fourier space pattern}
As mentioned at the course of \textit{Solid State Physics I and II} \cite{solidstate},
the fourier space of a \textit{Bravais} lattice is called the \textit{reciprocal} lattice.
In the particular case of a graphene lattice, as also the study of Dodd Gray \cite{doddgray} mentions,
its structure is hexagonal containing one atom per primitive cell.
This implies that computing the reciprocal lattice parametrisation, the resulting geometry is characterized by an hexagonal distribution
of the bragg peaks and the first \textit{Brillouin} zone has the same form (as shown in figure (\ref{lattice})).
The image in figure (\ref{grayscaledimg}) shows clearly a graphene-like pattern, then the reciprocal space has inherits the properties
described above, which is what's shown in figure ((\ref{transformedimg})).
\end{minipage}
\hspace{1cm}
\begin{minipage}{0.4\textwidth}
\includegraphics[width=\textwidth,keepaspectratio]{redaction/lattice.png}
\captionof{figure}{Reciprocal lattice of the graphene \cite{doddgray}}
\label{lattice}
\end{minipage}
\vspace{0.5cm}
\end{minipage}
\end{homeworkSection}
\begin{homeworkSection}{Filtering and measuring index dislocation}
\noindent
\begin{figure}[h!]
\centering
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth,keepaspectratio]{redaction/filter}
\caption{Filtered fourier space}
\label{filteredfourierimg}
\end{subfigure}
\hspace{1cm}
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth,keepaspectratio]{redaction/dislocation}
\caption{Filtered image}
\label{filteredimg}
\end{subfigure}
\caption{}
\end{figure}
\noindent
\begin{minipage}{\textwidth}
\begin{minipage}{0.6\textwidth}
In order to indentify the dislocation it's necessary to obscurate all the fourier space
except for a specific bragg peak.
As the figure (\ref{filteredfourierimg}) shows, the best peak to chose is the one located in
$(r = 0.077, \phi = \pi/2)$ in a polar coordinate system, where $r = 1$ reaches the top of the image.
The figure (\ref{filteredimg}) shows that, applying the correct filter, it's
possible to clearly see the dislocation.
Taking a look to the picture (\ref{filtered_zoomed}), between the green
and the blue spots there is a difference of 3 in right side, when in the left
side is 7. Then the measured dislocation is $N_{disl} = N_L - N_R = 4$.
\end{minipage}
\hspace{1cm}
\begin{minipage}{0.4\textwidth}
\includegraphics[width=\textwidth,keepaspectratio]{redaction/dislocation_zoomed}
\captionof{figure}{Dislocation region, zoomed}
\label{filtered_zoomed}
\end{minipage}
\vspace{0.5cm}
\end{minipage}
\end{homeworkSection}
\end{homeworkProblem}
\newpage
\section{Conclusion}
What's highlighted in this experience is that the fourier transform is one of the most powerful tools
for sound and image analysis and correction.
In general for every period signal, in any dimension, this tool allows to
separe onto the composing frequencies, detecting echo and dislocation.
It's handful too: even if not intuitive, it's possible to reduce the algorithm time
complexity in order to compute it faster. The case of the \textit{radix-2} is just
a limited example, because not any input is allowed, but there exist many other variants
of the fast fourier transform.
In conclusion the numerical application of this tool has given optimal results and
no particular instability.
\end{spacing}
\def\refname{D\MakeLowercase{ocumentation and sources}}
\begin{thebibliography}{99}
\bibitem{solidstate}
\url{https://edu.epfl.ch/coursebook/en/solid-state-physics-i-PHYS-309}
\bibitem{doddgray}
\url{http://community.wvu.edu/~miholcomb/graphene.pdf}
\end{thebibliography}
\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Event Timeline

c4science · Help