\chapter{Frequency dependent transfer functions} \label{fd_transfer_functions}

\section{Introduction}

Frequency dependent transfer functions are implemented in the Spice models using s-domain transfer functions, this includes the propagation corrections and transfer impedance terms. The coefficients of the s-domain transfer functions are required to provide the 'best fit' to a frequency dependent transfer function evaluated at a discrete set of frequencies. The software includes a process whereby the coefficients of an s-domain transfer function (taking the form of a rational function in s), which provides an approximation to a function, are determined. 

\section{Filter fitting process for s-domain transfer functions} \label{filter_fit}

The process implemented is based on that described in \cite{Dawson}.

The starting point for the process is a set of complex transfer function values $H'\left(j\omega_i\right)$ at a set of discrete frequencies $i=1..n$ and the aim is to find a rational function $H(j\omega)$ which best fits the available data points.

The transfer function takes the form:
%
\begin{equation} \label{eq:H_rational_function}
H\left(j\omega\right)=\frac{a_{0}+a_{1}\left(\frac{j\omega}{\omega_{0}}\right)+a_{2}\left(\frac{j\omega}{\omega_{0}}\right)^{2}+\dots+
a_{nn}\left(\frac{j\omega}{\omega_{0}}\right)^{nn}}
{1+b_{1}\left(\frac{j\omega}{\omega_{0}}\right)+b_{2}\left(\frac{j\omega}{\omega_{0}}\right)^{2}+\dots+
b_{nd}\left(\frac{j\omega}{\omega_{0}}\right)^{nd}}
=\frac{A(j\omega)}{1+B(j\omega)}
\end{equation}
%
where the order of the numerator and the denominator are to be specified and the coefficients $a_j$ and $b_j$ are to be determined.
If we assume that the numerator and denominator order are specified as $nn$ and $nd$ then we have unknown coefficients $a_j, j=0 .. nn$ and $b_j, j=1 .. nd$.

At a frequency $w_i$ the error between the filter fit model and the transfer function value is
\begin{equation} \label{eq:err1}
err(j\omega_i)=frac{A(j\omega_i)}{1+B(j\omega_i)}-H'\left(j\omega_i\right)
\end{equation}

multiplying equation \ref{eq:err1} by $1+B(j\omega_i)$ gives

\begin{equation} \label{eq:err2}
\left(1+B(j\omega_i)\right)err(j\omega_i)=err'(j\omega_i)=A(j\omega_i)-\left(1+B(j\omega_i)\right)H'\left(j\omega_i\right)
\end{equation}

The left hand side can be thought of as a weighted error function, $err'$.

Evaluating equation \ref{eq:err2} for each frequency $\omega_i$ gives a set of n equations 

\begin{multline} \label{eq:err3}
\begin{pmatrix} 
err'(j\omega_1) \\ err'(j\omega_2) \\ \vdots \\ err'(j\omega_n)
\end{pmatrix}
=
\begin{pmatrix} 
-H'(j\omega_1) \\ -H'(j\omega_2) \\ \vdots \\ -H'(j\omega_n)
\end{pmatrix}
+ \\[10pt]
\begin{bmatrix}
-\left(j\omega_1 \right) H'(j\omega_1)  & \hdots & -\left(j\omega_1 \right)^{nd} H'(j\omega_1) 
& 1 & \left(j\omega_1 \right)  & \hdots & \left(j\omega_1 \right)^{nn} \\ 
-\left(j\omega_2 \right) H'(j\omega_2)  & \hdots & -\left(j\omega_2 \right)^{nd} H'(j\omega_2) 
& 1 & \left(j\omega_2 \right)  & \hdots & \left(j\omega_2 \right)^{nn} \\ 
\vdots & \vdots & \vdots  & \vdots & \vdots & \vdots & \vdots & \vdots   \\ 
-\left(j\omega_n \right) H'(j\omega_n) &  \hdots & -\left(j\omega_n \right)^{nd} H'(j\omega_n) 
& 1 & \left(j\omega_n \right)  & \hdots & \left(j\omega_n \right)^{nn} 
\end{bmatrix}
\begin{pmatrix}  
b_1 \\  \vdots \\ b_{nd} \\ a_0 \\ a_1 \\ \vdots \\ a_{nn}
\end{pmatrix}
\end{multline}

The coefficients of the transfer function are required to be real. This can be enforced by augmenting the set of equations with an equal number of equations evaluated at negative frequency where we impose the condition $H(-j\omega)=H^{*}(j\omega)$.

\begin{multline} \label{eq:err4}
\begin{pmatrix}
err'(-j\omega_1) \\ err'(-j\omega_2) \\ \vdots \\ err'(-j\omega_n)
\end{pmatrix}
=
\begin{pmatrix}
-H'^{*}(j\omega_1) \\ -H'^{*}(j\omega_2) \\ \vdots \\ -H'^{*}(j\omega_n)
\end{pmatrix}
+ \\[10pt]
\begin{bmatrix}
-\left(-j\omega_1 \right) H'^{*}(j\omega_1)  & \hdots & -\left(-j\omega_1 \right)^{nd} H'^{*}(j\omega_1) 
& 1 & \left(-j\omega_1 \right)   & \hdots & \left(-j\omega_1 \right)^{nn} \\ 
-\left(-j\omega_2 \right) H'^{*}(j\omega_2)  & \hdots & -\left(-j\omega_2 \right)^{nd} H'^{*}(j\omega_2) 
& 1 & \left(-j\omega_2 \right)   & \hdots & \left(-j\omega_2 \right)^{nn} \\ 
\vdots & \vdots & \vdots & \vdots & \vdots  & \vdots & \vdots  & \vdots  \\ 
-\left(-j\omega_n \right) H'^{*}(j\omega_n)  & \hdots & -\left(-j\omega_n \right)^{nd} H'^{*}(j\omega_n) 
& 1 & \left(-j\omega_n \right)   & \hdots & \left(-j\omega_n \right)^{nn} 
\end{bmatrix}
\begin{pmatrix}
b_1 \\  \vdots \\ b_{nd} \\ a_0 \\ a_1 \\ \vdots \\ a_{nn}
\end{pmatrix}
\end{multline}

This combined set of equations is of the form

\begin{equation} \label{eq:err5}
\left( err' \right) \left(D\right) -\left[ X \right] \left( W \right)
\end{equation}

The mean square error is found as

\begin{equation} \label{eq:MSE_WH}
\begin{array}{c} 
\left\| err \right \| ^2 = \frac{1}{n} \Re \left\{ \left( err' \right)^{*} \left( err' \right) \right\} \\
=\frac{1}{n} \Re \left\{ \left(D\right)^{*} \left(D\right)+ \left(W\right)^{T} \left[X\right]^{*}\left[X\right] \left(W\right) -2\left(D\right) \left[X\right]^{*} \left(W\right) \right\}
\end{array}
\end{equation}

where $*$ denotes the conjugate transpose and $\Re$ denotes the real part 

If we define 

\begin{equation} 
\left[R\right]=\Re \left\{ \left[X\right]^{*}\left[X\right] \right\}
\end{equation}

and 

\begin{equation} 
\left(P\right)= \left[X\right]^{*}\left(D\right) 
\end{equation}

then

\begin{equation} \label{eq:MSE_WH2}
\left\| err \right \| ^2 = 
\frac{1}{n} \Re \left\{ \left(D\right)^{*} \left(D\right)+ \left(W\right)^{T} \left[R\right]\left(W\right) -2\left(P\right)^{T} \left(W\right) \right\}
\end{equation}

This expression is minimised with respect to the vector $\left( W \right)$ of unknown coefficients by differentiating with respect to $\left( W \right)$ and equating the result to zero i.e.

\begin{equation} \label{eq:MSE_WH3}
\frac{\partial\left\| err \right \| ^2 }{\partial \left(W \right)}= 
\frac{1}{n} \Re \left\{ 2 \left[R\right]\left(W\right) -2\left(P\right)^{T} \right\}=0
\end{equation}

This system of equations may be solved directly to give the unknown coefficient vector 

\begin{equation} \label{eq:MSE_WH4}
\left(W \right)=\left[R\right]^{-1}\left(P\right)
\end{equation}

An alternative approach to reaching equation \ref{eq:MSE_WH4} is to start from the set of equations \ref{eq:err3} and \ref{eq:err4}. If we set the error on the left hand side to zero then we have the overdetermined set of equations

\begin{equation} \label{eq:Moore_Penrose1}
\left(D \right)=\left[X\right]\left(W\right)
\end{equation}

A solution to this set of equations may be found by solving the 'normal equations' which are found by pre-multiplying equation \ref{eq:Moore_Penrose1} by $\left[X\right]^{*}$

\begin{equation} \label{eq:Moore_Penrose2}
\left[X\right]^{*}\left(D \right)=\left[X\right]^{*}\left[X\right]\left(W\right)
\end{equation}

The matrix $\left[X\right]^{*}\left[X\right]$ is a square matrix which may be inverted to give the solution

\begin{equation} \label{eq:Moore_Penrose_inverse}
\left(W\right)=\left[ \left[X\right]^{*}\left[X\right] \right]^{-1}\left[X\right]^{*}\left(D \right)
\end{equation}

The matrix $\left[ \left[X\right]^{*}\left[X\right] \right]^{-1} \left[X\right]^{*} $ is the Moore Penrose pseudo inverse and solves the original overdetermined system of equation in a least squares sense. 

\subsubsection{Stability of the transfer function}

The spice models are required to operate in the time domain which implies that all transfer functions in the model should be causal and that their impulse response should not diverge as $t \rightarrow \infty$. This implies that all poles of the transfer function $H(s)$ should lie in the left hand side of the s-plane. This condition is not naturally enforced by the method described above however experience has shown that if the function to be fitted is a physically reasonable then the filter fit process will produce a stable model.  


\subsubsection{Enforcing constraints on the transfer function}

It may be the case that the transfer function is subject to constraints. If this is the case then the process may be adapted such that the
constraints are naturally enforced. The filter fitting process is applied in a number of different aspects of the model. the filter fit required for the propagation correction function is subject to the constraint that $H\left(\omega=0\right)=1$. In this case
we calculate a model fit to $H'\left(\omega\right)=H\left(\omega\right)-1$ where $H'\left(\omega=0\right)=0$ thus $a'_0=0$ for the initial model. Once we have the coefficients for the initial model we can easily calculate the filter coefficients for $H\left(\omega\right)=1+H'\left(\omega\right)$.

\clearpage
