\documentclass{article}
\usepackage{amsmath,cite}
\usepackage{graphicx}

\usepackage{tikz}
\usepackage{tikz-qtree}

\usepackage{amssymb}

\usepackage{url}

\usepackage[framed,numbered,autolinebreaks,useliterate]{mcode}

\graphicspath{{figures/}}

\title{ScatNet : \\ a MATLAB Toolbox for Scattering Networks}

\author{Laurent Sifre, Joakim And\'{e}n, Michel Kapoko, \\ Edouard Oyallon, and Vincent Lostanlen}

\date{v0.2, october 2013}

\begin{document}
    
\maketitle

\section{Introduction}

\subsection{Overview}
ScatNet is a MATLAB scientific toolbox which provides tools for the analysis and classification of digital signals, such as sounds, images or time series. Mostly owing to its properties of group invariance and stability to deformations, it has shown to achieve state-of-the art results in the challenges of music genre recognition, image, texture classification, and fetal heart rate characterization.
Its core feature relies on the construction of a \emph{scattering network}, i.e. a stack of signal processing layers of increasing width. Each layer consists in the association of a linear filter bank \mcode{Wop} with a non-linear operator, namely the complex modulus.
The \emph{scattering transform} of an input signal \mcode{x} is defined as the set of all paths that \mcode{x} might take from layer to layer. In this sense, the architecture of a scattering network closely resembles a convolutional deep network.

\subsection{Getting started}
\subsubsection{Audio signal processing}
An online tutorial for the scattering representation of audio signals is available at the following page : \\
\hfill \url{http://www.di.ens.fr/data/software/scatnet/quickstart-audio/} \hfill

\subsubsection{Image processing}
An online tutorial for the scattering representation of image signals is available at the following page : \\
\hfill \url{http://www.di.ens.fr/data/software/scatnet/quickstart-image/} \hfill

\subsection{Outline}
In this document, the objects in ScatNet are first reviewed from the general to the specific : Sections \ref{sec:scattering-transform}, \ref{sec:wavelet-factories} and \ref{sec:filters} respectively cover the architectures of the scattering operator, of each layer in the network, and of each filter in a layer.
Afterwards, Section \ref{sec:utilities} presents some utilities that facilitate the manipulation and formatting of scattering transforms, for the purposes of visualization and discriminative learning.
Section \ref{sec:display} is dedicated to the display tools of ScatNet, both for one-dimensional and two-dimensional input signals.
Finally, Section \ref{sec:classification} thoroughly describes the classification pipeline of affine space and support vector classifiers for scattering representations.

\section{Scattering transform \label{sec:scattering-transform}}

\subsection{Prototype \label{sub:scattering-transform-prototype}}
Assuming that the cell array of linear operators \mcode{Wop} is already built, the computation of the scattering transform of some input signal \mcode{x} is a fast operation, supported by the function \mcode{scat}.
Interestingly, the array \mcode{x} may either be of dimensions:

\begin{enumerate}
\item $(N \times 1)$, for a time series or an audio file of length $N$,
\item $(N_1 \times N_2)$ , for an image of height $N_1$ and width $N_2$, or
\item $(N_1 \times 1 \times N_3)$, for a set of $N_3$ audio files, each of length $N_1$.
\end{enumerate}

In all cases, the network of scattering coefficients is obtained through the single command :

\begin{lstlisting}
S = scat(x, Wop);
\end{lstlisting}

\mcode{S} is a cell array whose length \mcode{M+1} is the same as \mcode{Wop}. The integer \mcode{M} is the maximal order of the scattering transform, and is noted $\overline{m}$ in academic papers. From a theoretical point of view, the order $m$ ranges from $0$ to $\overline{m}$ ; but since MATLAB is one-based, the corresponding layer index must range from \mcode{1} to \mcode{M}, hence the need of an offset of \mcode{1} for the variable \mcode{m}.

Each layer $S_m$ of scattering coefficients proceeds from an averaging of modulus coefficients $U_{m-1}$ computed at the previous layer. These modulus coefficients can be separately obtained with \mcode{scat}, as an optional second output argument :

\begin{lstlisting}
[S, U] = scat(x, Wop);
\end{lstlisting}

\subsection{Implementation of the scattering transform}
Each operator \mcode{Wop\{1+m\}} performs two actions, leading to separate outputs :
\begin{enumerate}
\item an energy \emph{averaging} according to the largest scale, by means of a low-pass filter $\phi$, and
\item an energy \emph{scattering} along all scales, by means of band-pass filters $\psi_j$.
\end{enumerate}

After initializing \mcode{U\{1+0\}} to \mcode{x}, \mcode{scat} executes the following loop :

\begin{lstlisting}
for m = 0:M-1
   [S{1+m}, V] = Wop{1+m}(U{1+m});
   U{1+(m+1)} = modulus_layer(V);
end
S{1+M} = Wop{1+M)}(U{1+(M-1)});
\end{lstlisting}

The function \mcode{modulus_layer} computes the modulus of all signals in \mcode{V} while preserving its layer structure.

At the last layer, it is not necessary to perform any scattering, since the coefficients will not be further processed. Therefore, \mcode{Wop\{1+M\}} is merely an averaging operator, with a single output argument.

\subsection{Format of a scattering layer}

The layer format, in which every \mcode{S\{1+m\}} and \mcode{U\{1+m\}} are encoded, consists of two fields, namely \mcode{signal} and \mcode{meta}. The former is a cell array of real-valued signals, while the latter is a structure containing meta-information. In the sequel, \mcode{S\{1+m\}.meta} is an array, whose number of columns is equal to the number of signals within the layer.

Any user-defined function with one input and two outputs in the previous format may become an element of the cell array \mcode{Wop}, through an instruction of this kind :

\begin{lstlisting}
Wop{1+m} = @mth_user_defined_operator;
\end{lstlisting}

Most often, however, it is not necessary to do so ; indeed, ScatNet provides efficient linear operators by default, especially adapted to the processing of images and sounds. Therefore, the whole cell array \mcode{Wop} can be defined in one line through built-in "factories", which rely on the theory of wavelet analysis. These factories are reviewed in the next section.


\section{Wavelet factories \label{sec:wavelet-factories}}
All five built-in factories are available in the \mcode{core} directory. After reviewing their respective implementations and prototypes, we document the customization of the scattering transform through the \mcode{scat_opt} structure.

\subsection{Scattering representations of audio signals and time series \label{sec:wavelet-factory-1d}}
The function \mcode{wavelet_1d} gives an efficient representation of one-dimensional signals, such as sounds or time series. It bears the following prototype :

\begin{lstlisting}
[Wop,filters] = wavelet_factory_1d(size(x),filt_opt,scat_opt);
\end{lstlisting}

The input arguments \mcode{filt_opt} and \mcode{scat_opt} respectively gather the parameters of the wavelet filter banks and the scattering network ; note that both can be omitted at first glance. Section \ref{sub:filters-options} explains how to create custom filter banks. As regards \mcode{scat_opt}, the parameter most subject to change is the integer \mcode{scat_opt.M}, which represents the maximal order of the scattering transform, and equals $2$ by default. The \mcode{scat_opt} structure is further specified in Section \ref{sub:scattering-options}.

The first output argument is a cell array of function handles, which fit the required format in the prototype of \mcode{scat} (see Section \ref{sub:scattering-transform-prototype}). The second output is not mandatory for the processing of \mcode{x}, as it is only useful for the visualisation of the filter bank themselves: its architecture is reviewed in Section \ref{sub:filters-architecture}.

\subsection{Translation-invariant representations of images \label{sub:wavelet-factory-2d}}

The cell array of linear operators produced by the \mcode{wavelet_factory_2d} built-in represents filter banks of two-dimensional wavelets, providing a scattering representation of images which is invariant up to a planar translation. The prototype of \mcode{wavelet_factory_2d} is exactly the same as in the one-dimensional case (see Section \ref{sec:wavelet-factory-1d}), and all previous remarks hold: 

\begin{lstlisting}
[Wop,filters] = wavelet_factory_2d(size(x),filt_opt,scat_opt);
\end{lstlisting}

In addition, ScatNet provides another wavelet factory for images, called \mcode{wavelet_factory_2d_pyramid}, and relying on the cascade algorithm.  While pursuing the same goal as \mcode{wavelet_factory_2d}, it happens to be computationally faster, and doesn't require to be given the size of the signal \mcode{x} as a first input; hence the following prototype:

\begin{lstlisting}
[Wop,filters] = wavelet_factory_2d_pyramid(size(x),filt_opt,scat_opt);
\end{lstlisting}

In spite of these assets, \mcode{wavelet_factory_2d_pyramid} is not as customizable as its counterpart, especially in terms of anti-aliasing and numerical approximations. As a matter of fact, these two functions mainly differ in the format of \mcode{filters} : while \mcode{wavelet_factory_2d}, like \mcode{wavelet_1d}, stores its filters in the Fourier domain, \mcode{wavelet_factory_2d_pyramid} uses their spatial support.

\subsection{Roto-translation invariant representations of images}
Finally, the function \mcode{wavelet_factory_3d} enables a scattering representation of images which is invariant to both translations and rotations \cite{sifre13}. This original feature dramatically improves the performance of texture recognition in comparison with the translation-invariant representation provided by \\
\mcode{wavelet_factory_2d}. 

\begin{lstlisting}
[Wop,filters,filters_rot] = wavelet_factory_3d(size(x), filt_opt, filt_rot_opt,scat_opt);
\end{lstlisting}

With \mcode{wavelet_factory_3d}, at the first order, the energy is scattered along the spatial dimensions of the image, similarly to \mcode{wavelet_factory_2d}. At higher orders, however, \mcode{wavelet_factory_3d} computes a one-dimensional wavelet transform along the orientations of the previous coefficients. Therefore, for a scattering transform of order $2$, each coefficient bank bears three indices — vertical scale, horizontal scale, and angle — hence the "3d" denomination.

Whereas the structure \mcode{filt_opt} contains the parameters of the filter bank at the first order, \mcode{filt_rot_opt} is related to the roto-translation wavelets of higher orders.

Likewise the translation-invariant case, ScatNet provides an alternative built-in factory for \mcode{wavelet_3d}, which relies on the cascade algorithm (see Section \ref{sub:wavelet-factory-2d}). All previous remarks hold.

\begin{lstlisting}
[Wop,filters,filters_rot] = wavelet_factory_3d_pyramid(filt_opt,filt_rot_opt,scat_opt);
\end{lstlisting}

Contrary to what their name might suggest, neither \mcode{wavelet_3d} nor \\
\mcode{wavelet_3d_pyramid} are adapted to an input array \mcode{x} of size $(N_1 \times N_2 \times N_3)$. Up to now, there is no built-in factory for the processing of three-dimensional signals, such as video clips or voxel-based measurements. Nevertheless, it is not difficult to derive a custom pipeline that matches the geometric dimension of the data, in a similar fashion to the available factories.

\subsection{Customized wavelet factory \label{sub:scattering-options}}
All built-in wavelet factories presented above may take a structure \mcode{scat_opt} as last argument, offering a high degree of flexibility for the user. In this section, we review the available fields in \mcode{scat_opt} according to the chosen factory. For instance, one may activate the single floating-point precision by the following instruction :

\begin{lstlisting}
scat_opt.precision = 'single';
\end{lstlisting}

If \mcode{scat_opt} is an empty data structure or is omitted, all these fields are set to their default values. If the user provides an additional field to \mcode{scat_opt}, or spell one of the available fields mistakenly, ScatNet throws an error.

\subsubsection{Maximal order of the scattering transform}
The positive integer \mcode{scat_opt.M} encodes the maximal order of the scattering transform (see Section \ref{sub:scattering-transform-prototype}), i.e. the depth of the associated scattering network. When set to \mcode{1}, the scattering transform is merely the modulus of a wavelet transform. By default, it is equal to \mcode{2},  whatever be the chosen built-in factory. In most cases, increasing \mcode{scat_opt.M} marginally improves classification results, yet at a great computational cost.

\subsubsection{Oversampling of the scattering coefficients}
The positive integer \mcode{scat_opt.oversampling} encodes the $log_2$ of the desired oversampling factor. By default, \mcode{scat_opt.oversampling} is equal to \mcode{1}, which means that the scattering coefficients are computed right at the critical sample rate, and that no oversampling is carried out at all. Although oversampling $S_x$ may be useful for post-processing or visualisation, note that the required memory increases exponentially with the integer \mcode{scat_opt.oversampling}.

\subsubsection{Resolution of the input signal}
The integer \mcode{scat_opt.x_resolution} encodes the $log_2$ of the resolution of the input signal \mcode{x}, with respect to the finest resolution of the filter bank. By default, \mcode{scat_opt.x_resolution} is equal to \mcode{0}, which means that increasing \mcode{scat_opt.x_resolution} is equivalent to an interpolation of \mcode{x}, whereas decreasing it below $0$ is equivalent to a downsampling operation. Likewise the parameter \mcode{scat_opt.oversampling} (see above), the required memory increases exponentially with the integer \mcode{scat_opt.oversampling}. In summary, relatively to the scattering transform, the integer \mcode{scat_opt.x_resolution} behaves as a prior interpolation, while \mcode{scat_opt.oversampling} behaves as a posterior interpolation.

\subsubsection{Floating-point precision}
The string \mcode{precision} encodes the numerical precision of the scattering coefficients. Its default value is \mcode{'double'}, which means that real numbers occupy 8 bytes (64 bits) of computer memory, in accordance with the default format in MATLAB. In order to alleviate computational load, the floating-point precision may be manually downgraded to 4 bytes by setting \mcode{scat_opt.precision } to \mcode{'single'}. However, note that this choice may damage the discriminative performance of the resulting machine learning algorithms.

\section{Filter banks \label{sec:filters}}
Filter banks correspond to predefined sets of filters. The nature of the used filters is highly important since it determines the result of the scattering transform. Filters bank can be created using the \mcode{wavelet_factory} set of functions, and have a similar skeleton in 1D, 2D or 3D.

\subsection{Architecture \label{sub:filters-architecture}}
\begin{tikzpicture} [scale=.6]
\centering
\tikzstyle{level 1}=[sibling distance=80mm] 
\tikzstyle{level 2}=[sibling distance=40mm] 
\tikzstyle{level 3}=[sibling distance=15mm] 

 \node{filters(struct)}  
   child{ node{phi.}
    child{node{meta.}
      child{node{J}}}
    child{node{filter.}
       child{node{coefft\{\}}}
       child{node{type}}
       child{node{N}}}
       }
   child{node{psi.}
      child{node{filter\{\}}
        child{node{coefft\{\}}}
        child{node{type}}
        child{node{N}}
        }
       
       child{node{meta.}
            child{node{j}}
            child{node{theta}}
            }}
            
            
    child{node{meta.}
       child{node{J(:)}}
       child{node{Q(:)}}
       child{node{...}}
       };
 \end{tikzpicture}
 
 In each of this case, \mcode{filters} have the same architecture. The first field corresponds to a low-pass filter, the second to the band-pass filters and a field \mcode{meta} contains the information for the whole bank-filters.
 The band-pass filters are stored as cells of filters.


\subsection{Options \label{sub:filters-options}}


As described in the previous section, filter parameters are specified in the \mcode{filt_opt} structure passed to \mcode{wavelet_factory_1d}, \mcode{wavelet_factory_2d} or \mcode{wavelet_factory_3d}. Depending on the type of filter used, this structure can contain different parameters, but the following can always be specified:
\begin{itemize}
	\item \mcode{filt_opt.filter_type}: The wavelet type, such as `morlet\_1d', `morlet\_2d', `spline\_1d', for example (default \mcode{'morlet_1d'} for 1D signals, \mcode{'morlet_2d'} for 2D signals).
	\item \mcode{filt_opt.precision}: The numeric precision of the filters. Either \mcode{'double'} or \mcode{'single'} (default \mcode{'double'}).
\end{itemize}
Since the \mcode{filter_type} field determines the type of filter bank to create, it determines which other options can be specified. The rest of this section covers the different types of filters and their parameters.
\subsection{1D signals}
For 1D signals, it is often useful to specify different filter parameters for different orders. This is done by using arrays instead of scalar values for parameters in \mcode{filt_opt}. In the case of numeric parameters, this means that they are regular arrays, whereas for string parameters, such as \mcode{filter_type}, they are cell arrays.

For example, setting \mcode{filt_opt.filter_type = \{'gabor\_1d','spline_1d'\}}, means that the first-order filter bank should be a Gabor filter bank whereas the second should be a spline wavelet filter bank. If more filter banks are needed than parameters are available in \mcode{filt_opt}, the last element in each array is repeated as necessary.

The \mcode{meta.resolution} field has one row, indicating the extent the signal has been subsampled with respect to the original length.

In the case of wavelet transforms used as linear operators, an important field is \mcode{meta.j}, which has a variable number of rows, describing the scales of the wavelets used to compute the coefficient. For a first-order coefficient $|x\star\psi_{j_1}|\star\phi(t)$, this will simply be one row containing $j_1$. For a second-order coefficient $||x\star\psi_{j_1}|\star\psi_{j_2}|\star\phi(t)$, the first row will contain $j_1$ while the second will contain $j_2$.

Let us consider the example output \mcode{S} from the \mcode{scat} function using wavelet transforms as linear operators. In this case, \mcode{S\{m+1\}} will contain the $m$th-order scattering coefficients, with \mcode{S\{m+1\}.signal\{p\}} being the $p$th signal among these. Its scales $(j_1,j_2,\ldots,j_m)$ are specified in \mcode{S\{m+1\}.meta.j(:,p)}.

\subsubsection{1D Morlet/Gabor Filters}
The Gabor wavelets consist of Gaussian envelopes modulated by complex exponentials to cover the entire frequency spectrum. Morlet wavelets are derived from these by subtracting the envelope multiplied by a constant such that the integral of the filter equals zero. As they are closely related, the Gabor and Morlet wavelet have identical sets of options.

These are obtained by setting \mcode{filt_opt.filter_type} to \mcode{'morlet_1d'} or \mcode{'gabor_1d'}, respectively. The Morlet filters are the default filters for 1D signals. 

The Morlet/Gabor filter bank has the following options:
\begin{itemize}
	\item \mcode{filt_opt.Q}: The number of wavelets per octave. By default $1$.
	\item \mcode{filt_opt.J}: The number of wavelet scales.
	\item \mcode{filt_opt.B}: The reciprocal octave bandwidth of the wavelets. By default \mcode{filt_opt.Q}.
	\item \mcode{filt_opt.sigma_psi}: The standard deviation of the mother wavelet in space. By default calculated from \mcode{filt_opt.B}.
	\item \mcode{filt_opt.sigma_phi}: The standard deviation of the scaling function in space. By default calculated from \mcode{filt_opt.B}.
\end{itemize}
The maximal wavelet bandwidth (in space) is determined by $2^{J/Q}$ times the bandwidth of the mother wavelet, which is proportional to \mcode{sigma\_psi}. If \mcode{sigma\_psi} is smaller than a certain threshold, a number of constant-bandwidth filters are added, linearly spaced, to cover the low frequencies.

Again, we can specify different filter banks by setting \mcode{filt_opt.Q} and \mcode{filt_opt.J}, etc. to arrays instead of scalars. This is often useful if the nature of the signal is different at different orders, which is usually the case in audio. For this reason, the \mcode{default_filter_options} for \mcode{'audio'} specify \mcode{filt_opt.Q = [8 1]}, whereas for \mcode{'dyadic'}, \mcode{filt_opt.Q} is set to \mcode{1}. A higher frequency resolution is needed in audio for the first-order filter banks due to the highly oscillatory structure of the signals.

To calculate the appropriate \mcode{J} necessary for a given window size \mcode{T} and filter parameters \mcode{filt_opt}, there is the conversion function \mcode{T_to_J}, which is called using:
\begin{lstlisting}
	filt_opt.J = T_to_J(T, filt_opt);
\end{lstlisting}
This will ensure that the maximum wavelet scale, and thus the averaging scale of the lowpass filter, will be approximately \mcode{T}. Note that the previous prototype, \mcode{T_to_J(T,Q,B,phi_bw_multiplier)} is no longer supported in this version.

\subsubsection{1D Spline Filters}
The spline wavelet filter bank is a unitary filter bank as defined in \cite{mallatbook}. As a result, the scattering transform defined using these filters has perfect energy conservation \cite{stephane}.

In addition to the parameters listed above, the spline filter bank has the following options:
\begin{itemize}
	\item \mcode{filt_opt.J}: The number of wavelet scales.
	\item \mcode{filt_opt.spline_order}: The order of the splines. Only linear (spline order $1$) and cubic (spline order $3$) are supported.
\end{itemize}
The maximal bandwidth is specified here by $2^J$.

\subsection{1D Selesnick Filters}
Selesnick wavelet filters are compactly supported complex wavelets, whose real and imaginary parts have four vanishing moments and are nearly Hilbert transform pairs.

In addition to the parameters listed above, the Selesnick filter bank has the following options:
\begin{itemize}
  \item \mcode{filt_opt.J}: The number of wavelet scales.
\end{itemize}
The maximal bandwidth is specified here by $2^J$.

The code for the Selesnick wavelet filters was written by Shihua Cai, Keyong Li, and Ivan Selesnick. It has been included in ScatNet with their permission. The standalone package, which includes more functions and additional documentation, can be found at:

\texttt{http://eeweb.poly.edu/iselesni/WaveletSoftware/}.



\subsection{2D filters}
\label{sec_filters_2d}
\subsubsection{Morlet filters bank}
The Gabor wavelets consist of Gaussian envelopes modulated by complex exponentials to cover the entire frequency spectrum. Morlet wavelets are derived from these by subtracting the envelope multiplied by a constant such that the integral of the filter equals zero. 

A 2D Morlet filter bank consists in a Gaussian window 
\begin{equation*}
\phi_J (u) =  2^{-2j/Q} \phi(2^{-j/Q} (u,v))
\end{equation*}
and in dilated and rotated Morlet filters.
\begin{equation*}
\psi_{j, \theta}(u) = 2^{-2j/Q} \psi(r_\theta 2^{-j/Q} (u,v)) 
\end{equation*}
with 
\begin{eqnarray*}
j & \in & [0,J-1] \\
\pi^{-1} L^{-1} \theta & \in & [0,L-1]
\end{eqnarray*}
The full filter bank can be obtained by calling 
\mcode{morlet_filter_bank_2d}. 
\begin{lstlisting}
size_in = [128, 128];
filt_opt.L = 6;
filt_opt.Q = 1;
filt_opt.J = 3;
filters = morlet_filter_bank_2d(size_in, filt_opt);
\end{lstlisting}
This example builds a filters bank adapted to \mcode{128x128} images with the following options :
\begin{itemize}
	\item \mcode{filt_opt.Q} the number of wavelets per octave. By default $1$.
	\item \mcode{filt_opt.J}  the number of wavelet scales. By default $4$.
	\item \mcode{filt_opt.L}  the number of wavelet orientations. By default $8$.
\end{itemize}	



The mother Morlet filter is defined by
\begin{equation*}
\label{eq_mother_morlet}
\psi(u,v) = e^{ \frac{u^2 + s^2 v^2}{2 \sigma^2}}(K - e^{i u \xi})
\end{equation*}
whose name of the parameters are
\begin{itemize}
\item $\sigma$ the spread of the gaussian envelope 
\item $\xi$ the frequency of the oscillatory exponential
\item $s$ the eccentricity of the elliptical gaussian enveloppe
\end{itemize}
Setting \mcode{Q} and \mcode{L} influences the default parameters of the mother Morlet filter $\sigma$, 
$\xi$ and $s$. For example, increasing $L$ also increases the angular resolution by decreasing $\xi$.
This adaptive behavior can be overridden by manually setting $\sigma$, $\xi$, $s$ in the \mcode{filt_opt} structure.
\begin{itemize}
	\item \mcode{filt_opt.sigma_phi} the spread of $\phi$
	\item \mcode{filt_opt.sigma_psi} the spread $\sigma$ in (\ref{eq_mother_morlet})
	\item \mcode{filt_opt.xi_psi} the frequency $\xi$ in (\ref{eq_mother_morlet})
	\item \mcode{filt_opt.slant_psi} the eccentricity $s$ in (\ref{eq_mother_morlet})
\end{itemize}	
\subsubsection{Shanon filters bank}
Shanon wavelets are separable wavelets with a a finite support in frequency. Thus they suffer from a bad localization in time since their Fourier transform is not regular. Here they are designed to show filter banks that respects the Littlewood paley condition. Shanon mother wavelet of the Shanon filters bank is defined by
\begin{equation}
\label{eq_shanon_filters_bank}
\hat\psi(u,v) = e^{2i\pi(u+v)}1_{1>|u|>\frac{1}{2},1>|v|>\frac{1}{2}}
\end{equation}

In our implementation, these wavelets are set so that the Littlewood paley condition is respected.

The parameters of this filters bank are given by:

\begin{itemize}
    \item \mcode{filt_opt.J} the scale of the filters
\end{itemize}	

However, they are not designed for classification purposes, yet only for testing.

\section{Post-processing operations on scattering networks \label{sec:utilities}}
The signal processing tools originally provided in ScatNet were designed to address computer vision and machine listening problems, and can be evaluated through standard challenges in database classification. Most of machine learning tools require vector inputs of the same size with specifical renormalization and averaging(or pooling). In this section, we describe a full set of functions to manipulate the output of the scattering, in order to satisfy any conditions needed by a classifier.

\subsection{Renormalisation of high-order coefficients}
As described in \cite{joakim}, second-order scattering coefficients can be renormalized by dividing each $Sx(t,j_1,j_2)$ by $Sx(t,j_1)$, which decorrelates the amplitudes at the second layer from the first layer. The renormalisation utility in ScatNet has the following prototype :

\begin{lstlisting}
S_renorm = renorm_scat(S, epsilon);
\end{lstlisting}

where \mcode{S} is the scattering transform to be renormalized, and \mcode{epsilon} is a regularization constant that prevents the fraction from diverging when the first order becomes zero.

\subsection{Renormalization of first order coefficients for one-dimensional scattering representation}
Sometimes, in audio, renormalizing first order scattering coefficient of the signal by the zeroth order coefficients of the scattering representation of its absolute value (the zeroth order coefficients of the scattering representation of x are almost nul). The intuition is that this way, our first order coeffiients will be invariant to amplitude modification. 
In ScatNet, this can be achieved with the function \mcode{rdivide_scat} whose prototype is 

\begin{lstlisting}
  X_ren =  rdivide_scat(X,Xn,layerNb,epsilon);
\end{lstlisting}

where:

\begin{itemize}

    \item \mcode{X} is the scattering transform S or U whose layer should be normalized.
	\item \mcode{Xn} is the scattering transform S or U whose layer would be used to normalize the one of X
	\item \mcode{layerNb} is the indice of the layer of Xn to use for normalization
    \item \mcode{epsilon}  is a regularization constant that prevents the fraction from diverging when the first order becomes zero. It is set to $2^(-20)$ by default.

\end{itemize}

You would then need to proceed similarly as indicated in the sample code below:

\begin{lstlisting}
  filt1_opt.filter_type = {'gabor_1d','morlet_1d'};
  filt1_opt.Q = [8 1];
  filt1_opt.J = T_to_J(T,filt1_opt);

% set the scattering options
sc1_opt = struct();
sc1_opt.M = 1;

Wop1 = wavelet_factory_1d(N1,filt1_opt,sc1_opt);

% function to compute the scattering transform of  the signal abs(x)
absx_scat = @(x)(scat(abs(x),Wop1));

filt2_opt = filt1_opt;
sc2_opt.M = 2;


Wop2 = wavelet_factory_1d(N1, filt2_opt, sc2_opt);

% function to compute the scattering transform of the signal x

x_scat = @(x)(scat(x,Wop2));

% Compute the order one renormalized scattering transform of x
Sx = rdivide_scat(x_scat(x),absx_scat(x),1);

\end{lstlisting}


\subsection{Full Renormalization of one-dimensional scattering representation}
When working with the second order scattering representation of an audio signal, the combination of the functions \mcode{renorm_scat} and \mcode{rdivide_scat} enables a fully renormalized second-order scattering representation. This makes the feature insensible to amplitude variation, resulting to a considerable improvement in some cases. Keeping the framework introduced in \mcode{rdivide_scat}, one may compute this feature using the lines below:


\begin{lstlisting}
 absx_scat = @(x)(scat(abs(x),Wop1));
 x_scat = @(x)(scat(x,Wop2));
 Sx = rdivide_scat(renorm_scat(x_scat(x)),absx_scat,1);

\end{lstlisting}


\subsection{Logarithmic transformation of a scattering network}
Some classification frameworks, such as the MFC in speech recognition, rely on log-power representation of the frequency distribution. ScatNet has the ability to reproduce this behaviour by computing the logarithm of scattering coefficients. Log-scattering coefficients are also useful for visualisation, as they shrink the range of values to a small interval of orders of magnitude.

\begin{lstlisting}
log_S = log_scat(S,epsilon);
\end{lstlisting}
The positive real number \mcode{epsilon} may be omitted : by default, it is set equal to $2^-(20)$, i.e. $10^-6$. It is a small constant taht is added to each coefficient, in order to keep finite results when the original coefficient is null.

\begin{lstlisting}
log_layer = log_scat(S{1+m},epsilon);
\end{lstlisting}
For several signals, it is sometimes necessary to compute the logarithm of scattering coefficients. This is achieved in ScatNet using the \mcode{log_scat} function. Its signature is

\begin{lstlisting}
S_log = log_scat(S, epsilon);
\end{lstlisting}
where \mcode{S} is the scattering transform whose logarithm we want to compute, and \mcode{epsilon} is a regularization constant that prevents the logarithm from going to negative infinity when the scattering coefficients approach zero.

\subsection{Time averaging of a one-dimensional scattering representation}

For several classification problems, it is useful to average scattering coefficients in time as it smoothes their distributions along longer time scales. This is achieved in ScatNet using the \mcode{average_scat} function, whose prototye is the following :

\begin{lstlisting}
  S_av =  average_scat(S,T,step,window_fun);
\end{lstlisting}

where \mcode{S} is the scattering transform we want to average, \mcode{T} is the length of the averaging window, \mcode{step} is the stepping of the successive 
windows (which default value is $T/2$ ), and finally \mcode{window_fun} represents the windowing function to use for averaging. If not provided, \mcode{window_fun} is set to  $@hanning\_standalone$ by default ; this function mimicks the behavior of the \mcode{hanning} window function, without requiring to install the MATLAB Signal Processing Toolbox. 

\subsection{Time aggregation of a one-dimensional scattering representation}

Some signals (mostly audio signals) have a somewhat obvious and particular structure (e.g. the attack sustain and release phases of a muscical note) which we might want 
to provide as a cue to classifiers. In that case, it would totally make sense to concatenate features computed over a certain number of successive time frames. In ScatNet,
the function \mcode{aggregate_scat} empowers you with such an ability. Its  signature is:

\begin{lstlisting}
  S_ag =  aggregate_scat(S,N);
\end{lstlisting}
where \mcode{S} is the scattering transform we want to aggregate and \mcode{N} is the length of the window on which we wish to aggregaate the features.

\subsection{Array-like formatting of a scattering network}

In order to use the scattering coefficients for classification and other tasks, it is often useful to convert the representation into a purely numeric format, separating 
out the metadata. This is the role of the \mcode{format_scat} function, whose signature is
\begin{lstlisting}
[S_table, meta] = format_scat(S, fmt);
\end{lstlisting}
where \mcode{S} is the scattering transform, \mcode{fmt} is the type of output we want, \mcode{S_table} is the formatted output, and \mcode{meta} contain the metadata 
from \mcode{S} associated with the signals in \mcode{S_table}.

There are three different formats available for \mcode{fmt}. The first is \mcode{'raw'}, which leaves the representation intact, setting \mcode{S_table = S}. The second is \mcode{'order_table'}, which creates
 a two-dimensional table for each scattering order, and sets \mcode{S_table} to the cell array of these tables, with \mcode{meta} being a cell array of the associated \mcode{S\{m+1\}.meta}. Finally, \mcode{'table'} (the default), concatenates the tables of all orders and does the same for the meta fields. In the case where \mcode{meta} fields have different number of rows for different orders, the smaller fields are extended with \mcode{-1}. Note that this last format is only possible if scattering coefficients of different orders are of the same resolution, which is most often the case.

\section{Display \label{sec:display}}

\subsection{1D}

\subsubsection{\mcode{scattergram}}

The \mcode{scattergram} function displays the scattering coefficients as a function of time and the last filter index with the remaining filter indices fixed. For first-order coefficients, this means that it is shown as a function of time and $j_1$. For second-order coefficients, it is shown as a function of time and $j_2$ for a fixed $j_1$, and so on.

For one display, a pair consisting of a scattering layer along with a scale prefix is needed. We can specify multiple displays by specifying multiple of these pairs, which will then be displayed one on top of the other.

Its signature is the following
\begin{lstlisting}
img = scattergram(S{m1+1}, jprefix_1, S{m2+1}, jprefix_2, ...);
\end{lstlisting}
where \mcode{S} is the scattering transform, \mcode{m1} is the order of the first display, \mcode{jprefix_1} is the first scale prefix, and so on. Note that is \mcode{m1 = 1} the layer is of first order so no scale prefix is necessary and we set \mcode{jprefix_1} to \mcode{[]}.

Instead of a scattering output \mcode{S}, we can also specify intermediate modulus coefficients \mcode{U}.

\subsubsection{\mcode{display_slice}}

The \mcode{display_slice} function displays the scattering coefficients of two consecutive scattering levels $m$ and $m+1$ (with $m = 1, 2, \ldots$) as a function of $j_m$ and $j_{m+1}$ for a fixed time and fixed $j_1, j_2, \ldots, j_{m-1}$. This display consists of a column that shows the scattering coefficients for all index values $j_m$ at time $t$ (with the fixed $j_1, j_2, \ldots, j_{m-1}$), and a matrix that shows the scattering coefficients for all index pair values $(j_m,j_{m+1})$ (i.e., rows indexed by $j_m$ and columns indexed by $j_{m+1}$) at time $t$ (with the fixed $j_1, j_2, \ldots, j_{m-1}$). Notice the matrix only displays coefficient values for $j_{m+1} > j_m$, and all other coefficients are displayed as $-\infty$ since they should not contain important information according to scattering properties.

The signature of \mcode{display_slice} is the following
\begin{lstlisting}
[sc1,sc2] = display_slice(S,t,scale);
\end{lstlisting}
where \mcode{S} is the scattering transform, \mcode{t} is time, and \mcode{scale} contains the fixed scale path $j_1, j_2, \ldots, j_{m-1}$. By default, when used without providing the \mcode{scale} list (or using an empty list), the function displays the first two levels (i.e., $j_1$ and $j_2$) of the scattering.

\subsection{2D}

The \mcode{image_scat} opens one figure per layer and displays all the path of a layer in the corresponding figure. Its signature is 
\begin{lstlisting}
Sx = scat(x, Wop);
% display scattering coefficients
image_scat(Sx)
\end{lstlisting}

To display all the coefficients of a particular layer, one can use 
\begin{lstlisting}
m = 2;
image_scat_layer(Sx{m+1})
\end{lstlisting}
This stacks all paths of the layer next to each other, renormalize them and overlay their meta information. Both these features can be turned off using
\begin{lstlisting}
m = 2;
image_scat_layer(S{m+1}, 0, 0);
\end{lstlisting}
Paths can also be arranged in a specific order along rows and column.
For example, the following arranges scattering coefficient of second order
by lexicographic order on $j$ along rows and lexicographic order on $\theta$ along columns.
\begin{lstlisting}
var_y{1}.name = 'j';
var_y{1}.index = 1;
var_y{2}.name = 'j';
var_y{2}.index = 2;
var_x{1}.name = 'theta';
var_x{1}.index = 1;
var_x{2}.name = 'theta';
var_x{2}.index = 2;
m = 2;
image_scat_layer_order(S{3},var_x, var_y,1);
\end{lstlisting}
The spatial coefficients of all filters of a filter bank can be displayed using
\begin{lstlisting}
display_filter_bank_2d(filters);
\end{lstlisting}

\section{Classification \label{sec:classification}}

ScatNet contains a classification pipeline for affine space and support vector classifiers. These let the user specify a feature transformation for computing features in batch form, then set up a training/testing mechanism for evaluating results.

\subsection{Batch Computation \label{sec:batch-computation}}

To compute scattering coefficients for a database of signals, we first define a source. A source is a set of files, each file containing a number of objects to be classified. Each object has a class and a location within the file.

The source can be created using \mcode{create_src(directory, extract_objects_fun)} (or one of its wrappers, such as \mcode{gtzan_src}, \mcode{phone_src}, \mcode{uiuc_src}, etc). Given a directory, this function recursively traverses it, looking for `.jpg', `.wav', or `.au' files. For each such file, it calls \mcode{extract_objects_fun} to determine its constituent objects and their classes. To add a new database, it suffices to write the corresponding \mcode{extract_objects_fun} function.

The \mcode{extract_objects_fun} function has the following signature
\begin{lstlisting}
[objects, classes] = extract_objects_fun(file);
\end{lstlisting}
Given a file, it returns a structure array \mcode{objects} of the objects it contains and a cell array of their class names \mcode{classes}. The structure of the \mcode{objects} variable is described below.

Most often, a file will only contain one object, so \mcode{objects} is a  structure of size one, but this can vary from database to database (for example in speech recognition a file corresponds to a recording of a phrase, which can contain different phones to be classified).

Similarly, in musical instrument recognition, you might want for application purpose, to recognise the instrument over a
frame of 2 seconds of duration, when the actual audio file in the database is about 1 min long.  For the creation
of a source in such a context, one could write the following code:

\begin{lstlisting}

function train_src =instrument_sounds_src(directory,T)

if nargin < 2
    T= 2^16;% at a sampling frequency of 32 KHz, this represents a duration of 2s. 
end

src = create_src(directory,@(file)(objects_fun(file,directory,T)));
end

function [objects,classes] = objects_fun(file,directory,T)

tmpsnd = data_read(file);
sig_sz = size(tmpsnd,1);
nb_objs = floor(sig_sz/T);

shortf_ind = file;
path_folders = regexp(directory,filesep,'split');

 objects = struct();
 classes = {};
 path_str = fileparts(file);
 path_parts = regexp(path_str,filesep,'split');
    
  cname = path_parts{end};
   
    if mod(slgth,T) ~= 0
        error('T must divide slgth and T must be a power of 2');
    end
     
 for n = 1:nb_objs
    objects(n).u1 = 1+(n-1)*T;
    objects(n).u2 = n*T;
    classes{n} = cname;
 end

end

\end{lstlisting}

The objects structure has the fields \mcode{objects.u1} and \mcode{objects.u2} which correspond to the lower-left and upper-right corners (start- and endpoints) of the object in a 2D (or 1D) signal. These bounds are inclusive. For the case where a file only contains one object, the bounds will simply be the bounds of the signal. If desired, additional fields can be specified that will be passed on and saved in the source structure.

Given a directory and this function, \mcode{create_src} returns a source structure \mcode{src}, which has the following fields
\begin{itemize}
	\item \mcode{src.classes}: A cell array of the class names.
	\item \mcode{src.files}: A cell array of the file names.
	\item \mcode{src.objects}: A struct array of the objects returned by applying \mcode{objects_fun} to each filename, but with two supplementary fields: \mcode{src.objects.ind}, corresponding to the index of its parent filename in \mcode{src.files} and \mcode{src.objects.class}, corresponding to the index of its class in \mcode{src.classes}.
\end{itemize}

The scattering coefficients of the objects in \mcode{src} can now be calculated using the function \mcode{prepare_database}. This function takes for input the \mcode{src} structure, a function handle \mcode{feature_fun}, which calculates the feature representation, and an options structure. If we want to calculate multiple feature representations and concatenate the results, their feature functions can be assembled into a cell array and given instead of \mcode{feature_fun}.

A ``feature function'' takes as input one or more signals and returns the their feature vectors. Its signature is thus
\begin{lstlisting}
features = feature_fun(signal);
\end{lstlisting}

When \mcode{signal} is a 1D signal, \mcode{features} should be a 2D matrix (potentially just a column vector) where the first dimension is the feature index and the second dimension is a time index. The \mcode{signal} argument can also be a collection of 1D signals, arranged as a columns of a 2D matrix. In this case, the output \mcode{features} is a 3D array, with the first two dimensions being the same, but the third dimension corresponding to the signal index in \mcode{signal}.

Similarly, for \mcode{signal} being a 2D signal, \mcode{features} is a 2D matrix with first dimension corresponding to feature dimension and the second dimension corresponding to a space index. Again, multiple 2D signals can be given to \mcode{feature_fun} by specifying \mcode{signal} to be a 3D array, with signals arranged along the third dimension. The output \mcode{features} will then also be a 3D array as in the 1D case.

Note that a feature function is not required to output only one feature vector for each signal. That is, it can leave a time/space index in the output feature representation. If more than one feature vector is associated to a single signal object, the classifier will classify each feature vector separately and aggregate the results, either through averaging the approximation error (for the affine space classifier) or through voting (for the SVM classifier). If this is not desired, the user has to make sure that the feature function only returns feature vectors with one time index per signal.

The \mcode{scat} function, combined with the \mcode{format_scat} function, fulfills the above criteria. As such, we can use it to define a feature function for \mcode{prepare_database}. For batch calculation of second-order audio scattering vectors for the GTZAN dataset, we have the following code.
\begin{lstlisting}
N = 5*2^17;
src = gtzan_src;

filt_opt = default_filter_options('audio', 8192);
scat_opt.M = 2;

Wop = wavelet_factory_1d(N, filt_opt, scat_opt);

feature_fun = @(x)(format_scat( ...
	log_scat(renorm_scat(scat(x, Wop)))));
	
database_opt = struct();
	
database = prepare_database(src, feature_fun, database_opt);
\end{lstlisting}
Observe the signal length \mcode{N}, which is equal to the length of the signals contained in \mcode{gtzan_src}. 

The \mcode{database_opt} contains options given to \mcode{prepare_database}. These include:
\begin{itemize}
	\item \mcode{database_opt.feature_sampling}: The factor by which the features should be downsampled. Useful for reducing computational complexity of the classifier training and testing (default \mcode{1} - no downsampling).
	\item \mcode{database_opt.file_normalize}: The norm to use when normalizing each file. Can be either \mcode{1}, \mcode{2}, \mcode{Inf} or \mcode{[]}, corresponding to $\mathbf{l}^1$, $\mathbf{l}^2$, $\mathbf{l}^\infty$ norms and no normalization, respectively (default \mcode{[]}).
	\item \mcode{database_opt.parallel}: If \mcode{1}, \mcode{prepare_database} will attempt to parallelize the code by using the \mcode{parfor} construct from the MATLAB Distributed Computing Toolbox. Note that if this is not available, performance may be worse since MATLAB will only use one core (default \mcode{1}).
\end{itemize}

The resulting \mcode{database} structure now contains three fields:
\begin{itemize}
	\item \mcode{database.src}: The original \mcode{src} structure from which the features were computed.
	\item \mcode{database.features}: The feature vectors, arranged in a 2D matrix, each feature vector forming a column.
	\item \mcode{database.indices}: The indices corresponding to each object in \mcode{database.src}. Specifically, \mcode{database.features(:,database.indices\{k\})} contains the features vector(s) calculated from \mcode{database.src.objects(k)}. If there is one feature vector per object, each element of \mcode{database.indices} is a scalar.
\end{itemize}

\subsection{Classifier Framework}
Once the database of feature vectors is constructed, models can be trained and tested. To create a train/test partition, the function \mcode{create_partition} is available. Given a source \mcode{src} and a ratio, it partitions the objects so that each class is divided evenly into the training and testing set according to the ratio. For example, the following code defines a train/test partition with $80\%$ of the objects in the training set:
\begin{lstlisting}
[train_set, test_set] = create_partition(src, 0.8);
\end{lstlisting}
Here, \mcode{train_set} will contain the indices in \mcode{src.objects} that correspond to the training instances, while \mcode{test_set} will contain those corresponding to testing instances.

Alternatively, a partition can be defined manually by traversing the \mcode{src.objects} array and recording the desired indices.

Given a partition, a model can be trained using the appropriate training function. We will denote the training function by \mcode{train}, but in reality it will be either \mcode{affine_train} or \mcode{svm_train}. The training is then done by:
\begin{lstlisting}
model = train(database, train_set, train_opt);
\end{lstlisting}
The structure \mcode{train_opt} contains various parameters for training the model, which will depend on the type of classifier used.

To test the model, the functions \mcode{affine_test} and \mcode{svm_test} are used. Here, we denote the testing function by \mcode{test}. The labels obtained for the testing instances specified by \mcode{test_set} are then obtained by:
\begin{lstlisting}
labels = test(database, model, test_set);
\end{lstlisting}
To calculate the classification error, the \mcode{classif_err} is used:
\begin{lstlisting}
err = classif_err(labels, test_set, src);
\end{lstlisting}
Note that the original \mcode{src} structure is necessary here to verify the class membership of the testing instances.

To calculate the average recognition rate, the \mcode{classif_recog} is used:
\begin{lstlisting}
    mean_recog_rate = classif_recog(labels, test_set, src);
\end{lstlisting}
Note that the original \mcode{src} structure is necessary here to verify the class membership of the testing instances.


To compute the average error rate, the \mcode{classif_mean_err_rate} is used:
\begin{lstlisting}
    mean_err_rate = classif_mean_err_rate(labels, test_set, src);
\end{lstlisting}
Note that, just like for the functions above, the original \mcode{src} structure is necessary here to verify the class membership of the testing instances.


To calculate the errors over a range of training parameters in \mcode{train_opt}, the functions \mcode{affine_param_search} and \mcode{svm_param_search} are available. Denoting them both by \mcode{param_search}, they have the following signature
\begin{lstlisting}
[err, param1, param2] = param_search( ...
	database, train_set, valid_set, train_opt);
\end{lstlisting}

Where \mcode{database} and \mcode{train_set} are the feature database and training set, respectively, \mcode{valid_set} is a validation set to calculate the errors on and \mcode{train_opt} contain the training parameters for the model. Multiple training parameters are specified by setting the field to an array instead of a scalar in \mcode{train_opt}. The classification error for different parameter combinations is then given in \mcode{err}, with the actual parameter values stored in \mcode{param1}, \mcode{param2}, etc.

The validation set can be taken to be \mcode{test_set}, but this would lead to overfitting of the parameters. For some tasks, a separate validation test is given for this purpose. Another alternative is to perform cross-validation on the training set \mcode{train_set}. This is obtained by setting \mcode{test_set} to \mcode{[]}. By default, \mcode{param_search} will then split \mcode{train_set} into five parts, train on four and test on the remaining one, then shuffle and repeat. The number of folds can be determined by the \mcode{train_opt.cv_folds} option. In this case, \mcode{err} is a 2D matrix, with the parameter set index along the first dimension and the fold index along the second. To obtain the best parameter set, an average can then be performed along the fold index.

\subsection{Affine Space Classifier}

The affine classifier is defined by the functions \mcode{affine_train} and \mcode{affine_test}. The former takes as an option the number of dimensions, \mcode{train_opt.dim} used for the affine space model of each class. 

The following code calculates the error for an affine space classifier dimension $40$:
\begin{lstlisting}
train_opt.dim = 40;
model = affine_train(database, train_set, train_opt);
labels = affine_test(database, model, test_set);
err = classif_err(labels, test_set, src);
\end{lstlisting}

For an affine space classifier, the model structure contains the dimension for which it is defined \mcode{model.dim}, the centers of the classes \mcode{model.mu} (a cell array) and the direction vectors defining the affine space \mcode{model.v} (a cell array).

In addition to the labels, \mcode{affine_test} can also output the approximation errors for each feature vector and object. This is done by the additional outputs \mcode{obj_err} and \mcode{feature_err}, as in the following call:
\begin{lstlisting}
[labels, obj_err, feature_err] = affine_test(database, model, test_set);
\end{lstlisting}
Here, \mcode{obj_err} is a 2D matrix with the first dimension corresponding to the class index of the affine space the object was approximated using and the second dimension corresponding to the object index. Similarly, \mcode{feature} is a 2D matrix with first dimension being the class index and the second dimension corresponding to the feature vector index, arranged according to the objects they belong to. That is, the feature vectors of \mcode{test_set(1)} come first, then those of \mcode{test_set(2)}, etc. Note that \mcode{obj_err} can be calculated from \mcode{feature_err} by averaging over the all the feature vectors for one object.

To find the optimal dimension of the affine space, \mcode{affine_param_search} is used, with a range of dimensions specified for \mcode{train_opt.dim}. Specifically, we call
\begin{lstlisting}
[err, dim] = affine_param_search(database, train_set, valid_set, train_opt);
\end{lstlisting}
where \mcode{database}, \mcode{train_set}, \mcode{valid_set}, \mcode{train_opt} are as described in the previous subsection. The output \mcode{dim} is equal to \mcode{train_opt.dim}, and \mcode{err} contains the errors for this range of dimensions.

\subsection{Support Vector Classifier}
ScatNet is compatible with the LIBSVM library (see \url{http://www.csie.ntu.edu.tw/~cjlin/libsvm/}), which provides tools for linear classification based on support vector machines. To execute the operations described here, you will first need to download this free software.

The support vector machine classifier is defined by the functions \mcode{svm_train} and \mcode{svm_test}. The training options consist of:
\begin{itemize}
	\item \mcode{train_opt.kernel_type}: The kernel type used in the SVM. Can be \mcode{'linear'} for a linear kernel $u^Tv$ or \mcode{'gaussian'} for a Gaussian kernel $e^{-\gamma\|u-v\|^2}$.
	\item \mcode{train_opt.C}: The slack factor $C$ used for training the SVM.
	\item \mcode{train_opt.gamma}: The regularity constant $\gamma$ used in the case of a Gaussian SVM kernel.
    \item \mcode{train_opt.reweight}: reweight $\in \{0,1\}$. If reweight is set to one, then for a 2-class SVM classication problem where the two classes are numbered $1$ and $2$, the function to minimize in the svm optimization problem will be :
    \begin{equation*}
        \operatorname*{arg\,min}_{w,b,\xi} \, \frac{1}{2}w^{T}w + \frac{C(N_1+N_2)}{N_1}\sum_{i=1}^{N_1}\xi_i +  \frac{C(N_1 + N_2)}{N_2}\sum_{i=1}^{N_2}\xi_i
    \end{equation*}
    
    where $N_1$ and $N_2$ are the total numbers of training features respectively for class $1$ and class $2$. This remains true for any K-class classification problem with $K \geq 2$. This means
    that an error on the class with the smallest number of features will be penalized more than an error
    on the class with the highest number of features.
  \end{itemize}  
    Sometimes, it happens that the database is very unbalanced, with some of the classes containing far more samples than others. In such case, one of the most convenient way to measure the perfomance of the classifier is the \mcode{average recognition rate} (this is achieved in ScatNet through the function \mcode{classif_recog} which computes the individual recognition rate of each class as well as their mean).   
    In such case we should tell the svm that we want it to give the same importance to the individual error rate of each class and this is achieved by setting the training option \mcode{reweight} to 1.\par The function \mcode{classif_mean_err_rate} is derived from \mcode{classif_recog} and computes the average error rate of a classification experiment.
  

The following code calculates the error for an SVM classifier with a linear kernel and $C = 8$:
\begin{lstlisting}
train_opt.kernel_type = 'linear';
train_opt.C = 8;
model = svm_train(database, train_set, train_opt);
labels = svm_test(database, model, test_set);
err = classif_err(labels, test_set, src);
\end{lstlisting}

If the feature vectors are very large, which can be the case for scattering coefficients, it is often useful to precalculate the kernel. This significantly speeds up training and testing, but at a cost of memory consumption. To precalculate the kernel, the function \mcode{svm_calc_kernel} is used on \mcode{database}, as in
\begin{lstlisting}
database = svm_calc_kernel(database, kernel_type);
\end{lstlisting}
where \mcode{kernel_type} is either \mcode{'linear'} or \mcode{'gaussian'}. Afterwards, when \mcode{svm_train} is called, it can take advantage of this kernel if the same kernel type is specified in \mcode{train_opt.kernel_type}.

To determine the optimal parameters, a parameter search function \mcode{svm_param_search} is provided, with the same syntax as described above. Its signature is
\begin{lstlisting}
[err, C, gamma] = svm_param_search(database, train_set, valid_set, train_opt);
\end{lstlisting}
where \mcode{database}, \mcode{train_set}, \mcode{valid_set}, \mcode{train_opt} are as described previously. We set the values of \mcode{C} and \mcode{gamma} that we wish to test in \mcode{train_opt}. All combinations of these parameters are then evaluated, with errors recorded in \mcode{err} and the particular values for each set stores in \mcode{C} and \mcode{gamma}.

Since parameters can vary over a large domain, and a fine-grained parameter search can be quite costly, there is also an adaptive parameter search function, \mcode{svm_adaptive_param_search}. This function calls \mcode{svm_param_search}, recenters the \mcode{C}-\mcode{gamma} parameter grid on the optimal parameter combination while increasing the grid resolution,  then reapplies \mcode{svm_param_search}, and so on. This is done \mcode{train_opt.search_depth} times, which by default equals $2$.

If a linear kernel is used, we can extract the discriminant vector $w$ and bias $\rho$ from the model using the function \mcode{svm_extract_w}. The function takes as input the database \mcode{database} and the SVM model \mcode{model}. It outputs the $w$s corresponding to each pair of classes in the model, arranged in the order $1$vs$2$,$1$vs$3$,\ldots$1$vs$N$,$2$vs$3$,\ldots,$2$vs$N$,\ldots,$(N-1)$vs$N$, where $N$ is the number of classes. The function is called as in:
\begin{lstlisting}
[w,rho] = svm_extract_w(database, model);
\end{lstlisting}

\section{Separable time and frequency scattering}
\label{sec:separable-tf}

For temporal signals, the regular scattering transform as implemented in \mcode{scat} using operators \mcode{Wop} from \mcode{wavelet_factory_1d} calculates a representation that is locally invariant to time-shifts and stable to time-warping deformation. However, for many tasks, this is not sufficient. Often an additional invariance to frequency transposition and stability to frequency warping is needed. For example, this is the case for speaker-independent phone segment identification, where a single phone can be pronounced at different pitches, and the formant peaks can be deformed in frequency.

A simple way of creating such an invariant representation out of the scattering transform is to average along the log-frequency axis $\lambda_1$. Since a frequency transposition corresponds to a translation along log-frequency, this will create the desired invariance and stability properties. However, it also loses a great amount of information on frequential structure.

To remedy this problem, we apply frequential scattering transform along the log-frequency axis of the original temporal scattering transform. Since the frequency axis is not entirely sampled on a logarithmic scale, but is linear at low frequencies, this is not strictly the same as scattering along log-frequency, but the difference is only visible on the linear regime of the filter bank, which is often below audible frequencies in many applications.

Having computed a temporal scattering transform of a signal \mcode{x} using
\begin{lstlisting}
filt_opt.Q = [8 1];
filt_opt.J = T_to_J(1024, filt_opt);

scat_opt.M = 2;

Wop = wavelet_factory_1d(size(x), filt_opt, scat_opt);

S = scat(x, Wop);
\end{lstlisting}
we now wish to calculate a scattering transform along the \mcode{j(1,:)} index of S, for all orders, which yields the desired invariance property.

To do so, we must first construct a one-dimensional wavelet operator defined on the \mcode{j(1,:)} axis. This is done as in the temporal case, but the signal length now corresponds to the maximum number of \mcode{j(1,:)} samples, which can be found by calling \mcode{length(S{2}.signal)}, which returns the number of first-order coefficients. We thus have
\begin{lstlisting}
Nfreq = length(S{2}.signal);

freq_filt_opt.J = ceil(log2(Nfreq));
freq_scat_opt.M = 1;

Wop_freq = wavelet_factory_1d(Nfreq, freq_filt_opt, freq_scat_opt);
\end{lstlisting}
The options structures \mcode{freq_filt_opt} and \mcode{freq_scat_opt} have the same parameters and function as their temporal counterparts \mcode{filt_opt} and \mcode{scat_opt}.

Note that \mcode{freq_filt_opt.J} now controls the degree of transposition invariance. Since \mcode{freq_filt_opt.Q = 1} by default, we have that the width of the averaging filter $\phi$ is given by $2^J$. Setting \mcode{freq_filt_opt.J} to \mcode{ceil(log2(Nfreq))} thus guarantees total transposition invariance. The degree of invariance necessary depends on the task at hand. If only one octave of invariance is needed, we can specify
\begin{lstlisting}
freq_filt_opt.J = filt_opt.Q;
\end{lstlisting}
where \mcode{filt_opt.Q} is the number of wavelets per octave in the orignal, temporal, filter bank.

Finally, we calculate the frequential scattering transform by calling the \mcode{scat_freq} function.
\begin{lstlisting}
[S_freq, U_freq] = scat_freq(S, Wop_freq);
\end{lstlisting}
Just like \mcode{S}, \mcode{S_freq} is a cell array of scattering layers, each layer corresponding to the original temporal scattering order. However, each temporal scattering layer now contains the coefficients of the frequential scattering transform applied to that layer. That is, \mcode{S_freq{2}} contains all the frequential scattering coefficients of \mcode{S{2}}, while \mcode{S_freq{3}} is the transform of \mcode{S{3}}.

The identity of the coefficients is given by the meta fields. For example \mcode{meta.fr_order} gives the order of the frequential scattering transform, while \mcode{meta.fr_j} gives the frequential scale, corresponding to the ``quefrency'' of the frequential filter. For example,
\begin{lstlisting}
coeff = S{2}.signal(S{2}.meta.fr_order==0);
\end{lstlisting}
extracts the zeroth-order frequential scattering coefficients of the first-order temproal coefficients. From the definition of the transform, these are simply the original temporal scattering coefficients averaged along log-frequency. The temporal scale of these coefficients is still stored in \mcode{meta.j}, so we can write
\begin{lstlisting}
j = S{2}.meta.j(:,S{2}.meta.fr_order==0);
\end{lstlisting}
to obtain the scales \mcode{j} corresponding to \mcode{coeff}. Note that since the scattering transform subsamples to the critical  bandwidth, the coefficients have been subsampled along frequency, which is reflected in \mcode{j}.

We can also extract the first-order frequential coefficients and the corresponding temporal scales \mcode{j} for a particular frequency scale \mcode{fr_j}, using the following code
\begin{lstlisting}
coeff = S{2}.signal(S{2}.meta.fr_order==1 & ...
   S{2}.meta.fr_j==0);
j = S{2}.meta.j(:,S{2}.meta.fr_order==1 & ...
   S{2}.meta.fr_j==0);
\end{lstlisting}
And finally, we can extract the same frequential coefficients, but for the second temporal order scattering coefficients. In this case, we also need to specify \mcode{meta.j(2,:)}. The code is then
\begin{lstlisting}
coeff = S{3}.signal(S{3}.meta.fr_order==1 & ...
   S{3}.meta.fr_j==0 & ...
   S{3}.j(2,:)==9);
j = S{3}.meta.j(:,S{3}.meta.fr_order==1 & ...
   S{3}.meta.fr_j==0 & ...
   S{3}.j(2,:)==9);
\end{lstlisting}

\bibliographystyle{plain}
\bibliography{doc_references}

\end{document}