\documentclass[conference]{IEEEtran}
\usepackage{amsfonts,amsmath}
\usepackage[dvips]{epsfig}
\DeclareGraphicsRule{.eps.gz}{eps}{.eps.bb}{`gunzip -c #1}
\DeclareGraphicsExtensions{.eps}
\graphicspath{{figures/}{eps/}{pdf/}{png/}{images/eps/}{images/pdf/}{images/png/}}

%% macros
\providecommand{\code}[1]{\textsf{#1}}
\providecommand{\up}[1]{\textsuperscript{#1}}

\begin{document}

\title{A fast OFDM frequency domain channel simulator for Matlab/Octave}
\author{%
  \IEEEauthorblockN{Marc de Courville, Laurent Mazet and V\'eronique Buzenac}
  \IEEEauthorblockA{Motorola-Labs Paris -- Parc Les  Algorithmes, Saint Aubin -- 91193 Gif sur Yvette Cedex -- France \\
    $\{$marc.de.courville, laurent.mazet, veronique.buzenac$\}$@motorola.com}
}


\maketitle

\begin{abstract}
\end{abstract}

\begin{IEEEkeywords}
Channel model, OFDM, Fast implementation
\end{IEEEkeywords}

\section{Introduction}
\label{sec:introduction}

The purpose of this contribution is to present a fast channel simulator
dedicated to OFDM systems.  This specific channel simulator has the
advantage to work directly in the frequency domain through parallel
scalar multiplications which enables a very low complexity
implementation. It implements 3 kind of Doppler: Jakes, Rician and
gaussian.  However this simulator has some restrictions:
\begin{itemize}
\item it assumes that the overall channel bandwidth is shorter than
  the cyclic prefix duration
\item the channel is assumed static during one OFDM baud (or block)
\item for genericity the convolution of the RX and TX filters is assumed
to be a rectangular window of duration the sampling period. In other words the
current channel simulator doesn't include TX and RX filter.
\end{itemize}
The channel simulator has been implemented in Matlab/Octave and the
usage is described in this document.

Moreover the simulator embeds natively all the ETSI BRAN 5.2GHz
channel models as defined in \cite{misc.bran98} originally defined
over a bandwidth of 100MHz and is capable of deriving the equivalent
tap model channel viewed at the receiver after sampling in order to be
able to process variable bandwidth signals.
% This makes it attractive for assessing the 4.9GHz MCLB CGISS system which due to the scarcity of the spectrum (only 50MHz allocated by FCC) operates over 20/10/5MHz channels.

Section \ref{sec:brief-how-to} describes briefly how to interpolate a 
20/10/5MHz bandwidth channel from a 100MHz bandwidth channel model
through basic windowed down-sampling.
Section \ref{sec:doppler-simulation} describes the implementation of the
Doppler spectrum.
Finally, section \ref{sec:how-use-matlaboctave} describes the channel
simulator code and how to interface it..

\section{Brief how-to for changing sampling resolution}
\label{sec:brief-how-to}
%\section{Channel projection-denoising and change of sampling resolution}
%\label{sec:chann-proj-deno}

\subsection{Introduction}

The goal of this section is to explain how the change in sampling is
implemented in the simulator going back to the filter-bank theory.
The ultimate objective being to derive from the 100MHz channel specification
the channel viewed at the receiver after sampling at 20/10/5MHz.

\subsection{Notations}
Subscripts $_t$ and $_f$ denotes respectively the time or frequency domain
coefficients. Integer upper-script represents the size of the vector considered.
$\mathbf{C}$ is the channel vector. $N$ is the total number of carriers over 
20/10/5MHz bandwidth, while $D$ denotes the number of samples contained in the
Cyclic Prefix and $K$ is equal to the number of carriers bearing information
(non set to zero).

\subsection{Modeling of the channel over a smaller bandwidth}

The problem addressed is the derivation of the discrete time domain
equivalent filter expression capturing both the channel filter $C(z)$
defined over a larger bandwidth $N\times B$ than the real one $B$ of the
digital system and also the TX and RX filters: $F_{\text RX}(z)$ and
$F_{\text TX}(z)$.  In the present case of interest $B=20$MHz and
$N=5$ implying that $C(z)$ has been measured over a 100MHz
bandwidth which is the case for the channels specified by ETSI BRAN
for 5GHz OFDM WLANs.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=.8\linewidth]{filtering}
  \caption{Channel filtering modeling with TX and RX filters}
  \label{fig:channel}
\end{figure}

The filtering process (c.f. figure \ref{fig:channel}) is summarized 
by the following set of equations:
\begin{align}
H(z) & := F_{\text RX}(z) H(z) F_{\text TX}(z)\\
Y(z) & =H(z)X(z^N)
\end{align}
where $X(z^N)$ represent the output of the DAC result of the $N$ fold expansion
of the signal $X(z)$ over the full $NB$ bandwidth.

In order to specify the expression of the output of the decimation, it is suitable
to recall the polyphase decomposition formalism \cite{book.vaidyanathan93} illustrated by figure
\ref{fig:polyphase} applied on $Y$ for order $N$:
\begin{align}
Y(z) &= \sum_{k=0}^{N-1}Y_k(z^N)z^{-k}\\
Y_0(z^N) & = \frac{1}{N}\sum_{k=0}^{N-1}Y\left(ze^{j\frac{2\pi}{N}k}\right) \\
Y_i(z^N) & = \lfloor Y(z)z^i \rfloor_0
\end{align}
$Y_i(z)$ is the $i$th polyphase component of $Y(z)$ gathering all its coefficient aligned with
index $i$ modulo $N$.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=.8\linewidth]{polyphaseI}
  \caption{Polyphase representation illustration}
  \label{fig:polyphase}
\end{figure}

Decimation by a factor $N$ of $Y(z)$ amounts to isolate the $0$th polyphase 
component of $Y(z)$: $Y_0(z)$ thus,
\begin{align}
Y_0(z^N) & = \frac{1}{N}\sum_{k=0}^{N-1}Y\left(ze^{j\frac{2\pi}{N}k}\right)\\
Y_0(z^N) & = \frac{1}{N}X(z^N)\sum_{k=0}^{N-1}H\left(ze^{j\frac{2\pi}{N}k}\right)\\
Y_0(z) & = X(z)H_0(z)
\end{align}
This result shows that in the digital domain over $B$ the filtering can be
represented by the direct filtering of $X(z)$ by the $0$th polyphase component
of the aggregate filter $H_0(z)$ (the $N$ fold version of $H(z)$).

\subsubsection{Rectangular windowing strategy for TX and RX filters}

%\begin{itemize}

Let assume that the overall $H(z)$ filter is a rectangular window:
\begin{equation}
F_{\text RX}(z)F_{\text TX}(z)=\sum_{k=0}^{N-1}z^k
\end{equation}

Note that as justified below when using this window the equivalent filter
$H_0(z)$ is equal to a basic shift and sum of the coefficients of $C(z)$
realigned on the time sampling grid corresponding to bandwidth $B$.
\begin{align}
C(z)&=\sum{i=0}^{N-1}C_i(z^N)z^{-i}\\
H(z)&=F_{\text RX}(z)C(z)F_{\text TX}(z)\\
&=\sum_{i=0}^{N-1}\sum{k=0}^{N-1}z^{-i+k}C_i(z^N)\\
H_0(z^N)&=\sum_{l=0}^{N-1}C_l(z^N)
\end{align}

%\item Hanning: minimize global $H(z)$ filter size
%\item Hamming: minimize global $H(z)$ filter size
%\item Prototype IEEE802.11a demo filter: minimize global $H(z)$ filter size
%\end{itemize}

Since $H_0(z)$ is simply the equivalent channel seen in the discrete domain after
sampling. Thus the channel simulator directly uses this trick to derive
from the time domain $H(z)$ over 100MHz the $H_0(z)$ after sub-sampling

In order to simulate TX and RX filters, the implementation we propose is:
\begin{enumerate}
\item compute the adequate IIR filters (e.g. Chebychev) that matche the analog desired responses
\item find good causal FIR approximations of these IIR filters: $F_{\text RX}(z)$ and $F_{\text TX}(z)$
\item at each update of the channel model $C(z)$:
  \begin{enumerate}
  \item compute the polynomial product (convolution) of $H(z)=F_{\text RX}(z)C(z)F_{\text TX}(z)$
  \item extract the $0$th polyphase of $H(z)$
  \item use $H_0(z)$ to filter the block
  \end{enumerate}
\end{enumerate}

\section{Some background information about projection denoising algorithm and links with change in
bandwidth resolution}

This material should be considered as background information but complements well
the previous section. Indeed it intends to describe how to operate a change in
the bandwidth resolution (i.e. the number of carriers in a given bandwidth)
for deriving the expression of the channel.

The channel projection denoising algorithm \cite{report.sandell96} aims
at exploiting the prior knowledge that the channel impulse response in
the time domain is by construction of the system shorter than the guard
interval duration.  It enables to remove significantly some estimation
noise when deriving the channel from the preambles and yields
to an enhanced estimator.

\subsection{The channel projection denoising algorithm}

When assuming that the time domain channel $\mathbf{C}_t^N$ is confined in the
CP duration,
the frequency domain channel that can be measured (i.e.: the $K$ carrier to
which we have access) $\mathbf{C}_f^K$ is expressed as:
\begin{align}
\mathbf{C}_f^K & =\mathbf{S}_f^K\mathbf{F}_N\mathbf{S}_t^D\mathbf{C}_t^N\\
\mathbf{M} & :=\mathbf{S}_f^K\mathbf{F}_N\mathbf{S}_t^D\\
\widehat{\mathbf{C}}_t^N & =\mathbf{M}^{\dag}\mathbf{C}_f^K
\end{align}
where $\widehat{\mathbf{C}}_t^N$ is the denoised channel in time for which the estimation noise
has been removed by imposing the time domain mask of the guard time on the
overall impulse response.
$\mathbf{C}S_t^K$ and $\mathbf{C}S_f^K$ denote respectively the time and frequency selecting 
matrices zeroing either the $N-D$ coefficients in time out of the guard time
duration or the $N-K$ carriers in frequency that should be null.
\begin{figure}[htbp]
  \centering
  \includegraphics[width=.8\linewidth]{fourier-truc}
  \caption{Truncation of the Fourier matrix by $S_f^K$ and $S_t^D$}
  \label{fig:fourier-truc}
\end{figure}
Since $\mathbf{M}$ is of full column rank: $\mathbf{M}^{\dag}=(\mathbf{M}^H\mathbf{M})^{-1}\mathbf{M}^H$:
\begin{align}
\widehat{\mathbf{C}}_t^N&=(\mathbf{M}^H\mathbf{M})^{-1}\mathbf{M}^H\mathbf{C}_f^K\\
\widehat{\mathbf{C}}_f^K&=\mathbf{S}_f^K\mathbf{F}_N(\mathbf{M}^H\mathbf{M})^{-1}\mathbf{M}^H\mathbf{C}_f^K
\end{align}

\subsection{Link with change of time/frequency resolutions grids}

Known result: when changing the sampling frequency the channel 
frequency response doesn't change. However the digital time domain one does.
What is the virtual time domain channel profile that would lead to the
same frequency domain response when using a lower sampling rate?
Be $L$ the number of taps that corresponds to the change in sampling
rate: i.e. in one block of same duration instead of $N$ samples we 
observe now $L<N$ ones. We thus have:
\begin{equation}
\widehat{\mathbf{C}}_t^L=\mathbf{F}^{-1}_L\mathbf{S}_f^N\mathbf{F}_N\mathbf{S}_t^D\mathbf{C}_t^N
\end{equation}

\section{Doppler simulation}
\label{sec:doppler-simulation}

The Doppler is applied on each of the time domain taps. The taps
are updated at each OFDM baud. The Doppler function is a classical
implementation based on an oversampling of the Doppler spectrum
for derivating the evolution of the channel taps through a cubic
interpolation.  The code supports gaussian, Jakes and Rician Doppler
models.

\section{How to use the Matlab/Octave channel simulator}
\label{sec:how-use-matlaboctave}

\subsection{General use: the \code{ofdmchannel} function}
\label{sec:general-use}
This function performs the filtering in the frequency domain
of several OFDM symbol as indicated its `help' message:
\begin{tiny}\begin{verbatim}
usage: [y, cir, cfr] = ofdmchannel (x, model, sig, [dspeed])

Filters several OFDM symbol by the channel named model.

x: matrix stacking in columns the frequency domain OFDM symbol vectors.
model: channel model. Admitted values are 'Awgn', 'Rayleigh', 'Bran-x'
       (for more information, see cirpowers)
sig: noise standard deviation.
dspeed: terminal velocity in m/s [default 0].

cir: matrix gathering the variations of the time domain channel impulse
     response tap vector with a granularity of the OFDM baud according
     to the Doppler (Channel Impulse Response).
cfr: corresponding frequency domain channel matrix (Channel Frequency
     Response).
y: output signal matrix resulting of the convolution of each x column
   by each cir column.

See also: genfading, cirpowers, doppler
\end{verbatim}\end{tiny}

The \code{x} parameter is a matrix of data to be transmitted into the
channel (each column corresponds to an OFDM symbol), the \code{model}
parameter defines the channel model, the \code{sig} parameter is the
noise standard deviation $\sigma$ defined such that
$\sigma^2 =10^{-\frac{\mathrm{SNR}}{10}}$ and the last and optional parameter,
denoted by \code{dspeed} is the terminal velocity. By default,
\code{dspeed} is set to 0.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=.8\linewidth]{filter}
  \caption{Channel filtering model}
  \label{fig:channel-filter}
\end{figure}

This function returns \code{y}, a matrix containing the
result of the filtering of each column of \code{x} by each
column of \code{cir} polluted by noise.
\code{cir} is a matrix capturing the channel impulse responses 
variations along time according the the Doppler spectrum at each
OFDM baud.

The \code{model} parameter, defining the channel model, follows this
set of rules:
\begin{itemize}
\item \code{'Agwn'} means that no filtering is applied: only
  a gaussian noise of variance $\sigma^2$ is added.
\item \code{'Rayleigh'} defines a channel with Rayleigh fading on each tap.
  Using a string such as \code{'Rayleigh($n$)'} limits the number of
  fading to $n$. One can also set the relative power of the channel
  and also the oversampling factor using a string as
  \code{'Rayleigh($n$)\_$p$\%$o$'} ($n$ is the number of fading, $p$ is
  the relative power, and $o$ is the oversampling factor). For
  example, a string such as \code{'Raleigh(16)\_-2\%5'} specifies 16
  taps of variance -2dB oversampled by a factor 5.
\item \code{'Bran-$x$'} defines the Bran $x$ channel model (where $x$
  belongs to the set $\{$\code{A}, \code{B}, \code{C}, \code{E}$\}$).
  Using a string such as \code{'Bran-E(5)'} specifies a Bran E model over a
  5MHz bandwidth. Possible bandwidths are 5MHz, 10MHz and 20MHz.
  Default bandwidth for Bran models is 20MHz.
\end{itemize}


\subsection{The fading generator: \code{genfading}}
\label{sec:fading-gener-codeg}

This function is internally called by \code{ofdmchannel}. This paragraph
is intended for developers only. The function `help' message is:
\begin{tiny}\begin{verbatim}
usage: [cir, oversamp] = genfadings (len, model, dspeed, [nbtaps, [oversamp]])

Generate channel impulse response

len: variation length
model: channel model. Admitted values are 'Awgn', 'Rayleigh', 'Bran-x' (where
       x is either A, B, C or E). For a Rayleigh channel, a string such as 
       'Raleigh(16)_-2%5' specifies 16 taps of variance -2dB over
       sampled by a factor 5.  For Bran channel, a string such as
       'Bran-E(5)'specifies a Bran E model over a band of 5 MHz.
       Possible bandwidths are 5 MHz, 10 MHz and 20 MHz. Default
       bandwidth for Bran models is 20 MHz.
dspeed: Doppler speed
nbtaps: number of taps (useful for awgn channel only) [64]
oversamp: over sampling factor (useful for awgn channel only) [1]

cir: channel impulse response (nbtaps x len)
oversamp: oversampling factor.

See also: cirpower, doppler
\end{verbatim}\end{tiny}

\subsection{The channel impulse response profile: \code{cirpower}}
\label{sec:chann-impulse-resp}

This function is internally called by \code{genfading}. This paragraph
is intended for developers only. The function `help' message is:
\begin{tiny}\begin{verbatim}
usage: [tapvars, oversamp] = cirpowers (model, [nbtaps, [oversamp]])

This function returns the taps variances of a static multipath channel
and the oversampling factor.

model: channel model. Admitted values are 'Awgn', 'Rayleigh', 'Bran-x' (where
       x is A, B, C or E). For Rayleigh channel, a string such as 
       'Raleigh(16)_-2%5' specifies 16 taps of variance -2dB over
       sampled by a factor 5.  For Bran channel, a string such as
       'Bran-E(5)'specifies a Bran E model over a band of 5 MHz.
       Possible bandwidths are 5 MHz, 10 MHz and 20 MHz. Default
       bandwidth for Bran models is 20 MHz.
nbtaps: for Rayleigh channel, you can specify the number of taps
        (by default 64)
oversamp: also for Rayleigh channel, you can specify an oversampling
          factor (by default 1).

tapvars: taps variances.
oversamp: oversampling factor.
doptype: Doppler type.
\end{verbatim}\end{tiny}

\subsection{The Doppler effect generator: \code{code}}
\label{sec:doppl-effect-gener}

This function is internally called by \code{genfading}. This paragraph
is intended for developers only. The function `help' message is:
\begin{tiny}\begin{verbatim}
usage: [dop] = doppler(len, fd, fs, model)

Generate normalize Doppler variations

len: variation length or [number taps, variation length]
fd: maximum Doppler frequency
fs: update signal frequency
model: Doppler model ('gauss', 'jakes', 'none', 'rice') ['jakes']

dopvar: normalize Doppler variation
\end{verbatim}\end{tiny}

\section{Conclusion}


This simulator is intended for quick assessment of the 4.9GHz
system performance through direct frequency domain processing.

\bibliographystyle{IEEEtran}
\bibliography{abbrev,alien,standard}
\end{document}

