我们假设所建模的系统是线性且时不变的，输入为 $u(t)\in {\mathbb{R}}^m$ ，输出为 $y(t)\in {\mathbb{R}}^q$。由于线性和时不变性，输出可以表示为输入 $u(t)$ 和系统的冲激响应 $h(t)\in {\mathbb{R}}^{q\times m}$ 的卷积：

$$y(t)={\int }_{-\infty }^{+\infty }h(t-\tau )u(\tau )d\tau (1.1)$$

将两边进行拉普拉斯变换，得到

$$Y(s)=H(s)U(s)(1.2)$$

其中 $s=\sigma +\text{j}\omega$ 为复频率。在式 (1.2) 中， $U(s)\in {\mathbb{C}}^m$ 和 $Y(s)\in {\mathbb{C}}^q$ 分别为 $u(t)$ 和 $y(t)$ 的拉普拉斯变换，而 $H(s)\in {\mathbb{C}}^{q\times m}$ 为系统的传递函数。VF 算法解决以下问题：给定传递函数的 $k$ 个测量值

$$H_k=H(\text{j}{\omega }_k),k=1,\dots ,k,(1.3)$$

确定一个有理函数 $$\tilde{H}(s)$$ 来逼近给定的测量值

$$\tilde{H}(\text{j}{\omega }_k)\approx H_k\forall k=1,\dots ,k.(1.4)$$

在 VF 中，选择 $\tilde{H}(s)$ 为一个有理函数。有理函数是通用的逼近器，因此可以在任意精度下逼近广义的函数。此外，由于集总系统的传递函数本质上是有理的，因此这是建模动态系统的自然选择。最后，有理函数可以表示为状态矩阵、极点-残差形式、一组微分方程、等效电路等多种形式。这种灵活性有助于将简化模型集成到现有的计算数学和系统仿真软件中。

问题(1.4)是一个需要通过数值方法解决的频率响应估计问题。早在1950年代，Levy、Sanathanan 和 Koerner等人就开始了这方面的工作，这为后来的 VF 算法提供了基础。为了简化问题，文章首先考虑一个 单输入单输出系统（即输入和输出各只有一个信号），后面会讨论更一般的情况。

选择合适的模型

为了解决频率响应估计问题，我们首先需要选择一个合适的 数学模型 来描述系统的行为。这个模型通常是传递函数$\tilde{H}(s)$，也就是系统的输入输出关系。为了建立模型，最常见的选择是使用两个多项式的比值形式，即：

$$\tilde{H}(s)=\frac{n(s)}{d(s)}=\frac{{\sum }_{n=0}^N{a}_n{s}^n}{{\sum }_{n=0}^N{b}_n{s}^n}（1.5）$$

其中，$a_n$和$b_n$是待估计的系数，$n$是多项式的阶数。为了简化计算，通常可以将其中一个系数标准化（比如设定$b_n=1$）。这个选择的多项式阶数$n$是为了保证系统的传递函数在频率$s\to \infty$时是有界的。

在某些应用中（例如被动电路的建模），传递函数的阶数可能会随频率$s$增长，这时候需要对模型的阶数做出调整（比如增加分子的阶数）。

最小化误差

一旦确定了模型的形式，接下来的任务就是通过给定的样本数据来 估计模型参数（即$a_n$和$b_n$的值）。为了做到这一点，我们需要选择一个适当的 误差度量，也就是我们要最小化的目标函数。文章中选择的是 二范数（$l_2$范数），即最小化样本响应$H_k$和模型响应$\tilde{H}(j{\omega }_k)$之间的差异。

具体来说，我们希望最小化的目标函数是：

$$e^2=\frac{1}{K}\sum _{k=1}^K{\left|H_k-\tilde{H}(j{\omega }_k)\right|}^2（1.6）$$

这是一个标准的最小二乘问题，其中$H_k$是样本数据，$\tilde{H}(j{\omega }_k)$是根据模型计算出的频率响应。

非线性最小二乘问题

问题的复杂之处在于，分母中的系数$b_n$是未知的，这使得目标函数是一个 非线性最小二乘问题。虽然可以使用非线性优化算法来求解这个问题，但这种方法通常比较耗时，且容易陷入局部最小值。

线性化求解

为了解决这个问题，提出了一种更加高效的方法：将非线性问题 线性化 成一个线性最小二乘问题。这可以通过 QR分解 等算法来高效求解。线性化的过程将使得求解过程更为简洁，并且避免了非线性优化中可能出现的陷阱。

我们首先将(1.6)重写为：

$$e^2=\frac{1}{K}\sum _{k=1}^K{\left|\frac{H_k{\sum }_{n=0}^n{b}_n(ȷ{\omega }_k{)}^n-{\sum }_{n=0}^n{a}_n(ȷ{\omega }_k{)}^n}{{\sum }_{n=0}^n{b}_n(ȷ{\omega }_k{)}^n}\right|}^2.（1.7）$$

Levy 提出通过简单地忽略分母来将(1.7)线性化，并最小化以下目标函数：

$$(e_L{)}^2=\frac{1}{k}\sum _{k=1}^K{\left|H_k\sum _{n=0}^n{b}_n(ȷ{\omega }_k{)}^n-\sum _{n=0}^n{a}_n(ȷ{\omega }_k{)}^n\right|}^2.（1.8）$$

最终，这个过程简化为求解一个线性方程组的最小二乘问题。不幸的是，这个简单的技巧通常无法提供(1.4)的准确解。事实上，误差泛函(1.7)和(1.8)仅在

$$\sum _{n=0}^n{b}_n(ȷ\omega {)}^n$$

为常数时等价，而这种情况一般很少发生。

为了解决这个问题，Sanathanan和Koerner提出了一种迭代过程来提高解的精度。在第一次迭代（即$i=1$）中，最小化 Levy 泛函(8.8)，得到模型系数的首次估计，我们将其记为$a_n^{(1)}$和$b_n^{(1)}$。在随后的迭代中（即$i\ge 2$），最小化以下(8.7)的线性化形式：

$$(e_{SK}^{(i)}{)}^2=\frac{1}{k}\sum _{k=1}^K{\left|\frac{H_k{\sum }_{n=0}^n{b}_n^{(i)}(ȷ{\omega }_k{)}^n-{\sum }_{n=0}^n{a}_n^{(i)}(ȷ{\omega }_k{)}^n}{{\sum }_{n=0}^n{b}_n^{(i-1)}(ȷ{\omega }_k{)}^n}\right|}^2（1.9）$$

从而得到模型系数$a_n^{(i)}$和$b_n^{(i)}$的新估计。可以看出，在(8.9)中，使用了上一轮迭代中的系数$b_n^{(i-1)}$来近似(8.7)中的“非线性”项。由于未知数$a_n^{(i)}$和$b_n^{(i)}$只出现在分子中，Sanathanan–Koerner方法只需要求解线性最小二乘问题。

可以看到，在(1.9)中，分母中的项${\sum }_{n=0}^n{b}_n^{(i-1)}(ȷ{\omega }_k{)}^n$起到了频率相关的权重作用。该权重旨在逐步消除(1.6)线性化过程中引入的偏差。对于离散时间系统，Sanathanan–Koerner方法的对应方法由Steiglitz和McBride提出。