3. Calculating the HTFs in Matlab

For this simple differential equation, the HTFs can also be computed in Matlab using the methods described in [1]. To do so, first, we describe the time-variation of the $a(t)$-coefficient as a Fourier series:

To obtain the HTFs at a frequency $f$, we construct the following large matrix

with

and

where $s=j2\pi f$. This $\boldsymbol{\mathcal{A}}$-matrix is a band matrix. In our case, the multisine has $10$ excited frequency lines, so the $\boldsymbol{\mathcal{A}}$-matrix has $10$ non-zero numbers above and below its diagonal. Every other number in the matrix is zero. Theoretically, the $\boldsymbol{\mathcal{A}}$-matrix should be infinitely large, in practice, we will use a large square matrix of size $N\times N$. When $N$ is significantly larger than the amount of HTFs that needs to be calculated, the error made by truncating $\boldsymbol{\mathcal{A}}$ is limited [1]. We need to calculate $201$ HTFs, so $N$ was set to $401$ to minimise errors.

The inverse of the $\boldsymbol{\mathcal{A}}$-matrix is the Wereley matrix for this system [1] [2], which will be denoted as $\boldsymbol{\mathcal{G}}\left(s\right)$

The Wereley matrix contains the HTFs of the circuit in the following arrangement:

so by selecting the right column out of $\boldsymbol{\mathcal{G}}\left(s\right)$, the HTFs of the system are obtained.