% This file will cover some aspects of time-discrete signals, filtering etc

\chapter{Time-discrete signals}
Time-discrete signals are typically used to approximate time-continuous
signals, and are commonly found in modern digital signal handling.

Typically, a time-continuous signal $f(t)$ is sampled at discrete, equally
spaced times $t$ that are separated by a fixed sampling interval $\Delta t$:
\begin{equation}\label{eq:discretization}
  f(t)\rightarrow f_k = f(t=k \Delta t)\,,\quad k\in \mathbb{Z}\,.
\end{equation}
The signal is now represented by a discrete sequence of numbers $f_k$, that
map the original function in are ``reasonable'' manner. Intuitively, if the
interval $\Delta t$ is chosen small enough, the features in the original
signal $f(t)$ is captured sufficiently well. The sequence of numbers $f_k$ is
an infinite set of numbers, in the same way as the time argument $t$ in the
original function is not limited.

\begin{figure}
 \centerline{ \includegraphics{timediscrete/sampling1}}
  \caption{A time-continuous function $f(t)$ gets sampled at discrete times
    $t=k\Delta t$, leading to samples $f_k$.}
\end{figure}

\section{Sampling and the Dirac comb function}
A useful concept to formalize this transition from a continuous to a sampled
signal makes use of an array of delta functions spaced in equal time intervals
$\Delta t$. We refer to this as the Dirac comb function
\begin{equation}\label{eq:combfunction}
 \comb_{\Delta t}(t)=\sum\limits_{k=-\infty}^{\infty}\delta(t-k\Delta t)\,.
\end{equation}

\begin{figure}
  \centerline{ \includegraphics{timediscrete/sampling2}}
  \caption{\label{fig:samplerepresentation}
    Multiplication of the original function $f(t)$ with a Dirac comb
    $\comb{\Delta t}(t)$ leads to a representation $f_R(t)$ of the sampled
    function that can be used to understand its spectral properties.}
\end{figure}

To represent the time-discrete sampled signal in a time-continuous way, we
can use a similar concept to the comb function above, but modify the weight of
each $\delta$ peak. This can be done by simply multiplying $f(t)$ by the comb
function $\comb_{\Delta t}(t)$,
\begin{equation}\label{eq:samplingexpression}
  f_R(t)=f(t)\cdot \comb_{\Delta t}(t)=\sum\limits_{k=-\infty}^{\infty}f_k \delta(t-k \Delta t)\,,
\end{equation}
as the Kronecker-$\delta$s in $\comb_{\Delta t}(t)$ sample the original
function at the correct positions (see Fig.~\ref{fig:samplerepresentation}).

This new function $f_R(t)$ has, for $T\gg\Delta t$ and for a
sufficiently smooth original function $f(t)$, approximately the same 
integral over $T$ as the original function:
\begin{equation}
\int\limits_{0}^{T}\! f(t)dt \approx
\int\limits_{0}^{T}\!f_R(t)dt
=\int\limits_{0}^{T}\!\sum\limits_{k=-\infty}^{\infty}f_k \delta(t-k \Delta t)dt
\approx\Delta t\sum\limits_{k=0}^{T/\Delta t}f_k\,.
\end{equation}

\subsection{Spectral properties of sampled functions}
To understand quantitatively the effect of sampling, it is helpful to
compare the spectral content of $f_R(t)$ to the one of the continuous signal
$f(t)$, given by the Fourier-transformed signal:
\begin{equation}
  \tilde{f}(\omega)=
  {\cal{F}}[f(t)]=\int\limits_{-\infty}^{+\infty}\!f(t)e^{-i\omega t}\,dt\,. 
\end{equation}
%Strictly speaking, this expression only holds if $f(t)$ is a square-integrable
%function, but most practical signals can be made to have this property. 
The Fourier transformation of $f_R(t)$ can then be calculated (see details in
appendix \ref{sec:FTdiscrete}):
\begin{equation}\label{eq:discreteFourier}
  \tilde{f}_R(\omega)=\sum\limits_{k=-\infty}^{\infty}f_k\,e^{-i\omega \Delta t\cdot k}\,,
\end{equation}
This is a  {\em continuous} function in $\omega$, even though the underlying
time series $f_k$ is discrete; this is referred to as a Fourier series. As
such,  $\tilde{f}(\omega)$ is a periodic function:
\begin{equation}
  \tilde{f}(\omega)=\tilde{f}(\omega+k\cdot2\pi/(\Delta
  t))\quad\text{for}\quad k\in\mathbb{Z}\,.
\end{equation}

A special case of such a Fourier series is one where $f_k=1$ for all $k$,
corresponding to $f_R(t)=\comb_{\Delta t}(t)$, the Dirac comb. In this case,
the result is
\begin{equation}
  {\cal F}\left[\comb_{\Delta t}(t)\right]=\sum\limits_{k=-\infty}^{\infty}\,e^{-i\omega \Delta
                       t\cdot k}
\end{equation}
The Dirac comb function can be represented by a Fourier series,
%{\bf likely correct pre-factor}
\begin{equation}
  \comb_{\Delta t}(t)=\sum\limits_{k=-\infty}^{\infty}\delta(t-k\Delta t)=
  \frac{1}{\Delta t}\sum\limits_{n=-\infty}^{+\infty}e^{i\frac{2\pi}{\Delta
      t}nt}\,,
 \end{equation}
so the Fourier transform of the Dirac comb function is again a Dirac
comb. Proper care of pre-factors (see appendix \ref{appendix:DiracComb}) gives

\begin{equation}
  {\cal F}[\comb_{\Delta t}(t)] = \omega_S
  \sum\limits_{k=-\infty}^{\infty}\delta(\omega-k\omega_S)
  \quad\text{with}\quad \omega_S=\frac{2\pi}{\Delta t}\,.
\end{equation}

\begin{figure}
  \centerline{\includegraphics{timediscrete/dirac_comb.pdf}}
  \caption{\label{fig:diraccomb} A Dirac comb $\comb_{\Delta t}(t)$ in time
    with a separation of $\Delta t$ between the peaks  has a frequency
    spectrum that is also a Dirac comb, with a frequency separation $\omega_S$
  between the different components.}
\end{figure}

This allows for a useful interpretation of the spectral content of a sampled
system presented in (\ref{eq:discreteFourier}), making use of convolution
theorem (\ref{eq:convolution1}):
%{\bf Doubts about pre-factor}
\begin{eqnarray}
  \tilde{f}_R(\omega)&=&{\cal{F}}[f(t)\cdot\comb_{\Delta t}(t)] =
                         \tilde{f}(\omega)\ast{\cal F}[\comb_{\Delta t}(t)](\omega)
                         \nonumber\\
                     &=&\int\limits_{-\infty}^{\infty}\!\tilde{f}(\omega')\cdot
                         {\cal F}[\comb_{\Delta t}t)](\omega-\omega')d\omega'
                         \nonumber\\
                     &=&\int\limits_{-\infty}^{\infty}\!\tilde{f}(\omega')\left[
                         \omega_S\!\sum\limits_{k=-\infty}^{\infty}\!\delta
                         (\omega-\omega'-k\omega_S)\right]d\omega'\nonumber\\
                     &=&\omega_S
                         \!\sum\limits_{k=-\infty}^{\infty}
                         \tilde{f}(\omega-k\omega_S)\label{eq:sampledspec1}
\end{eqnarray}
This means that the spectrum $\tilde{f}_R(\omega)$ of a (reconstructed)
sampled signal  contains a superposition of copies the original signal
spectrum, $\tilde{f}(\omega)$, shifted by integer multiples of the sampling
(angular) frequency $\omega_S$. This has two important consequences, namely a
relationship between useful bandwidth of a signal and the sampling
frequency, and a method of shifting frequency components simply by sampling.

\subsection{Sampling theorem}
We first return to the intuitive insight that, in order to capture the
essential features of a continuous signal $f(t)$ in a sampling process, the
sampling time interval $\Delta t$ has to be significantly shorter than the
time scale of changes in the original signal $f(t)$.
To quantify the smoothness of the original signal, we assume that the
spectral amplitudes $\tilde{f}(\omega)$ decrease significantly above a certain
frequency $\omega_C$, referred to as a cutoff frequency or the bandwidth of a
signal. Such a situation is e.g. found in the range of audible signals, where
the bandwidth corresponds to a few kHz, but just about any signal that is
processed in some way has such an upper frequency $\omega_C$.

\begin{figure}
  \centerline{\includegraphics{timediscrete/baseband.pdf}}
  \caption{\label{fig:baseband}Fourier distributions $\tilde{f}(\omega)$ (or
  spectra) of typical signals, with a representative cutoff frequency
  $\omega_C$, indicating a region $-\omega_C\ldots\omega_C$ in which most of
  the signal is contained. The left sample indicates a realistic spectrum, the
  right side a symbolic distribution as often used when representing
  signals. The distributions are symmetric around $\omega=0$, a property of
  all real-valued signals $f(t)$. Both distributions would be a assigned a
  bandwidth $B=\omega_c/2\pi$ (when referenced in frequencies rather than
  angular frequencies).}
\end{figure}

We now consider a situation where a bandwidth-limited signal $f(t)$ is
sampled with a time interval $\Delta t$. The condition of sampling with a
sufficiently high rate compared to the time scale of changes in $f(t)$ is
reflected in a condition $\omega_S\gg\omega_C$, where $\omega_S=2\pi/\Delta t$
is the sampling (angular) frequency, and $\omega_C$ a cutoff frequency (or
bandwidth) of the bandwidth-limited signal.

The resulting spectrum of the reconstructed sampled signal $f_R(t)$ through
(\ref{eq:sampledspec1}) now is a sum of copies of the original spectrum
$\tilde{f}(\omega)$, separated by multiples of the sampling frequency
$\omega_S$ (see Figure~\ref{fig:sampledspectrum}).

\begin{figure}
  \centerline{\includegraphics{timediscrete/sampledspec1.pdf}}
  \caption{\label{fig:sampledspectrum}
    Spectrum $\tilde{f}_R(\omega)$ of the reconstructed sampled
  signal, sampled with $\omega_S$. Copies of the baseband spectrum
  $\tilde{f}(\omega)$ appear at integer multiples of $\omega_S$.}
\end{figure}

To reconstruct the original spectrum $\tilde{f}(\omega)$ from
$\tilde{f}_R(\omega)$, a low pass filter can be applied to the sampled
spectrum to remove the higher order copies of the baseband spectrum. This
works as long as spectral components of the first order copy around $\omega_S$
does not overlap with the baseband spectrum. Quantitatively, this means that
\begin{equation}\label{eq:samplingtheorem}
  \omega_C<\omega_S/2\,,
\end{equation}
or in other words, the highest frequency $\omega_C$ in the signal has to be
smaller than half the sampling frequency $\omega_S$ in order to be uniquely
reconstructible. This relationship is referred to as the {\em sampling
  theorem} or {\em Nyquist theorem}.
An example is the time-discrete sampling of an audio signal in the
CD standard, which uses $\omega_S=2\pi\cdot44.1$\,kHz, which in theory would
be able to reproduce audio signals with an upper frequency around 20\,kHz, the
perception limit of most (young) humans.

\subsection{Nyquist zones}
For a given sampling frequency $\omega_S$, the whole spectrum can be segmented
into different zones, separated by multiples of half the sampling
frequency (see Figure~\ref{fig:nyquistzones}). Then, the sampling theorem
states that in order for a sampled signal to be uniquely reconstructible, its
spectrum has to reside in the first ``Nyquist'' zone. Higher order zones,
e.g. the second zone covering the interval $[\omega_S/2,\omega_S]$ and its
negative counterpart $[-\omega_S,-\omega_S/2]$ contain copies of the original
signal spectrum. For even zones, they contain the mirrored original spectrum,
while the odd zones contain the original spectrum.

\begin{figure}
  \centerline{\includegraphics{timediscrete/nyquistzones.pdf}}
  \caption{\label{fig:nyquistzones}Segmentation of a spectrum in Nyquist zones
    spaced by half the sampling frequency $\omega_S$. Copies of the original
    spectrum appear in the odd zones, whereas copies of mirrored original
    spectra appear in the even zones.}
\end{figure}


\subsection{Frequency shifting by sampling}
The appearance of several copies of the original spectrum $\tilde{f}(\omega)$
in a sampled signal at multiples of the sampling frequency can be conveniently
used to shift the signal of a spectral band. Modern software-defined radio
receivers and similar signal processors make extensively use of this idea, as
a narrow spectral segment at a high carrier frequency can be converted to a
lower frequency band by adequate sampling.

If the spectrum of a real-valued signal $f(t)$ is contained in a small
band around a carrier frequency $\omega_0$, its Fourier transform is actually
non-zero in two bands at positive and negative frequencies.
%
The sampling process then generates multiple copies of both bands at positive
and negative frequencies (see figure~\ref{fig:freqshift}).

\begin{figure}
  \centerline{\includegraphics{timediscrete/freqshift.pdf}}
  \caption{\label{fig:freqshift}A signal with a spectrum centered around a
    frequency $\omega_0$, when sampled with frequency $\omega_S$, leads to a
    reconstructed spectrum $|\tilde{f}_R(\omega)|$ with a frequency spectrum
    shifted by multiples of $\omega_S$. In the example shown above, the
    spectrum centered around $\omega_0$ in the 3. Nyquist zone gets moved to a
    new center frequency $\omega_0-\omega_S$ in the first zone.}
\end{figure}

Now, the condition for unique reconstructability requires that the signal has
to be contained in a single Nyquist zone to avoid overlap with other mirrored
or displaced copies of the spectrum generated in the sampling process.

A restriction of the reconstructed sampled signal to the first Nyquist zone,
as an appropriate filtering of a time-discrete filtering will provide, will
lead to a spectrum of the reconstructed sampled signal $\tilde{f}_R(\omega)$
that is frequency-shifted from its original center frequency by multiples of
the sampling frequency.

When the shifting is from an even Nyquist zone to the first one, the frequency
content gets mirrored. There, higher frequency components in the sampled signal
closer to a multiple of the sampling frequency will be transferred closer to
frequency 0, i.e., to lower frequencies. When the shifting is from an odd
Nyquist zone to the first zone, the complete spectrum is just shifted by
multiples of $\omega_S$.

To come back to the practical application of a higher-order Nyquist sampling
in a software-defined radio, the sampling of the original signal should happen
in a way that is well-represented by a multiplication of the original signal
by a Dirac $\delta$ function according to (\ref{eq:samplingexpression}). In
practice, this means that the sampling mechanism should be significantly
shorter than the sampling interval $\Delta t$.

A consequence of this process is that, if a time-continuous signal $f(t)$ has
spectral components in higher Nyquist zones, their spectral content will be
mapped into the first Nyquist zone, and interfere with the spectral content
there. Therefore, contributions form other Nyquist zones than that of interest
need to be suppressed before sampling.

\section{Time-discrete filters}
Filters are among the most important signal processing elements. To formally
describe these elements for time-discrete signals, it is useful to recall how
conventional (i.e., time-continuous) filters are described. Filters transform
a signal $f(t)$ into another signal $g(t)$. They are are
typically time-invariant in their action,
\begin{equation}
  f(t+T)\rightarrow g(t+T)\quad\text{for all}\quad T\,,
\end{equation}
linear in their response,
\begin{equation}
  a\cdot f(t)\rightarrow a\cdot g(t)\,,
\end{equation}
which implies that sums of signals are transformed into the sum of their
filtered versions, and can most conveniently described as a
multiplication with a frequency-dependent gain $h(\omega)$ in the Fourier
domain:
\begin{equation}\label{eq:filteractioncont}
  \tilde{g}(\omega)=\tilde{h}(\omega)\cdot\tilde{f}(\omega)\,.
\end{equation}
Using the convolution theorem, the filter action can also be written as
\begin{equation}\label{eq:filteractionconvolut}
  g(t)=h(t)\ast f(t)= \int\limits_{-\infty}^{\infty}\! h(t-t')\cdot f(t')dt'\,,
\end{equation}
with $h(t)={\cal F}^{-1}[\tilde{h}(\omega)]$ referred to the {\em impulse
  response} of the filter. It is also referred to as the {\em Green function}
of a filter.

Even for real-valued signals $f(t)$ leading to real-valued results $g(t)$, the
filter gain $\tilde{h}(\omega)$ can be complex-valued. However, the impulse
response must obey
\begin{equation}
  \tilde{h}^*(\omega)=\tilde{h}(-\omega)\,.
\end{equation}

Another common property of filters implemented in real time is {\em
  causality}, namely, that the response of the filter at any given time can
not depend on signals in the future, the impulse response can only have
contributions for positive times:
\begin{equation}
  h(t)=0\quad\text{for}\quad t<0\,.
\end{equation}

To describe filters for time-discrete signals, one one can simply replace the
time-continuous signals $f(t)$ and $g(t)$ in the expressions above with
weighted Dirac combs (\ref{eq:samplingexpression}). As the Fourier transforms
$\tilde{f}(\omega), \tilde{g}(\omega)$ are then periodic in $\omega$, the
filter gain $\tilde{h}(\omega)$ must be as well,
\begin{equation}
\tilde{h}(\omega+k\omega_S)=\tilde{h}(\omega)\,,\quad k\in\mathbb{Z}\,,
\end{equation}
but the filter action description (\ref{eq:filteractioncont}) remains unchanged.
The periodicity of $\tilde{h}(\omega)$ implies that the impulse response
$h(t)$ of the filter can also represented by a time-discrete stream $h_k$,
\begin{equation}
  h(t)=\sum\limits_{k=-\infty}^{\infty}h_k \delta(t-k \Delta t)\,.
\end{equation}

In time domain, the action of a filter on a can be described as
\begin{eqnarray}
  g(t)&=&h(t)\ast f(t)=\int\limits_{-\infty}^{\infty}h(t-t')\cdot f(t')dt'
          \nonumber\\
      &=&\int\limits_{-\infty}^{\infty}\left[\sum_{k=-\infty}^{\infty}\!h_k
          \delta(t-t'-k\Delta t)\right]\cdot
          \left[\sum_{l=-\infty}^{\infty}\!f_l\delta(t'-l\Delta t)\right]dt'
          \nonumber\\
      &=&\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}h_k f_l
          \delta(t-(l+k)\Delta t)\nonumber\\
  &=&\sum_{k=-\infty}^{\infty}\sum_{u=-\infty}^{\infty}h_k
      f_{u-k}\delta(t-u\Delta t)\nonumber\\
  &=&\sum_{u=-\infty}^{\infty}\underbrace{\sum_{k=-\infty}^{\infty}h_kf_{u-k}}_{=:g_u}\delta(t-u\Delta t)\,.
\end{eqnarray}
So the convolution for time-discrete signals takes simply the form of a sum,
\begin{equation}\label{eq:discretefiltertimedomain}
  g_k = \sum_{l=-\infty}^{\infty}h_{l}f_{k-l}\,.
\end{equation}
For causal filters, $h(l<0)=0$, so that the summation only need to be carried
out over positive $l$.

This form also  suggests to distinguish two classes of filters. If the sum in
(\ref{eq:discretefiltertimedomain}) has to be carried over all $l$, the filter
shows an {\em infinite impulse response}. If coefficients $h_l$ vanish for
some index $l>L$, the filter has a {\em finite impulse response}. 

\subsection{Time delay blocks}
\begin{figure}
  \centerline{\includegraphics{timediscrete/storage.pdf}}
  \caption{\label{fig:storageelement}Storage element, implemented by a single
    register that stores a discrete signal $f_k$ each sampling cycle, leading
    to a delayed copy $g_k$. This can be represented in Fourier domain by
    multiplication with a complex value $z^{-1}$.}
\end{figure}
  
A basic element for implementing time-discrete filters is a storage
element. It is typically implemented by a register, clocked by the sampling
frequency. It's action on the discrete input stream $f_k$ can be described as
\begin{equation}
  g_k = f_{k-1}\,.
\end{equation}
Its action in Fourier domain can simply be calculated from this defining
expression:
\begin{eqnarray}
  \tilde{g}(\omega)&=&\sum_{k-=-\infty}^{\infty}g_ke^{-i\omega\Delta tk}
                       =\sum_{k-=-\infty}^{\infty}f_{k-1}e^{-i\omega\Delta
                       tk}\nonumber\\
                   &=&\sum_{k-=-\infty}^{\infty}f_{k-1}e^{-i\omega\Delta
                       t(k-1)}e^{-i\omega \Delta t}\nonumber\\
                   &=&z^{-1}\cdot\tilde{f}(\omega)\,,\quad\text{with}\quad
                       z=e^{i\omega \Delta t}
\end{eqnarray}
A time delay by $\Delta t$, here one sampling cycle, can be represented by a
filter function $\tilde{h}(\omega)=z^{-1}$, a pure phase factor with a phase
linear in $\Delta t$. This is exactly the same filter function that represents
a time delay in time-continuous signals.  As the corresponding impulse
response in time domain is
\begin{equation}
  h_k=\delta_{k,1}\,
\end{equation}
the delay element belongs to the class of finite impulse response filters.

\subsection{Floating average filter}
Another building block consists of a floating average filter. For an averaging
over $N$ samples, its formal defining action can be written as
\begin{equation}
  g_k=\frac{1}{N}\sum_{n=0}^{N-1}f_{k-n}\,,
\end{equation}
corresponding to its impulse response function
\begin{equation}
  h_k=\begin{cases}1/N,&0\le k< N\,,\\0&\text{otherwise}.\end{cases}
\end{equation}
The frequency-dependent filter gain $\tilde{h}(\omega)$ can be directly
calculated from this by Fourier transformation,
\begin{equation}\label{eq:transferaverage}
  \tilde{h}(\omega)=\sum_{n=0}^{N-1}\frac{1}{N}e^{i\omega\Delta t n} =
  \frac{1}{N}\frac{1-e^{-i\omega\Delta t\cdot N}}{1-e^{-i\omega\Delta t}}\,,
\end{equation}
were the key step is simply using an expression for a geometric sum. It is
instructive to re-write this expression:
\begin{equation}
  \tilde{h}(\omega)=\frac{1}{N}e^{-i\omega\Delta
    t(N-1)/2}\frac{\sin(\omega\Delta T N/2)}{\sin(\omega\Delta T/2)}\,.
\end{equation}
A global phase factor, linear in frequency, corresponding to a time delay of
half of the averaging time, followed by a real-valued ration of two sine
functions in frequency. The latter provides strict zeros at frequencies
$\omega_z=2\pi k/(N\cdot\Delta t)$ for $k>0$. In the limit of $\Delta
t\rightarrow 0$, the latter ratio converges against a sinc function, as
expected for the Fourier transformation of a square pulse. The frequency
response of the modulus of the filter is shown in figure~\ref{fig:sincfilter}.
The filter gain $\tilde{h}(\omega)$ has the characteristic of a low pass
filter, reflecting the intuitive impression that averaging over several
samples suppresses faster changes in the sampled signal.

\begin{figure}
  \centerline{\includegraphics{timediscrete/sincfilter.pdf}}
  \caption{\label{fig:sincfilter}Filter transmission of a moving average
    filter with $N$ averaging points for a sampled signal, in linear and
    logarithmic scales. The filter transmission has several zeros, and is
    symmetric in the first and second Nyquist zone.}
\end{figure}


Using the notation $z=e^{i\omega\Delta t}$, the transfer function can also be
written in a compact form that is often found in signal processing literature:
\begin{equation}
  \tilde{h}(\omega)=\frac{1}{N}\sum_{n=0}^{N-1}z^{-n}=\frac{1}{N}\frac{1-z^{-N}}{1-z^{-1}}
\end{equation}
In principle, the discrete coefficients $h_k$ of the impulse response of a
filter with a desired frequency-dependent gain $\tilde{h}(\omega)$  could be
obtained by an inverse  Fourier transformation of the periodic
$\tilde{h}(\omega)$:
\begin{equation}\label{eq:discreteIFT}
  h_k=\frac{1}{2\pi}\int\limits_{\omega_S/2}^{\omega_S/2}\!\tilde{h}(\omega)
  e^{ik\omega\Delta t}\,,\quad\text{for}\quad k\in \mathbb{Z}\,.
\end{equation}
However, such a filter is not automatically causal, nor has it a finite impulse
response. Even the standardized filters (see ...) have an infinite
response.
\subsection{Filters with finite impulse response (FIR)}
Implementing a filter numerically through direct convolution of a signal with
a impulse response through (\ref{eq:discretefiltertimedomain}) requires a
finite length $L$, 
\begin{equation}
  g_k = \sum_{l=-\infty}^{L-1}h_{l}f_{k-l}\,.
\end{equation}
Finding a finite set $\{h_l\}$ for a desired filter characteristic is a complex
process, and can often be only solved numerically and approximately.

One consideration is that coefficients $h_l$ obtained by
(\ref{eq:discreteIFT}) are typically symmetric around $l=0$. Truncating the
impulse response to $-L/2<k<L/2$ still results in an acausal filter, so the
impulse response function typically gets shifted by $L/2$, corresponding to a
time delay of $L\Delta t/2$.

(More to come...)

\subsection{Filters with infinite impulse response (IIR)}
A filter with infinite impulse response that is very easy to implement is a
first order low pass filter with a cutoff frequency $\omega_C$, with a transfer function
\begin{equation}
  \tilde{h}(\omega)=\frac{1}{1+i\omega/\omega_C}
\end{equation}
with impulse response for a continuous time filter 
\begin{equation}
  h(t)=\begin{cases}e^{-\omega_C t}\,,&t>0\,,\\0\,,&t\le0\end{cases}\,.
\end{equation}
When sampled at intervals $t=k\Delta t$, the response function becomes a
geometric series,
\begin{equation}
  h_k\propto a^k\,,\quad\text{with}\quad a=e^{-\omega_C\Delta t}=
  e^{-\omega_C\frac{2\pi}{\omega_S}}\quad\text{for}\quad k>0\,.
\end{equation}
Such a response can be easily implemented by subtracting a fraction of the
content of a register from itself every time step $\Delta t$, with an
iteration rule $h_k=a\cdot h_{k-1}$. The iteration rule can be easily expanded
to add a sampled signal $f_t$, leading to a recursion relation implementing
the complete filter:
\begin{equation}
  g_k = a\cdot g_{k-1}+f_k\,,\quad\text{with}\quad0<a<1
\end{equation}
For small cutoff frequencies $\omega_C\ll\omega_S$,  $a=1-\epsilon$ is close
to 1, with $\epsilon\approx2\pi\omega_C/\omega_S$. 
This recursion relation can be properly normalized to a unit stationary gain:
\begin{equation}
  g_k = a\cdot g_{k-1}+\frac{1}{1-a}f_k=(1-\epsilon)g_{k-1}+\epsilon f_k\,.
\end{equation}
A formal implementation, and its Fourier domain signal flow, is shown in
figure \ref{fig:lowpass1scheme}.
\begin{figure}
  \centerline{\includegraphics{timediscrete/lowpass1.pdf}}
  \caption{\label{fig:lowpass1scheme}Symbolic implementation of a first order
    low pass filter as an example of an infinite impulse response (IIR)
    filter, with its associated symbolic signal flow scheme in the Fourier
    domain.} 
\end{figure}

Such a filter is particularly efficient to implement for coefficients
$\epsilon=2^{-k}$, as multiplication by such a coefficient can be implemented by
a simple bit shift (see figure~\ref{fig:lowpass1impl}). The implementation
requires only bit shifts (a cheap resource), two additions, and a
register. The normalizing factor in figure~\ref{fig:lowpass1impl} has been
moved to the output, but can be omitted when working with integers. As the
output is extracted after the  storage element, it also adds an additional
delay by one sampling time. 
\begin{figure}[!b]
  \centerline{\includegraphics{timediscrete/lowpass1impl.pdf}}
  \caption{\label{fig:lowpass1impl}Implementation of a first order low
    pass filter with coefficients $\epsilon=2^{-k}$. } 
\end{figure}
The nominal cutoff frequency is
\begin{equation}
  \omega_C=-\frac{\omega_s}{2\pi}\ln(1-2^{-k}) \approx \frac{2^{-k}}{2\pi}\omega_S\,.
\end{equation}

The exact filter gain $\tilde{h}(\omega)$ for such a time-discrete filter can be
obtained from the recursion equation in the Fourier domain from
figure~\ref{fig:lowpass1scheme}:
\begin{equation}
  \tilde{g}=(1-\epsilon)\cdot z^{-1}\tilde{g}+\epsilon\tilde{f}\,,
  \quad\text{or}\quad
  \tilde{h}(\omega)=\frac{\epsilon}{1-(1-\epsilon)\cdot
    z^{-1}}\quad\text{with}\quad z^{-1}=e^{-i\omega\Delta t}
\end{equation}

The filter gain and phase for a few values of $k$ is shown in a Bode plot
figure~\ref{fig:lowpass1gain}.
%
\begin{figure}
  \centerline{\includegraphics{timediscrete/firstlowpass.pdf}}
  \caption{\label{fig:lowpass1gain}Bode diagram of the filter gain
    $\tilde{h}(\omega)$ of the first order low pass filter with feedback
    coefficients $\epsilon=2^{-k}$.} 
\end{figure}
%
While for the lower frequency parts, the filter gain is very similar to a
corresponding continuous time first order low pass filter, with a unity gain
at low frequencies, and a slope in the stop band of 6\,dB/octave. However,
there are a few notable differences:
\begin{itemize}
  \item The gain magnitude is minimal at half of the sampling frequency, and
    increases again in the second Nyquist zone. This is reflecting the
    periodicity in the filter gain $\tilde{h}(\omega)$. If the gain is plotted
    linear in frequency, it appears mirrored at half the sampling frequency.
  \item The phase lag of the filtered signal behaves significantly different
    from its continuous-time counterpart: the phase shift becomes
    exactly zero at half of the sampling frequency, and then reverses sign
    in the second Nyquist zone. This is because the second Nyquist zone
    reflects the negative frequency part of the first zone, where the phase
    shift changes sign. And unless the cutoff frequency is significantly
    smaller than the sampling frequency, it does not reach the asymptotic high
    frequency limit of -90$^\circ$.  
  \end{itemize}

(More to come on higher order filters)
\section{Time-discrete mixers}
In time-continuous signal processing, the frequency spectrum of a signal
$f(t)$ can be shifted by multiplying a signal with a ``local'' harmonic
oscillation at a frequency $\omega_L$: 
\begin{equation}
  g(t)=f(t)\cdot m(t)\quad\text{with}\quad m(t)=\cos(\omega_Lt)
\end{equation}
The effect of this process on the spectrum can be understood by directly
calculating the spectrum of the product using the convolution theorem:
\begin{equation}
  \tilde{g}(\omega)={\cal F}\left[g(t)\right] = {\cal F}\left[f(t)\cdot
    m(t)\right] = \frac{1}{2\pi}
  \left[ \tilde{f}(\omega)\ast\tilde{m}(\omega)\right]\,.
\end{equation}
With the spectrum of the ``local oscillator'',
%{\bf (Having doubts about 2 pi)}
\begin{equation}
  \tilde{m}(\omega)={\cal F}
  \left[\frac{1}{2}(e^{i\omega_Lt}+e^{-i\omega_Lt})\right]=
  \pi\left(\delta(\omega-\omega_L)+\delta(\omega+\omega_L)\right)\,,
\end{equation}
this results in a simple expression for the spectrum of the product,
%{\bf(more doubts about 2pi)}
\begin{eqnarray}
  \tilde{g}(\omega)&=&\frac{1}{2\pi}\left[\tilde{f}(\omega)\ast
  \pi\left(\delta(\omega-\omega_L)+\delta(\omega+\omega_L)\right)\right]
  \nonumber\\
  &=&  \frac{1}{2}\left(\tilde{f}(\omega-\omega_L)+\tilde{f}(\omega+\omega_L)\right)\,,
\end{eqnarray}
which is just a sum of two copies of the original spectrum, displaced by the
local oscillator frequency $\omega_L$ in each direction. For signal of limited
spectral distribution, this transformation is visualized in
Figure~\ref{fig:contmixspect}.
\begin{figure}
  \centerline{\includegraphics{timediscrete/mixing1.pdf}}
  \caption{\label{fig:contmixspect} Transformation of a spectrum
    $\tilde{f}(\omega)$ through multiplication of $f(t)$ with a local
    oscillator $m(t)$ at frequency $\omega_L$.}
\end{figure}

For a time-discrete signal $f_k$, the mechanism works in exactly the same
way, as the derivation of the spectrum $\tilde{g}(\omega)$ just uses the
convolution theorem. However, the spectrum of a time-discrete signal contains
periodic copies of the spectral component of interest, with a period of the
sampling frequency $\omega_S$. The spectrum of the product and the local
oscillator contains more components that need to be taken care of.

For the simplest case, consider a local oscillator frequency that is a quarter
of the sampling frequency, $\omega_L=\omega_S/4$. This has the advantage that
the time-domain signal of the local oscillator has a very simple form. Possible
choices (shown in figure~\ref{fig:localosc1}) are:
\begin{figure}
  \centerline{\includegraphics{timediscrete/localosc1.pdf}}
  \caption{\label{fig:localosc1} Sampled local oscillators $m(t)$ and $m'(t)$
    at frequency $\omega_L=\omega_S/4$.}
\end{figure}
%
\begin{equation}
  m_k=\begin{cases}+1\,&\text{if}\quad k\!\mod4 =0\\
        -1\,&\text{if}\quad k\!\mod4 =2\\
        0&\text{otherwise}\end{cases}
      \qquad\text{or}\qquad
  m'_k=\begin{cases}+1\,&\text{if}\quad k\!\mod4 =0,1\\
        -1\,&\text{if}\quad k\!\mod4 =2,3\end{cases}
\end{equation}
This first signal corresponds to a sampled version of a cosine function,
\begin{equation}\label{eq:localosc}
  m(t)=\cos(\omega_L t)\cdot\comb_{\Delta
    t}(t)=\sum\limits_{k=-\infty}^{\infty}m_k\delta(t-k\cdot\Delta t)\,,
\end{equation}
the second one to a sampled cosine function with a phase shift:
\begin{equation}
  m'(t)=\sqrt{2}\cos(\omega_L t-\frac{\pi}{4})\cdot\comb_{\Delta
    t}(t)=\sum\limits_{k=-\infty}^{\infty}m'_k\delta(t-k\cdot\Delta t)
\end{equation}
The corresponding spectrum of the local oscillator can directly be obtained
from (\ref{eq:localosc}):
%{\bf (more 2pi doubts)}
\begin{eqnarray}
  \tilde{m}(\omega)&=&{\cal F}\left[\cos(\omega_L t)\cdot\comb_{\Delta
                       t}(t)\right]\nonumber\\
                   &=&\frac{1}{2\pi}\left(\pi(\delta(\omega-\omega_L) +
                       \delta(\omega+\omega_L))\right) \ast
                       \left(\omega_S\sum\limits_{k=-\infty}^{\infty}
                       \delta(\omega-k\omega_S)\right) \nonumber\\
                   &=&\omega_S\sum\limits_{k=-\infty}^{\infty}\frac{1}{2}\left(
                       \delta(\omega-k\omega_S+\omega_L)+
                       \delta(\omega-k\omega_S-\omega_L)\right)\,,
\end{eqnarray}
which is a Dirac comb in frequency with peaks at odd multiples of $\omega_L$.
%
The expression for the resulting spectrum is again obtained by
convolution of the local oscillator spectrum with the signal spectrum:
%{\bf (more 2pi doubts...)}
\begin{eqnarray}
  \tilde{g}(\omega)&=&\tilde{f}(\omega)\ast \tilde{m}(\omega)
                   =\int\limits_{-\infty}^{\infty}\!\tilde{f}(\omega-\omega')\cdot
                       \tilde{m}(\omega')d\omega' \nonumber\\
                   &=&\int\limits_{-\infty}^{\infty}\!\tilde{f}(\omega-\omega')\cdot
                       \left[\omega_S\sum\limits_{k=-\infty}^{\infty}\frac{1}{2}\left(
                       \delta(\omega'-k\omega_S+\omega_L)+
                       \delta(\omega'-k\omega_S-\omega_L)
                       \right)\right]d\omega' \nonumber\\
                   &=&\omega_S\sum\limits_{k=-\infty}^{\infty}\frac{1}{2}\left(
                       \tilde{f}(\omega+\omega_L-k\omega_S)+
                       \tilde{f}(\omega-\omega_L-k\omega_S)\right)
\end{eqnarray}

The multiplication of the time-discrete signal $f_k$ with the local oscillator
$m_k$ then generates copies of the original spectrum, displaced by
$\pm\omega_L$. 
%
Assuming that the original signal has its spectral content in a
small band around the local oscillator frequency, that band will be moved to
frequency 0 (baseband), and another copy of it around $2\omega_L=\omega_S/2$,
just at the border between Nyquist  zones 1 and 2 (see figure
\ref{fig:discmixspect}).
%
\begin{figure}
  \centerline{\includegraphics{timediscrete/mixing2.pdf}}
  \caption{\label{fig:discmixspect} Transformation of a spectrum
    $\tilde{f}(\omega)$ by multiplying with a time-discrete local oscillator
    at a quarter of the sampling frequency $\omega_L$.}
\end{figure}
The spectrum $\tilde{g}(\omega)$ is again a series of small spectral bands as
$\tilde{f}(\omega)$, but this time appearing at even multiples of $\omega_L$.
This looks like a spectrum that could be sampled at half the sampling
frequency, without loosing any information contained in the product spectrum,
as long as the interesting bandwidth in the original spectrum is narrower than
$\omega_L$. Such a process is called {\em decimation}. If the original
sampling rate is to be retained, the contributions near $\omega_S/2$ have to
be suppressed by proper filtering.

The original signal spectrum in a practical system
$\tilde{f}(\omega)$ will also contain components near $\omega=0$, be it a
simple constant offset, or very slowly varying contributions to the sampled
signal. These undesired spectral components will be moved to regions around
$\omega_L+k\omega_S$, and $-\omega_L+k\omega_S$ for all integer
$k\in\mathbb{Z}$, populating regions at odd multiples of
$\omega_L$. Therefore, appropriate filtering techniques have to be applied in
order to suppress these contributions.

So far, a mixing scenario was considered with a local oscillator frequency
$\omega_L=\omega_S/4$. This results in the largest bandwidth that can be
transferred to the baseband, and is particularly easy to implement
numerically, as the multiplication process is restricted to values of
$m_k=\pm1$, which just requires a sign change. However, if other
frequencies are desired, other sinusoidal local oscillator values $m_k$
can be used in a similar way. Alternatively, other local oscillator waveforms,
e.g. rectangular signals with different frequencies  with values $m_k\pm1$ can
be used. Such approaches can be useful to implement lock-in amplifiers or
receiver designs where the mixing process is carried out numerically at a
fixed sampling rate in a signal processor for different target
frequencies. The derivation of the resulting spectral shift can be carried out
in the same way, with a local oscillator spectrum $\tilde{m}(\omega)$ that may
contain higher harmonics of $\omega_L$. 

\section{Decimation of time-discrete signals}
If the information of a time-discrete signal $f_k$ is contained in a
sufficiently small spectral band in $\tilde{f}(\omega)$, there is a redundancy
in the samples. This suggests that only a subset of the samples need to be
taken. In the following, we consider to a situation where the sampling rate is
reduced by a factor of $M$ by keeping only every $M$-th sample, and skipping
the $M-1$ samples in between. The decimated time-discrete signal $f'_k$ will
have a new sampling frequency $\omega'_S=\omega_S/M$, resulting in a more
densely spaced spectrum, with a lower periodicity:
\begin{equation}
  \tilde{f}'(\omega)=\tilde{f}'(\omega+k\cdot\omega_S/M)\,,\quad k\in\mathbb{Z}
\end{equation}
The sampling can be understood by introducing a decimating signal
$d_k$:
\begin{equation}
  f'_k = f_k\cdot d_k\,,\quad\text{with}\quad
  d_k=\begin{cases}1&\text{for}\quad k\mod M=0\,,\\
        0&\text{otherwise}\,.\end{cases}
\end{equation}
This decimating signal $d_k$ can be expressed through its discrete Fourier
transformation,
\begin{equation}
  d_k=\frac{1}{M}\sum\limits_{l=0}^Me^{2\pi i\frac{l}{M}}\,,
\end{equation}
which can be used to evaluate the spectrum of the decimated signal:
\begin{eqnarray}
  \tilde{f}'_R(\omega)&=&\sum\limits_{k=-\infty}^{\infty}
                          f'_ke^{-i\omega\Delta t k} =
                          \sum\limits_{k=-\infty}^{\infty}
                          f_kd_ke^{-i\omega\Delta  t k}\nonumber\\
                      &=&\sum\limits_{k=-\infty}^{\infty}f_k
                          \frac{1}{M}\sum\limits_{l=0}^Me^{2\pi i\frac{l}{M}}
                          e^{-i\omega\Delta t k}\nonumber\\
                      &=&\frac{1}{M}\sum\limits_{l=0}^M
                          \sum\limits_{k=-\infty}^{\infty}f_k
                          e^{-i(\omega-\frac{\omega_S}{M}l)\Delta tk}\nonumber\\
                      &=&\frac{1}{M}\sum\limits_{l=0}^M
                          \tilde{f}_R(\omega-\frac{\omega_S}{M}l)
\end{eqnarray}
The reconstructed spectrum of the decimated signal $f'_k$ is $M$ copies of the
spectrum $\tilde{f}_R(\omega)$ of the original signal, displaced by the
multiples of the new sampling frequency $\omega_S/M$.
As the size of the Nyquist zones are reduced, an issue arises when some
spectral content of the original stream overlaps with the new Nyquist
zones corresponding to the decimated sampling frequency, leading to undesired
artifacts in the decimated stream. Therefore, the spectral content of the
original signal in these areas needs to be reduced with an adequate low pass
filter to acceptable levels, in the same way as sampling time-continuous
signals in the first place.
\begin{figure}
  \centerline{\includegraphics{timediscrete/downsampling.pdf}}
\caption{\label{fig:downsampling}Decimating a time-discrete signal generates a
set of copies of the original spectrum $\tilde{f}(\omega)$. Shown here the
spectral change for decimating by $M=3$, i.e., $\omega'_S=\omega_S/3$.}
\end{figure}
\section{Noise sources in time-discrete systems}
Noise on a signal is refered to a part of the signal that is not reproducible
(or correlated with any useful signal part). A noisy signal is treated as
the sum of a clean signal $f(t)$, and a stochastical noise function $n(t)$:
\begin{equation}
  f_N(t)=f(t)+n(t)
\end{equation}
The stochastic nature of the noise means requires that an ensemble average has
to be carried out for the noise part. Noise is then characterized by
statistical properties, such as a standard deviation, and a mean. Typically,
\begin{equation}
  \langle n(t)\rangle=0\,,\quad\langle n^2(t)\rangle=:\sigma_n^2\neq0\,.
\end{equation}
Here, $\langle\cdot\rangle$ indicates an ensemble average over different noise
possibilities, and $\sigma_n$ is the standard deviation of the noise
amplitude.

For a mathematical description, consider first the definition of a power $P$
contained in a continuous-time signal $f(t)$, or discrete time sampled $f_k$,
\begin{equation}
  P=\frac{1}{T}\int\limits_{0}^{T}f^2(t)dt\quad\text{and}\quad
  P=\frac{1}{N}\sum\limits_{k=0}^{N-1}f_k^2
\end{equation}
for some averaging time intervall $T$ or sample number $N$. This corresponds
to the ``root mean square'' of the signal function $f(t)$ or time-discrete set
$f_k$. For a noise signal $n(t)$, this quantity can also be evaluated:
\begin{equation}
  P=\frac{1}{T}\int\limits_{0}^{T}\left\langle n^2(t)\right\rangle dt = \sigma_n^2
\end{equation}
  
% As noise may have a particular structure in time or frequency, the
% quantity of an autocorrelation $C(t)$ function for a function $f(t)$ is
% introduced,
% \begin{equation}
%   C(t)=\int\limits_{-\infty}^{\infty}f^*(t')f(t'+t)dt'\,,
% \end{equation}
% where $f*$ denotes the complex conjugate. The Wiener-Kinchine theorem now
% states that the Fourier transformation of $C(t)$ is given by
% \begin{equation}
%   {\cal F}[C(t)](\omega) = \int\limits_{-\infty}^{\infty}C(t)e^{-i\omega t}dt=
%   \left|\tilde{f}^2(\omega)\right|\,,
% \end{equation}
% The autocorrelation function also exists for a noise function $n(t)$,
% \begin{equation}
%   C(t)=\int\limits_{-\infty}^{\infty}\left\langle n^*(t')n(t'+t)\right\rangle
%   dt'\,,
% \end{equation}
% and helps to characterize noise in systems with a particular spectral
% transmission.




\section{Noise propagation in time-discrete systems}

% 

%%% Local Variables:
%%% mode: latex
%%% TeX-master: "../eclass"
%%% End:
