HarmonicBalanceText.tex

%
% Copyright © 2014 Peeter Joot.  All Rights Reserved.
% Licenced as described in the file LICENSE under the root directory of this GIT repository.
%
\section{Introduction.}
\IEEEPARstart{T}{he}
Harmonic Balance (HB) method is an approach to determine the steady state response of linear and nonlinear circuits to periodic sources.
It is assumed that all the outputs are also periodic, and that the frequency spectrum of all output signals will only contain frequencies that are integer multiples of a fundamental frequency (harmonics) common to all input signals.
It is also assumed that a finite number \( 2 N + 1 \) of these harmonics will be sufficient to accurately represent all the output signals.
As only the steady state response of the circuit is of interest, the state of the circuit only needs to be determined over one period.
The input sources are sampled at a set of \( 2 N + 1 \) discrete evenly spaced times, where these points in time are assumed to be close enough that the inputs at intermediate times are not significantly different from the sampled values.
Because the sources and responses are all assumed to be periodic, a Fourier series representation of the signals can be used.
The system equations are expressed in terms of the Fourier components, using the Discrete Fourier Transform (DFT).  It will be shown how the results from the Modified Nodal Analysis (MNA) procedure, an automated method of assembling the equations for a set of linear circuit elements, can be used to generate the HB system equations.  If the MNA procedure results in a set of \( R \) equations, then the HB equations will be a set of \( R \times (2 N + 1) \) equations.
For linear systems, HB can be interpreted as an independent circuit for each frequency, with nonlinearities in the system coupling these circuits.

The objective of HB is to represent the system in the form
\begin{equation}\label{eqn:ece1254projectReport:560}
\BI(\BV) = \BY \BV,
\end{equation}
where the (frequency domain) current \( \BI \) includes contributions associated with the linear sources (voltage and current sources), as well as any nonlinear contributions due to elements such as diodes.
The admittance matrix \( \BY \) has a block diagonal structure, with each block having a frequency dependency.
The vector \( \BV \) contains the Fourier coefficients of all the voltage and current unknowns in the system, including a term for each unknown at each frequency included in the Fourier series.
When the system includes nonlinear elements, iterative methods will be required to solve for \( \BV \).

\mysubsection{Modifications to the netlist syntax.}
\index{netlist syntax}
The Matlab based MNA implementation developed during the course assignments only had explicit support for DC sources, although it was possible to use the resulting MNA equations for non-constant sources with a suitable re-interpretation of the returned matrix equations.  For this HB work it was worthwhile to explicitly extend the netlist syntax to constant and oscillatory sources, allowing netlist lines of the form

\begin{center}
\textbf{Ilabel \(n_1\) \(n_2\) AC \(m\) \(\nu\) \([\phi]\)}
\end{center}
\begin{center}
\textbf{Vlabel \(n_1\) \(n_2\) AC \(m\) \(\nu\) \([\phi]\)}.
\end{center}
These specify current and voltage sources of the form \( m \cos\lr{ 2 \pi \nu - \phi } \), with currents directed from node \( n_1 \) to \( n_2 \), and voltages measured at \( n_1 \) relative to \( n_2 \).
The phase \( \phi \) defaults to zero if unspecified.

The netlist parser was also extended to allow the specification of diode lines of the form
\begin{center}
\textbf{Dlabel \(n_1\) \(n_2\) \(I_0\) \(V_{\textrm{T}}\)},
\end{center}

This represents a diode modeled as a current source \( i_{\textrm{d}}(t) = I_0 \lr{ e^{ \lr{v^{(n_1)}(t) - v^{(n_2)}(t)}/V_{\textrm{T}} } - 1 } \).  The constant portion of this model is added into the linear portion of the resulting matrix equations.  The nonlinear contributions of each diode line are returned from NodalAnalysis() as a list of structures with fields (io, type, vt, vp, vn), with type = 'exp'.

In addition to the diode model, two nonlinear power-law circuit elements were implemented.  The first is a nonlinear conductance \( g(t) = I_0 \lr{ \lr{v^{(n_1)}(t) - v^{(n_2)}(t)}/V_{\textrm{T}}}^k \) that is specified with lines of the form
\begin{center}
\textbf{Plabel \(n_1\) \(n_2\) \(I_0\) \(V_{\textrm{T}}\) \( k \)}.
\end{center}
Any number of these may be specified in parallel to form an arbitrary Taylor expansion.

The second nonlinear power-law circuit element that was implemented was a nonlinear voltage controlled voltage source
\( v_{n_1} - v_{n_2} = g \lr{ \lr{v_{n_1}^{(\textrm{c})}(t) - v_{n_2}^{(\textrm{c})}(t)}/V_{\textrm{T}}}^k \)
that is specified with lines of the form
\begin{center}
\textbf{Elabel \(n_1\) \(n_2\) \(n^{(\textrm{c})}_1\) \(n^{(\textrm{c})}_2\) \(g\) \(V_{\textrm{T}}\) \( k \)}.
\end{center}

Each of the power-law nonlinear elements is returned from NodalAnalysis() as a list of structures with fields (io, type, vt, vp, vn, exponent), with type = 'power'.
\section{Background.}
\mysubsection{Discrete Fourier Transform.}
\index{discrete Fourier transform}
An outline of DFT procedure is presented in \citep{giannini2004NonlinearMicrowaveCircuitDesign} \S A.4.  Due to some errors and inconsistencies in that text, it is worthwhile to restate the underlying ideas in their entirety.  This may not follow the convention that the authors of that text may have intended to use.

A periodic bandwidth limited signal with period \( T = 2 \pi/\omega_0 \), may be described by the Fourier series transform pair
%\begin{equation}\label{eqn:ece1254projectReport:60}
\boxedEquation{eqn:ece1254projectReport:60}{
\begin{aligned}
x(t_k) &= \sum_{n = -N}^N X_n e^{ j n \omega_0 t_k} \\
X_n &= \inv{2 N + 1} \sum_{k = -N}^N x(t_k) e^{-j n \omega_0 t_k},
\end{aligned}
}
%\end{equation}
where the signal is evaluated are evaluated only at the Nykvist times
\index{Nykvist time instant}
\( t_k \),
\begin{equation}\label{eqn:ece1254projectReport:20}
t_k = \frac{T k}{2 N + 1}, \qquad k \in [-N, \cdots N]
\end{equation}
all contained within the interior of the \( [-T/2, T/2] \) range of interest.
A proof of this inversion relationship is given in \cref{appendix:discreteFourierInversion}.

A matrix representation of the transform pair is desired.
Let \( x_k = x(t_k) \), and
substitute \( t_k \) from \cref{eqn:ece1254projectReport:20}
\begin{equation}\label{eqn:ece1254projectReport:140}
\begin{aligned}
x_k &= \sum_{n = -N}^N X_n e^{ 2 \pi j n k /(2 N + 1)} \\
X_n &= \inv{2 N + 1} \sum_{k = -N}^N x_k e^{- 2 \pi j n k /(2 N + 1)}.
\end{aligned}
\end{equation}

With 
\begin{equation}\label{eqn:HarmonicBalanceText:401c}
\Bx =
{\begin{bmatrix}
x_{-N} &
\cdots &
x_0 &
\cdots &
x_{N}
\end{bmatrix}}^\T,
\end{equation}
\begin{equation}\label{eqn:HarmonicBalanceText:401b}
\BX =
{\begin{bmatrix}
X_{-N} &
\cdots &
X_0 &
\cdots &
X_{N}
\end{bmatrix}}^\T,
\end{equation}
\begin{equation}\label{eqn:HarmonicBalanceText:401a}
W = e^{ 2 \pi j /(2 N + 1) },
\end{equation}
and \( \BF \) defined as
\begin{equation}\label{eqn:ece1254projectReport:401}
\small
\begin{bmatrix}
 W^{\lr{-N} \lr{-N} } &  W^{\lr{1-N} \lr{-N} }  & . &  W^{ \lr{N-1} \lr{-N} }  &  W^{\lr{ N} \lr{-N} } \\
 W^{\lr{-N} \lr{1-N} } &  W^{\lr{1-N} \lr{1-N} }  & . &  W^{ \lr{N-1} \lr{1-N} }  &  W^{\lr{ N} \lr{1-N} } \\
 .              &  .                      & .      & .                           &  .               \\
 W^{\lr{-N} \lr{N-1} } &  W^{\lr{1-N} \lr{N-1} }  & . &  W^{\lr{N-1} \lr{N-1} }  &  W^{\lr{N} \lr{N-1} } \\
 W^{\lr{-N} \lr{N} } &  W^{\lr{1-N} \lr{N} }  & . &  W^{\lr{N-1}\lr{ N} }  &  W^{\lr{N}\lr{ N} } \\
\end{bmatrix},
\end{equation}
the transform pair has the matrix form
\begin{subequations}
\begin{equation}\label{eqn:ece1254projectReport:402}
\Bx = \BF \BX
\end{equation}
\begin{equation}\label{eqn:ece1254projectReport:403}
\BX = \inv{2 N + 1} \overbar{\BF} \Bx,
\end{equation}
\end{subequations}

where \( \overbar{\BF} \) is the complex conjugate of \( \BF \).

\mysubsection{Harmonic Balance equations.}
\index{harmonic balance}
The time domain MNA equations for an RLC network with nonlinear sources can be put into the form
\begin{equation}\label{eqn:ece1254projectReport:420}
\BG \Bx(t) + \BC \dot{\Bx}(t) = \sum_{i = 1}^M \Bb_i u_i(t) + \sum_{i = 1}^S I_i \Bd_i w_i(t).
\end{equation}

With \( R \) equal to the number of unknowns, \( S \) equal to the number of nonlinear elements, \( M \) equal to the number of (positive, negative and zero) frequency sources, the dimensions of these matrices are \( \BG, \BC \in \text{\R{R \times R}} \), and \( \Bb_i, \Bd_i, \Bx(t) \in \text{\C{R \times 1}} \).

The linear sources can be grouped into a vector.
For example, with a constant source, and one source at the twice the fundamental frequency, this would be
\begin{equation}\label{eqn:ece1254projectReport:440}
\Bu(t) =
\begin{bmatrix}
e^{ -j \omega_0 2 t } \\
1 \\
e^{ j \omega_0 2 t }
\end{bmatrix},
\end{equation}

In an implementation of code that generates MNA equations it is convenient to store just the frequencies associated with these sources, which carries the same information without needing to encode the time dependence explicitly, as in
\begin{equation}\label{eqn:ece1254projectReport:460}
\Bu =
\begin{bmatrix}
- 2 \omega_0 \\
0 \\
2 \omega_0
\end{bmatrix}.
\end{equation}

For the nonlinear sources, the magnitudes \( I_i \) have been factored out, so that the vectors \( \Bd_i \) contain only the values \( 0, \pm 1 \).
The reasons for this will become clear in the example to follow.
\mysubsection{Frequency domain representation of MNA equations.}
\index{frequency domain}

Assuming a bandwidth limited periodic representation of the unknowns vector \( \Bx \) in \cref{eqn:ece1254projectReport:420}

\begin{equation}\label{eqn:ece1254projectReport:500}
x^{(i)}(t) = \sum_{n=-N}^N X_n^{(i)} e^{j n \omega_0 t},
\end{equation}
the derivative has the form

\begin{equation}\label{eqn:ece1254projectReport:520}
\dot{x}^{(i)}(t) = j \omega_0 \sum_{n=-N}^N X_n^{(i)} n e^{j n \omega_0 t}.
\end{equation}

The Fourier component representation of the MNA equations can be written as
\singleAndDoubleColumnVariations
{
\begin{equation}\label{eqn:ece1254projectReport:600}
\begin{aligned}
0 &= \sum_{n = -N}^N e^{j n \omega_0 t}
\Biglr{
   \lr{ \BG + j n \omega_0 \BC }
\begin{bmatrix}
X_n^{(1)} \\
X_n^{(2)} \\
\vdots \\
\end{bmatrix} \\
&\quad
-
[ \Bb_1 \cdots \Bb_M ]
\begin{bmatrix}
U_n^{(1)} \\
U_n^{(2)} \\
\vdots \\
\end{bmatrix}
- \sum_{i= 1}^S I_i \Bd_i
W_n^{(i)}
}.
\end{aligned}
\end{equation}
}
{
\begin{equation}\label{eqn:ece1254projectReport:600}
\begin{aligned}
0 &= \sum_{n = -N}^N e^{j n \omega_0 t}
\Biglr{
   \lr{ \BG + j n \omega_0 \BC }
\begin{bmatrix}
X_n^{(1)} \\
X_n^{(2)} \\
\vdots \\
\end{bmatrix} \\
&-
[ \Bb_1 \cdots \Bb_M ]
\begin{bmatrix}
U_n^{(1)} \\
U_n^{(2)} \\
\vdots \\
\end{bmatrix}
- \sum_{i = 1}^S I_i \Bd_i
W_n^{(i)}
}.
\end{aligned}
\end{equation}
}

Given the
assumption of periodicity, each of these equations must separately equal zero for each \( ( n, i ) \) pair, or
\singleAndDoubleColumnVariations
{
\begin{equation}\label{eqn:ece1254projectReport:540}
\begin{aligned}
\lr{
\BG + j n \omega_0 \BC
}
\begin{bmatrix}
V_n^{(1)} \\
V_n^{(2)} \\
\vdots \\
\end{bmatrix}
&=
[\Bb_i \cdots \Bb_M]
\begin{bmatrix}
U_n^{(1)} \\
U_n^{(2)} \\
\vdots \\
\end{bmatrix}
+ \sum_{i = 1}^S I_i \Bd_i
W_n^{(i)}
.
\end{aligned}
\end{equation}
}
{
\begin{equation}\label{eqn:ece1254projectReport:540}
\begin{aligned}
\lr{
\BG + j n \omega_0 \BC
}
\begin{bmatrix}
V_n^{(1)} \\
V_n^{(2)} \\
\vdots \\
\end{bmatrix}
&=
[\Bb_i \cdots \Bb_M]
\begin{bmatrix}
U_n^{(1)} \\
U_n^{(2)} \\
\vdots \\
\end{bmatrix} \\
&+ \sum_{i = 1}^S I_i \Bd_i
W_n^{(i)}
.
\end{aligned}
\end{equation}
}

%An implementation may group the vectors \( \Bd_i \) into a matrix, or even compute it on the fly
It is desirable to use an ordering of the vector components given by
\begin{equation}\label{eqn:ece1254projectReport:580}
\BI =
\begin{bmatrix}
I^{(1)}_{-N} \\
I^{(2)}_{-N} \\
%I^{(3)}_{-N} \\
\vdots \\
I^{(1)}_{1-N} \\
I^{(2)}_{1-N} \\
%I^{(3)}_{1-N} \\
\vdots \\
I^{(1)}_{N-1} \\
I^{(2)}_{N-1} \\
%I^{(3)}_{N-1} \\
\vdots \\
I^{(1)}_{N} \\
I^{(2)}_{N} \\
%I^{(3)}_{N} \\
\vdots \\
\end{bmatrix},
\qquad \BV =
\begin{bmatrix}
V^{(1)}_{-N} \\
V^{(2)}_{-N} \\
%V^{(3)}_{-N} \\
\vdots \\
V^{(1)}_{1-N} \\
V^{(2)}_{1-N} \\
%V^{(3)}_{1-N} \\
\vdots \\
V^{(1)}_{N-1} \\
V^{(2)}_{N-1} \\
%V^{(3)}_{N-1} \\
\vdots \\
V^{(1)}_{N} \\
V^{(2)}_{N} \\
%V^{(3)}_{N} \\
\vdots \\
\end{bmatrix}
\end{equation}

The superscript index of each element corresponds to the physical circuit node it is related to (i.e.
current flowing to the node, or the node voltage).
Those unknown physical circuit elements are generated using a modified nodal analysis procedure.
The subscript corresponds to the harmonic frequency, so the end result is as if there is a
node voltage for each harmonic frequency.

Ordering the vectors this way allows a block diagonal admittance matrix to be generated (\citep{giannini2004NonlinearMicrowaveCircuitDesign}
 \S 1.85)
, consisting of the admittance matrix that would be seen for each frequency.
%Peeter: I believe your code would give us the G and C matrices that arise from LMS formulations. This is
%great as we can generate each block Y(n?0) = G - jn?0C.
%Using the ordering of \cref{eqn:ece1254projectReport:580},
The structure of this equation is
%\begin{equation}\label{eqn:ece1254projectReport:1340}
\boxedEquation{eqn:ece1254projectReport:1340}{
\BY \BV = \BI + \calI(\BV),
}
%\end{equation}
where \( \BY \) is block diagonal
\begin{equation}\label{eqn:ece1254projectReport:1360}
\small{
\begin{bmatrix}
\BG + j \omega_0 (-N) \BC &        &                                   & \\
                            \qquad\ddots &                                   & \\
                                   &  \BG + j \omega_0 (0) \BC         & \\
                                   &                            \qquad\ddots & \\
                                   &                                   & \BG + j \omega_0 (N) \BC
\end{bmatrix}
}
\end{equation}
Construction of the HB stamps for the constant current \( \BI \) and nonlinear current \( \BI(\BV) \) is nicely illustrated by the example to follow.
\mysubsection{Example.  RC circuit with a diode.}
\index{rc circuit}
\index{diode}
\imageFigureHere{../figures/ece1254-multiphysics/diode-RC-current-source.pdf}{Simple nonlinear circuit.}{fig:diodeWithResistorAndCapacitorFig1}{3in}{proj:makefigures.m}
A diode with an RC load and AC current source is illustrated in
\cref{fig:diodeWithResistorAndCapacitorFig1}.
With
\begin{equation}\label{eqn:ece1254projectReport:640}
\begin{aligned}
\BG &=
\begin{bmatrix}
0 & 0 \\
0 & Z \\
\end{bmatrix}, \quad
\BC =
\begin{bmatrix}
0 & 0 \\
0 & C \\
\end{bmatrix}, \\
\Bv(t) &=
\begin{bmatrix}
v^{(1)}(t) \\
v^{(2)}(t) \\
\end{bmatrix}, \quad
\Bd =
\begin{bmatrix}
1 \\
-1
\end{bmatrix}, \quad
\Bb =
\begin{bmatrix}
1 \\
0
\end{bmatrix},
\end{aligned}
\end{equation}
the time domain equations for this circuit are
\singleAndDoubleColumnVariations
{
\begin{equation}\label{eqn:ece1254projectReport:1320}
\BG \Bv(t)
+\BC \Bv'(t)
=
\begin{bmatrix}
\Bb & -I_0 \Bd
\end{bmatrix}
\begin{bmatrix}
i_{\textrm{s}}(t) \\
1
\end{bmatrix}
+
I_0 \Bd
e^{ (v^{(1)}(t) - v^{(2)}(t))/V_{\textrm{T}}}.
\end{equation}
}
{
\begin{equation}\label{eqn:ece1254projectReport:1320}
\begin{aligned}
\BG \Bv(t)
+\BC \Bv'(t)
&=
\begin{bmatrix}
\Bb & -I_0 \Bd
\end{bmatrix}
\begin{bmatrix}
i_{\textrm{s}}(t) \\
1
\end{bmatrix} \\
&+
I_0 \Bd
e^{ (v^{(1)}(t) - v^{(2)}(t))/V_{\textrm{T}}}.
\end{aligned}
\end{equation}
}
For the exponential, let the transform pair be
\begin{equation}\label{eqn:ece1254projectReport:1300}
\epsilon(t) =
e^{ (v^{(1)}(t) - v^{(2)}(t))/V_{\textrm{T}}}
\simeq
\sum_{n=-N}^N E_n e^{ j n \omega_0 t }.
\end{equation}
This is an approximation at any times other than the Nykvist sampling times \( t_k = T k/(2 N + 1) \).

The frequency domain representation is
\begin{equation}\label{eqn:ece1254projectReport:820}
\lr{ \BG + j n \omega_0 \BC }
\begin{bmatrix}
V^{(1)}_n \\
V^{(2)}_n \\
\end{bmatrix}
=
-I_0 \Bd \delta_{n 0}
+
\Bb I^{(\textrm{s})}_n
+ I_0 \Bd E_n.
\end{equation}

To illustrate the stamping procedure for the constant current consider the \( N = 1 \) case.
The block matrix form for the constant current can be seen to be
\begin{equation}\label{eqn:ece1254projectReport:840}
\BI
=
\begin{bmatrix}
\Bb I^{(\textrm{s})}_{-1} \\
\Bb I^{(\textrm{s})}_{0} - I_0 \Bd \\
\Bb I^{(\textrm{s})}_{1} \\
\end{bmatrix}
\end{equation}
The nonlinear current \( \calI(\BV) \) needs to be examined further.
How much of this can be precomputed, and what is the simplest way to compute the Jacobian?

Continuing to use \( N = 1 \) to illustrate, a vector of Fourier components for the exponentials can be formed
\begin{equation}\label{eqn:ece1254projectReport:880}
\BE =
\begin{bmatrix}
E_{-1} \\
E_{0} \\
E_{1} \\
\end{bmatrix}, \quad
\Bepsilon =
\begin{bmatrix}
\epsilon_{-1} \\
\epsilon_{0} \\
\epsilon_{1} \\
\end{bmatrix},
\end{equation}
which allows the nonlinear current to be expressed as
\begin{equation}\label{eqn:ece1254projectReport:900}
\begin{aligned}
\calI &=
I_0
\begin{bmatrix}
\Bd E_{-1} \\
\Bd E_{0} \\
\Bd E_{1} \\
\end{bmatrix}
=
I_0
\begin{bmatrix}
\Bd \begin{bmatrix} 1 & 0 & 0 \end{bmatrix} \BE \\
\Bd \begin{bmatrix} 0 & 1 & 0 \end{bmatrix} \BE \\
\Bd \begin{bmatrix} 0 & 0 & 1 \end{bmatrix} \BE
\end{bmatrix} \\
&=
I_0
\begin{bmatrix}
\Bd & 0 & 0 \\
0 & \Bd & 0 \\
0 & 0 & \Bd
\end{bmatrix}
\BF^{-1} \Bepsilon
\end{aligned}
\end{equation}
In the last step \( \BE = \BF^{-1} \Bepsilon \) has been factored out (in its inverse Fourier form).
With
\begin{equation}\label{eqn:ece1254projectReport:920}
\BD =
\begin{bmatrix}
\Bd & 0 & 0 \\
0 & \Bd & 0 \\
0 & 0   & \Bd \\
\end{bmatrix},
\end{equation}
the current can now be written in a compact form
\boxedEquation{eqn:ece1254projectReport:960}{
\calI(\BV) =
I_0 \BD \BF^{-1} \Bepsilon(\BV),
}
%
however, an appropriate form for \( \Bepsilon \) must still be found.
\begin{equation}\label{eqn:ece1254projectReport:980}
\begin{aligned}
\Bepsilon &=
\begin{bmatrix}
\epsilon(t_{-1}) \\
\epsilon(t_{0}) \\
\epsilon(t_{1}) \\
\end{bmatrix}
=
\begin{bmatrix}
\exp\lr{ \lr{ v^{(1)}_{-1} - v^{(2)}_{-1} }/V_{\textrm{T}} } \\
\exp\lr{ \lr{ v^{(1)}_{0} - v^{(2)}_{0} }/V_{\textrm{T}} } \\
\exp\lr{ \lr{ v^{(1)}_{1} - v^{(2)}_{1} }/V_{\textrm{T}} }
\end{bmatrix} \\
\end{aligned}
\end{equation}
Each of the time domain voltage differences can be cast into vector form, and then expressed in the frequency domain.
For example
\begin{equation}\label{eqn:ece1254projectReport:1380}
\begin{aligned}
v^{(1)}_{-1} - v^{(2)}_{-1}
&=
\begin{bmatrix}
1 & 0 & 0
\end{bmatrix}
\lr{ \Bv^{(1)} - \Bv^{(2)} } \\
&=
\begin{bmatrix}
1 & 0 & 0
\end{bmatrix}
\BF
\lr{ \BV^{(1)} - \BV^{(2)} }
\end{aligned}
\end{equation}

This difference of Fourier component vectors needs to be expressed in terms of the big vector of unknown Fourier components \( \BV \)
\begin{equation}\label{eqn:ece1254projectReport:1000}
\begin{aligned}
\BV^{(1)} - \BV^{(2)}
&=
\begin{bmatrix}
V^{(1)}_{-1} - V^{(2)}_{-1} \\
V^{(1)}_{0} - V^{(2)}_{0} \\
V^{(1)}_{1} - V^{(2)}_{1} \\
\end{bmatrix}
\\ &=
\begin{bmatrix}
1 & -1 & 0 & 0 & 0 & 0 \\
0 &  0 & 1 & -1 & 0 & 0 \\
0 &  0 & 0 & 0 & 1 & -1 \\
\end{bmatrix}
\begin{bmatrix}
V_{-1}^{(1)} \\
V_{-1}^{(2)} \\
V_{0}^{(1)} \\
V_{0}^{(2)} \\
V_{1}^{(1)} \\
V_{1}^{(2)} \\
\end{bmatrix}.
\end{aligned}
\end{equation}
This is just
\begin{equation}\label{eqn:ece1254projectReport:1400}
\BV^{(1)} - \BV^{(2)}
=
\begin{bmatrix}
\Bd^\T & 0 & 0 \\
0 & \Bd^\T & 0 \\
0 & 0 & \Bd^\T \\
\end{bmatrix}
\BV
= \BD^\T \BV.
\end{equation}
With
%\begin{equation}\label{eqn:ece1254projectReport:1040}
\boxedEquation{eqn:ece1254projectReport:1040}{
\BH
=
\BF \BD^\T /V_{\textrm{T}}
=
\begin{bmatrix}
\Bh_1^\T \\
\Bh_2^\T \\
\Bh_3^\T
\end{bmatrix},
}
%\end{equation}
a compact matrix representation of the nonlinear current is now complete
\boxedEquation{eqn:ece1254projectReport:1080}{
\Bepsilon(\BV)
=
\begin{bmatrix}
e^{\Bh_1^\T \BV} \\
e^{\Bh_2^\T \BV} \\
e^{\Bh_3^\T \BV} \\
\end{bmatrix}.
}
%
\mysubsection{Other nonlinear sources.}
\index{voltage controlled voltage source!nonlinear}
The procedure above can be adapted without significant modification to other nonlinear current sources that are dependent on the difference in the voltages across their nodes (such as the power-law current sources specified by Pnnnnn lines in the HB netlist parser.)  A small change is required for voltage controlled sources.  Recall that the \( \BG \) stamp for such a source included a source line of the form
\begin{equation}\label{eqn:ece1254projectReport:1121}
\kbordermatrix{
        & n_{+} & n_{-} & nc_{+} & nc_{-} &  i_{n_{+},n_{-}} \\
\setlr{n_{+},n_{-}}        &    -1 & 1 & g & -g   & 0 \\
}.
\end{equation}
This must be changed to a pair of stamps, one for \( \BG \)
\begin{equation}\label{eqn:ece1254projectReport:1123}
\kbordermatrix{
        & n_{+} & n_{-} & nc_{+} & nc_{-} &  i_{n_{+},n_{-}} \\
\setlr{n_{+},n_{-}}        &    -1 & 1 & 0 & 0   & 0 \\
},
\end{equation}
and one for the corresponding line of \( \calI(\BV) \).  If that line is \( m \), let
\begin{equation}\label{eqn:ece1254projectReport:1122}
\Bd_{\textrm{e}} =
{[\delta_{i,m}]}_i,
\end{equation}
so that nonlinear current for this source is
\begin{equation}\label{eqn:ece1254projectReport:1124}
\calI(\BV) = -g
\begin{bmatrix}
\Bd_{\textrm{e}} &  &  \\
   & \Bd_{\textrm{e}} & \\
   &  & \ddots
\end{bmatrix}
\BF^{-1}
\begin{bmatrix}
 \lr{\Bh_1 \BV}^k \\
 \lr{\Bh_2 \BV}^k \\
 \vdots
\end{bmatrix},
\end{equation}
where the \( \Bh_i \) are constructed with \cref{eqn:ece1254projectReport:920} as done for exponential nonlinear currents.
\mysubsection{Jacobian.}
\index{Jacobian}
With a compact matrix representation of the nonlinear current, attention can now be turned to the Jacobian of the nonlinear current.
Let \( \BA = I_0 \BD \BF^{-1} = [ a_{ij} ]_{ij} \).
The coordinates of the nonlinear current (with summation implied) are
\begin{equation}\label{eqn:ece1254projectReport:1120}
\begin{aligned}
\calI_i &= 
a_{ik} \epsilon_k= a_{ik} \exp\lr{ \Bh_k^\T \BV }.
\end{aligned}
\end{equation}
so the Jacobian components are
\begin{equation}\label{eqn:ece1254projectReport:1140}
\begin{aligned}
[\BJ^{\calI}]_{ij}
&= a_{ik} \epsilon_k
\\ &= a_{ik} \PD{V_j}{} \exp\lr{ \Bh_k^\T \BV }
\\ &= a_{ik} h_{kj} \exp\lr{ \Bh_k^\T \BV }.
\end{aligned}
\end{equation}

Back in matrix form, after factoring out the \( \BA \), the Jacobian is
\begin{equation}\label{eqn:ece1254projectReport:1141}
\BJ^{\calI}
=
\BA
{\begin{bmatrix}
h_{ij}
\exp\lr{ \Bh_i^\T \BV }
\end{bmatrix}}_{ij}.
\end{equation}
The second factor has all the rows of \( \BH \), but each multiplied by an exponential, so the Jacobian for the nonlinear current
is now completely determined
\boxedEquation{eqn:ece1254projectReport:1200}{
\BJ^{\calI}( \BV ) =
I_0 \BD \BF^{-1}
\begin{bmatrix}
\Bh_1^\T \exp\lr{ \Bh_1^\T \BV } \\
\Bh_2^\T \exp\lr{ \Bh_2^\T \BV } \\
\vdots
\end{bmatrix}.
}

A quick sanity check of dimensions seems worthwhile, and shows that all is well
\begin{itemize}
\item \( \BA \) : \( R(2 N + 1) \times 2 N + 1 \)
\item \( \BU \) : \( 2 N + 1 \times R(2 N + 1) \)
\item \( \BJ^{\calI} \) : \( R(2 N + 1) \times R(2 N + 1) \)
\end{itemize}

This derivation was done exclusively with the exponential nonlinearity that is found in the diode model, but is easy to generalize.
If the nonlinearity is of the form
\( I_0 g((v^{(1)} - v^{(2)})/V_{\textrm{T}}) \), then the nonlinear current and its Jacobian will have the same structural form
\begin{subequations}
\begin{equation}\label{eqn:ece1254projectReport:1440}
\calI(\BV)
=
I_0 \BD \BF^{-1}
\begin{bmatrix}
\Bh_1^\T g(\Bh_1^\T \BV) \\
\Bh_2^\T g(\Bh_2^\T \BV) \\
\vdots
\end{bmatrix}
\end{equation}
\begin{equation}\label{eqn:ece1254projectReport:1420}
\begin{aligned}
\BJ^{\calI}( \BV ) &=
I_0 \BD \BF^{-1}
\begin{bmatrix}
\Bh_1^\T \evalbar{ g'(x)}{x= \Bh_1^\T \BV } \\
\Bh_2^\T \evalbar{ g'(x)}{x= \Bh_2^\T \BV } \\
\vdots
\end{bmatrix}.
\end{aligned}
\end{equation}
\end{subequations}
Note that for voltage controlled voltage sources the leading \( I_0 \BD \BF^{-1} \) factor must be modified as specified in \cref{eqn:ece1254projectReport:1124}.
\mysubsection{Newton's method solution.}
\index{Newton's method}
All the pieces required for a Newton's method solution are now in place.
The goal is to find a value of \( \BV \) that provides the zero
\begin{equation}\label{eqn:ece1254projectReport:1220}
f(\BV) = \BY \BV - \BI - \calI(\BV).
\end{equation}
Expansion to first order around an initial guess \( \BV^0 \), gives
\begin{equation}\label{eqn:ece1254projectReport:1240}
f( \BV^0 + \Delta \BV ) = f(\BV^0) + \BJ(\BV^0) \Delta \BV \approx 0,
\end{equation}
where the full Jacobian of \( f(\BV) \) is
\begin{equation}\label{eqn:ece1254projectReport:1260}
\BJ(\BV) = \BY - \BJ^{\calI}(\BV).
\end{equation}

The Newton's method refinement of the initial guess follows by inversion
\begin{equation}\label{eqn:ece1254projectReport:1280}
\Delta \BV = -\lr{ \BY - \BJ^{\calI}(\BV^0) }^{-1} f(\BV^0).
\end{equation}

The nonlinearities of some of the circuits that were solved as tests required source stepping, as the Jacobian grew beyond the representable range of the double matrix representation in some cases and became uninvertible.
A heuristic that was found effective was to half the step size for each Newton's method iteration that the Jacobian became ill-conditioned, and to double it for each iteration where convergence was achieved.
\mysubsection{Alternative handling of the nonlinear currents and Jacobians.}
Two methods for the calculation of the nonlinear currents and Jacobians were developed while coding this HB implementation.
Our first implementation used
a ``big \( \BF \)'' DFT matrix, of size \( ( R (2 N + 1 ) \times R (2 N + 1) ) \).
With the elements of \( \BV \) grouped by frequency, a DFT matrix that operated on \( \BV \) was formed where
every value in each frequency group was multiplied by the same value in the DFT matrix.
The nonlinear currents \( \calI(\Bv) \), produced at each time sample to be determined were calculated.
Using the
``big'' DFT of these currents, the discrete frequency
components of the nonlinear currents were obtained.
Each nth-frequency component of this current would be added
to the equation of the nth-set of MNA equations corresponding to the nodes the element is connected
to.

This nonlinear current calculation was later refined to work completely in the frequency domain, precalculating the matrices \( \BF \BD^\T/V_{\textrm{T}} \) and \( I_0 \BD \BF^{-1} \) for each nonlinear source, where \( \BF \) was the (standard) DFT matrix with dimensions \( (2 N + 1 ) \times (2 N + 1 ) \), as described above.
\section{Results.}
\mysubsection{Low pass filter.}
\index{low pass filter}
%\subfloat[Parameters]{%
%\small
%\begin{tabular}[b]{|l|r|} \hline
%  Parameter & Value \\ \hline\hline
%  \( \nu_0 \) & 1 MHz \\ \hline
%  \( V_1 \) & 10 V (peak) \\ \hline
%  \( R_1 \) & 1 k \(\Omega\) \\ \hline
%  \( C_1 \) & 1 \(\mu\) F \\ \hline
%  \( D_1, I_0 \) & 10 pA \\ \hline
%  \( D_1, V_{\textrm{T}} \) & 25 mV \\ \hline
%\end{tabular}
%\label{tab:lowPassFilter:lowPassFilterParameters}}
\imageFigureHere{../figures/ece1254-multiphysics/lowPassFilter.pdf}{Low pass filter.}{fig:lowPassFilterFig1}{3.0in}{proj:makefigures.m}
A low pass filter with a cut off frequency of 5 MHz is shown in
\cref{fig:lowPassFilterFig1}.
With
inputs
ranging from 1 to 8 Mhz, this serves as a basic sanity test of the HB linear network solver.
The time domain input signal and the corresponding filtered output, with a visible reduction in higher frequency components,
is shown in \cref{fig:lowPassFilterFig2}.
\imageFigureHere{../figures/ece1254-multiphysics/report/lowPassFilterSourceAndOutputVoltagesFig2.pdf}{Source and output.}{fig:lowPassFilterFig2}{3.0in}{proj:makefigures.m}

The magnitudes of the frequency domain components in the input signal are shown in \cref{fig:lowPassFilterFig3}.
\imageFigureHere{../figures/ece1254-multiphysics/report/lowPassFilterInputVoltageSpectrumFig3.pdf}{Input voltage spectrum.}{fig:lowPassFilterFig3}{3.0in}{proj:makefigures.m}

The output spectrum, seen in \cref{fig:lowPassFilterFig4}, shows attenuation at all frequencies.
However, all the higher frequency contributions are cut down significantly, as expected.
\imageFigureHere{../figures/ece1254-multiphysics/report/lowPassFilterOutputVoltageSpectrumFig4.pdf}{Output voltage spectrum.}{fig:lowPassFilterFig4}{3.0in}{proj:makefigures.m}
\mysubsection{Half wave rectifier.}
\index{half wave rectifier}
The half wave rectifier of \cref{fig:typicalRectifierCircuitFig1}, was suggested (with \( C = 0 \)) as a first test of the nonlinear HB system.
%\imageFigureHere{../figures/ece1254-multiphysics/halfWaveRectifier}{Half wave rectifier circuit.}{fig:HalfWaveRectifierFig1}{3in}{proj:makefigures.m}
\imageFigureHere{../figures/ece1254-multiphysics/halfWaveRectifierWithCap.pdf}{Half wave rectifier circuit.}{fig:typicalRectifierCircuitFig1}{3.0in}{proj:makefigures.m}
Many of the circuits considered in this report were solved using source stepping of the modified system
\begin{equation}\label{eqn:ece1254projectReport:1460}
f(\BV, \lambda) = \BY \BV - \lambda \BI - \calI(\BV),
\end{equation}
with \( \lambda \in [0, 1] \), each time feeding the previous solution of \( \BV \) as the initial guess with the new source step value \( \lambda \).
That did not work for this system, which is ill-conditioned for \( \lambda = 0 \).
To handle this sort of ill conditioning the Newton's method procedure was modified to use \( \lambda \in [d\lambda, 1] \) if the Jacobian for \cref{eqn:ece1254projectReport:1460} was ill conditioned or NaN at \( \lambda = 0 \).
With that refinement of the source stepping algorithm, HB simulation of this half wave rectifier converges.
%produces the expected results, as plotted in
Adding a capacitor to the half wave rectifier, adds an imaginary component to the admittance matrix \( \BY \), while still only requiring a single nonlinear component, where the use of a large enough capacitor value produces a DC rectified signal from an AC input.
Plots of the simulation results with
capacitor values of \( 0, 1 \, \text{nF}, 10 \, \text{nF}, 1000 \text{nF} \) are plotted in \cref{fig:HalfWaveRectifierFig2}.
\imageFigureHere{../figures/ece1254-multiphysics/report/halfWaveRectifierRangeOfCapValuesMultiSource.pdf}{Half wave rectifier voltages.}{fig:HalfWaveRectifierFig2}{3in}{proj:makefigures.m}
%\mysubsection{AC to DC conversion}
%\cref{fig:typicalRectifierCircuitFig3}, and
%\cref{fig:typicalRectifierCircuitFig4} respectively.
%Further reduction of the capacitor size would eventually produce the original half wave rectifier results.
%\imageFigureHere{../figures/ece1254-multiphysics/report/typicalRectifierCircuitSourceAndOutputVoltagesFig2}{DC rectifier voltages at \( 1 \mu \) F.}{fig:typicalRectifierCircuitFig3}{3.0in}{proj:makefigures.m}
%\imageFigureHere{../figures/ece1254-multiphysics/report/typicalRectifierCircuitSmallerCapSourceAndOutputVoltagesFig2}{DC rectifier voltages at \( 1 \) nF.}{fig:typicalRectifierCircuitFig4}{3.0in}{proj:makefigures.m}
\mysubsection{Bridge rectifier.}
\index{bridge rectifier}
The bridge rectifier circuit of
\cref{fig:bridgeRectifierCapFilterFig1}, in configurations with zero and non-zero capacitance values,
was used as a test system that included multiple nonlinear elements.
%\imageFigureHere{../figures/ece1254-multiphysics/bridgeRectifier}{Bridge rectifier.}{fig:bridgeRectifierCircuitDiagramFig1}{2.5in}{proj:makefigures.m}
\imageFigureHere{../figures/ece1254-multiphysics/bridgeRectifierWithCap.pdf}{Bridge rectifier.}{fig:bridgeRectifierCapFilterFig1}{2.5in}{proj:makefigures.m}
Current will flow through only one pair of the diodes at any given point.
In the zero-capacitance configuration, using a three frequency source, inversion of any negative source voltages is observed in \cref{fig:bridgeRectifierFig2} at the output node.
\imageFigureHere{../figures/ece1254-multiphysics/report/bridgeRectifierSourceAndOutputVoltagesFig2.pdf}{Source and output voltages.}{fig:bridgeRectifierFig2}{3.0in}{proj:makefigures.m}
%\begin{figure}[!h]
%%\subfloat[Parameters]{%
%%\small
%%\begin{tabular}[b]{|l|r|} \hline
%%  Parameter & Value \\ \hline\hline
%%  \( \nu_0 \) & 1 MHz \\ \hline
%%  \( V_1 \) & 10 V (peak) \\ \hline
%%  \( R_1 \) & 1 k \(\Omega\) \\ \hline
%%  \( C_1 \) & 1 \(\mu\) F \\ \hline
%%  \( D_1, I_0 \) & 10 pA \\ \hline
%%  \( D_1, V_{\textrm{T}} \) & 25 mV \\ \hline
%%\end{tabular}
%\end{figure}
\mysubsection{Cpu time and error vs N.}
The bridge rectifier circuit, with \( C1 = 0 \) was also used for
measurements of error and cpu times for range of values of \( N \).
These are plotted as a log-log plots in
\cref{fig:bridgeRectifierFig5}.
Here error is measured as the biggest absolute difference between the solution for a given value of \( N \) vs
a ``sufficiently large'' value of N (100).
\imageFigureHere{../figures/ece1254-multiphysics/report/bridgeRectifierErrorAndCpuTimesFig5.pdf}{Error and Cpu vs N.}{fig:bridgeRectifierFig5}{3.0in}{proj:makefigures.m}

Observe that the relations are both nearly linear in the log-log plots, indicating that the error and cpu both have the form \( a N^m \).
Measuring the respective slopes (using the matlab function polyfit), the cpu time is found to be \( O(N^{1.5}) \), and the error is \( O(1/N^2) \), at least for this particular circuit.
%\begin{figure*}[!h]
%\subfloat[Parameters]{%
%\tiny{
%\begin{tabular}[b]{|l|r|} \hline
%  Parameter & Value \\ \hline\hline
%  \( \nu_0 \) & 1 MHz \\ \hline
%  \( I_0 \) & 10 mA (peak) \\ \hline
%  \( R_1 \) & 1 k \(\Omega\) \\ \hline
%  \( D_1, I_0 \) & 10 pA \\ \hline
%  \( D_1, V_{\textrm{T}} \) & 25 mV \\ \hline
%\end{tabular}
%}
%\subfloat[Output.]{
%   \includegraphics[width=1.5in]{../figures/ece1254-multiphysics/report/simpleRectifierCircuitVoltageFig2}%
%%   \def\svgwidth{1.5in}
%%   \input{../figures/ece1254-multiphysics/report_svg/simpleRectifierCircuitVoltageFig2.pdf_tex}%
%\label{fig:simpleRectifierCircuitFig2}}
%\end{figure*}
\mysubsection{Power law nonlinearities.}
As a test for the implementation's support of \textbf{Plabel} form conductance models (voltage differences scaled and raised to a power), two of the diodes in the bridge rectifier circuit were replaced with an 8 term Taylor approximation
\begin{equation}\label{eqn:ece1254projectReport:1480}
\begin{aligned}
i_{\textrm{d}}(t) &= I_0 \lr{ e^{(v^{(n_1)} - v^{(n_2)})/V_{\textrm{T}}} - 1}
\approx
I_0 \sum_{k =1}^8
\inv{k!} \lr{(v^{(n_1)} - v^{(n_2)})/V_{\textrm{T}}}^k
\end{aligned}
\end{equation}
The results are plotted in \cref{fig:bridgeRectifierPowSourceAndOutputVoltagesFig2} for a single frequency source.
\imageFigureHere{../figures/ece1254-multiphysics/report/bridgeRectifierPowSourceAndOutputVoltagesFig2.pdf}{Source and output voltages for partial approximate bridge rectifier.}{fig:bridgeRectifierPowSourceAndOutputVoltagesFig2}{3.0in}{proj:makefigures.m}

This is not a numerically smart approximation as computation of the nonlinear currents and the Jacobian requires much more work.
Those calculations also require the multiplication of pairs of very large and small numbers (the exponents and the inverse factorials) introducing more chance for round off error than using the exponential directly.
When possible, as with the exponential, encoding direct support for the functional form of a nonlinear conductance in the current and Jacobian calculation code is likely always superior.
%Despite that, there are likely some conductance models that have a good approximation as a set of power series terms.

As a test for the nonlinear voltage controlled \textbf{Elabel} elements, the circuit of \cref{fig:squareAndCubeAmp} containing a pair of simple nonlinear amplifiers was simulated with a single frequency source, with results plotted in \cref{fig:PowerLawAmplifiersvoltageFig1}.
\imageFigureHere{../figures/ece1254-multiphysics/report/squareAndCubeAmp.pdf}{Non-linear power-law amplifiers.}{fig:squareAndCubeAmp}{3.0in}{proj:makefigures.m}
\imageFigureHere{../figures/ece1254-multiphysics/report/PowerLawAmplifiersvoltageFig1.pdf}{Non-linear power-law amplifiers voltages.}{fig:PowerLawAmplifiersvoltageFig1}{3.0in}{proj:makefigures.m}

Observe that the time domain output of the squaring-amplifier
has the expected doubling of the fundamental frequency, and is shifted into the \( v > 0 \) region consistent with
\begin{equation}\label{eqn:HarmonicBalanceText:17}
\cos^2(\omega_0 t) = \lr{\cos(2 \omega_0 t) + 1}/2.
\end{equation}
\mysubsection{Stiff systems.}
\index{stiff}
For the low pass and rectifier circuits above, use of source stepping was sufficient (and required) to ensure convergence.
However, introduction of a non-zero value for the capacitor of the bridge rectifier circuit \cref{fig:bridgeRectifierCapFilterFig1} is enough to cause a number of Matlab messages ``\matlabText{Matrix is close to singular or badly scaled}'' along with associated very small conditioning numbers.
%\begin{figure}[!h]
%\subfloat[Parameters]{%
%\tiny{
%\begin{tabular}[b]{|l|r|} \hline
%  Parameter & Value \\ \hline\hline
%  \( \nu_0 \) & 1 MHz \\ \hline
%  \( V_1 \) & 3.5 V (peak) \\ \hline
%  \( R_1 \) & 1 k \(\Omega\) \\ \hline
%  \( C_1 \) & 1 \(\mu\) F \\ \hline
%  \( D_i, I_0 \) & 10 pA \\ \hline
%  \( D_i, V_{\textrm{T}} \) & 25 mV \\ \hline
%\end{tabular}
%}
%\label{tab:bridgeRectifierCapFilterParameters}}
%\end{figure}
The circuit simulation produces the expected a DC rectified signal as plotted in \cref{fig:bridgeRectifierCapFilterFig2}.
\imageFigureHere{../figures/ece1254-multiphysics/report/bridgeRectifierCapFilterSourceAndOutputVoltagesFig2.pdf}{Source and output voltages.}{fig:bridgeRectifierCapFilterFig2}{3.0in}{proj:makefigures.m}

Artifacts of the conditioning problems show up in the simulated outputs for the voltages through the diodes as plotted in \cref{fig:bridgeRectifierCapFilterFig3} for a single frequency source.
\imageFigureHere{../figures/ece1254-multiphysics/report/bridgeRectifierCapFilterDiodeVoltagesFig4.pdf}{Diode voltages.}{fig:bridgeRectifierCapFilterFig3}{3.0in}{proj:makefigures.m}

The doubling and halfing heuristic used for the step sizes is insufficient for this circuit.
Even after
reducing the step size by a many orders of magnitude, and capping the step doubling at the initial step size, the solutions for this circuit stubbornly refused to converge without Matlab numerical conditioning diagnostics.
For this circuit, and surely others, it appears that additional strategies beyond source stepping are required to compute a solution without hitting numerical instabilities.
\section{Conclusion.}
When the steady state solution of a uniformly discretized linear or nonlinear system with periodic sources is desired, the Harmonic Balance method is a natural way to express the system equations, and once solved provides both the frequency and time domain representation of the system response.  Some of the nonlinear circuits for which HB solutions were sought proved troublesome, even with source stepping and fine mesh sizes, so it is clear that developing a generic and numerically robust HB system would require additional study and careful implementation.
\appendices
\myappendix{Discrete Fourier Transform inversion}
\label{appendix:discreteFourierInversion}
To find \( X_n \) evaluate the sum
\begin{equation}\label{eqn:ece1254projectReport:80}
\begin{aligned}
\sum_{k = -N}^N &x(t_k) e^{-j m \omega_0 t_k} \\
&=
\sum_{k = -N}^N
\lr{
\sum_{n = -N}^N X_n e^{ j n \omega_0 t_k}
}
e^{-j m \omega_0 t_k} \\
&=
\sum_{n = -N}^N X_n
\sum_{k = -N}^N
e^{ j (n -m )\omega_0 t_k}
\end{aligned}
\end{equation}
This interior sum has the value \( 2 N + 1 \) when \( n = m \).
For \( n \ne m \), and
\( a = e^{j (n -m ) \frac{2 \pi}{2 N + 1}} \), this is
\begin{equation}\label{eqn:ece1254projectReport:100}
\begin{aligned}
\sum_{k = -N}^N
e^{ j (n -m )\omega_0 t_k}
&=
\sum_{k = -N}^N
e^{ j (n -m )\omega_0 \frac{T k}{2 N + 1}}
\\ &=
\sum_{k= -N}^N a^k
\\ &=
a^{-N} \sum_{k= -N}^N a^{k+ N}
\\ &=
a^{-N} \sum_{r= 0}^{2 N} a^{r}
\\ &=
a^{-N} \frac{a^{2 N + 1} - 1}{a - 1}.
\end{aligned}
\end{equation}
Since \( a^{2 N + 1} = e^{2 \pi j (n - m)} = 1 \), this sum is zero when \( n \ne m \).
This means that
\begin{equation}\label{eqn:ece1254projectReport:120}
\sum_{k = -N}^N
e^{ j (n -m )\omega_0 t_k} = (2 N + 1) \delta_{n,m}.
\end{equation}
Substitution back into \cref{eqn:ece1254projectReport:80} proves the Fourier inversion relation \cref{eqn:ece1254projectReport:60}.