Non-parametric frequency response identification for linear systems leverages two crucial facts.

1. Any reasonable signal \(y(t)\) which is \(T\)-periodic (meaning \(y(t) = y(t+T)\)) can be expressed as a linear combination of the complex exponentials \(\varphi_k(\w_o) = e^{j\w_okt}\). Moreover, the \(\varphi_k(\w_o)\) form an orthogonal basis over \(\mathbb{T}_{T_o}\).

1. For a linear system with impulse response \(h(t)\), \(\varphi_k\) is an eigenfunction. The eigenvalue associated with \(\varphi_k(\w_o)\) is the systems frequency response at \(\w_o\), \(H(j\w_o)\).

It is rather involved to prove that the \(\varphi_k(\w_o)\) form a basis for \(\mathbb{T}_T\) ¹. It is, however, pretty easy to show that they are orthogonal, by which we mean that \(<\varphi_k, \varphi_{\ell}\rangle = \delta_{k\ell}\). Or, explicitly using the inner product over \(\mathbb{T}_T\), we have

\begin{align} \langle \varphi_k, \varphi_{\ell}\rangle &= \frac{1}{T}\int_{T} e^{-j\w_okt}e^{j\w_o\ell t} dt\\ &= \frac{1}{T}\int_{T} e^{j\w_ot (\ell - k)} dt\\ &= \begin{cases} 1 & k=\ell\\ 0 & k\neq\ell. \end{cases} \end{align}

Now, let's first think about the idea that \(\varphi_k(\w_o)\) is an eigenfunction of \(h(t)\). If we drive \(h(t)\) with a complex exponential of frequency \(\w_o\), the output is given by the convolution

\begin{align} h(t)*e^{j\w_ot} &=\int_{\mathbb{R}} h(\tau)e^{j\w_o(t-\tau)}d\tau\\ &= H(j\w_o)e^{j\w_ot}. \end{align}

Since \(h(t)\) is a linear system, and convolution is a linear operator, we can represent this as the linear mapping \(\mathcal{A}:~\varphi_k\mapsto \Hjw \varphi_k\) or

\begin{equation} \mathcal{A}e^{j\w_ot} = H(j\w_o)e^{j\w_ot}. \end{equation}

In system identification, our goal to to find \(H(j\w_o)\). To do this, we experimentally let our system operate on the eigenfunction \(e^{j\w_ot}\) and try to determine the eigenvalue \(H(j\w_o)\). Of course, we can't measure \(H(j\w_o)\) directly. What we measure instead is the system output

\begin{equation} y(t) = e^{j\w_ot} H(j\omega).\label{eqn:out1} \end{equation}

How do we get \(\Hjw\) from \eqref{eqn:out1}? We recognize that \(\ejw\) is a member of the Fourier basis. It follows that the only non-zero coefficient in the Fourier series expansion of \(y(t)\) is \(y_1=\Hjw\). In other words,

\begin{align} y(t) & = \sum_{k=-\infty}^{\infty} y_k e^{j\w_okt} \\ y_k= &\begin{cases} H(\pm j\w_o), & k=\pm1\\ 0, & \text{else} \end{cases} \end{align}

So our job is calculate the first Fourier coefficient of \(y(t)\), which is given by

\begin{align} y_1 = \langle e^{j\w_o 1 t}, y(t)\rangle = \frac{1}{T}\int_{T} e^{-j\w_o t } y(t)dt.\label{eqn:calcy1} \end{align}

Of course, in practice, we only measure \(y(t)\) at discrete points, so we will approximate \eqref{eqn:calcy1} with a sum, which we will discuss further in Section No description for this link.

Let's take another look at \eqref{eqn:out1}. Pretend for a minute that these are finite dimensional vectors. To aide in our suspension of disbelief, call \(y(t)\), \(y\), and \(\Hjw\) as \(\lambda\) and \(\ejw\) as \(v\). Then \eqref{eqn:out1} looks like

\begin{equation} y = \lambda v. \end{equation}

We want to solve for \(\lambda\). We know how to do that, we just multiply both sides by \(v^T\).

\begin{align} v^Ty &= v^Tv \lambda\\ \lambda &=\frac{1}{v^Tv}v^Ty. \end{align}

Or in terms of the fancy inner product,

\begin{equation} \lambda = \frac{1}{\langle v, v \rangle} \langle v, y \rangle. \end{equation}

This is exactly what we did with the Fourier integral in \eqref{eqn:calcy1}. It just so happens that because our inner product on \(\T_T\) is normalized by \(1/T\), that \(\langle v,v\rangle = 1\), (the \(e^{j\w k t}\) are an ortho*/normal*/ basis).

Now, let's complicate things slightly. Suppose that we change the amplitude and phase of our input signal from \(\ejw\) to \(M e^{j(\w t - \phi}\), \ie, change with an amplitude and phase shift. This is just multiplication by a complex a number

\begin{equation} u(t) = M e^{j(\w t - \phi)} = (M e^{j\phi}) \ejw \end{equation}

Now, because \(\mathcal{A} = h(t) * u(t)\) is linear, we will find that

\begin{align} h(t)*u(t) & = (M e^{j\phi}) \Hjw \ejw\\ & = y(t)\\. \end{align}

Following the same argument as before, we see that

\begin{align} y(t) &= y_1\ejw\\ y_1 &= \langle\ejw, y(t)\rangle\\ & = (ME^{j\phi}) \Hjw. \end{align}

In practice, we may not know precisely what \(Me^{j\phi}\) is, especially when we are doing things in closed loop. But we do know that

\begin{align} u(t) & = \sum_{k=-\infty}^{\infty} u_k e^{j\w_okt} \\ & u_1 = Me^{j\phi}\\ & u_{k\neq 1} = 0. \end{align}

Thus \(y_1 = u_1 \Hjw\), so that

\begin{equation} \Hjw = \frac{y_1}{u_1} = \frac{ \langle\ejw, y(t)\rangle}{\langle\ejw, u(t)\rangle}. \end{equation}

Given a sequence of measured input/output data, we must approximately compute Fourier integral. The two most straightforward methods are to approximate the integral with a square rule or a trapezoid rule. I will show here that if we have a \(T\)-periodic signal sampled at regular intervals \(T_s\), and if \(T/T_s \in \mathbb{N}\), then the trapezoid approximation rule is exactly equal to a square approximation. First, considering integrating a signal, not necessarily periodic, from \(t=0\) to \(t=T\). This is illustrated in Figure 2. Assume the signal has been sampled with sampling interval \(T_s\) and that \(T/T_s\) is an integer.

Approximating the integral with a square rule requires \(N\) samples while the trapazoidal approximation requires requires \(N+1\) samples. Then the trapezoid approximation is given by

\begin{align} \int_0^T y(t) dt &\approx \sum_{k=1}^{N} \frac{T_S}{2} (y_k + y_{k-1})\nonumber\\ &= \frac{T_s}{2}\sum_{k=0}^N y_k + \frac{T_s}{2} \sum_{\ell=0}^{N}y_{\ell} - \frac{1}{2}(y_0+y_N)\nonumber\\ &= T_S\sum_{k=0}^N y_k - \frac{T_s}{2}(y_0+y_N) \label{eqn:trapzsimp} \end{align}

Now, suppose that \(y(t)\) is periodic with \(T\), so that for all \(t\), \(y(t) = y(t+T)\), and in particular, \(y(0)=y(T)\). Then \eqref{eqn:trapzsimp} becomes

\begin{align*} T_S\sum_{k=0}^N y_k - \frac{T_s}{2}(y_0+y_N) &=T_s \sum_{k=0}^{N-1} y_k + T_sy_N - \frac{T_s}{2}(y_N+y_N)\\ &= T_s\sum_{k=0}^{N-1} y_k. \end{align*}

This is precisely the sum we would have if we approximated the integral with a square rule instead of a trapezoid rule. Putting this together with the Fourier integral, we have

\begin{align} y_k &= \frac{1}{T} \int_0^T e^{-j\w k t}y(t) dt \approx \frac{T_s}{T}\sum_{\ell=0}^{N-1} e^{-j\w k T_s\ell}y_k\nonumber\\ &=\frac{1}{N}\sum_{\ell=0}^{N-1} e^{-j\w k T_s\ell}y_k.\label{eqn:dft_omega} \end{align}

Unsurprisingly, this is precisely the definition of the Discrete Fourier Transform (DFT). The form of \eqref{eqn:dft_omega} can be put into the standard form of the DFT by noting that \(\w T_s = \frac{2\pi}{NT_s}T_s\), which gives

\begin{align} y_k &=\frac{1}{N}\sum_{\ell=0}^{N-1} e^{\frac{-j 2\pi k}{N} \ell}y_k.\label{eqn:dft_normal} \end{align}

It is rather fascinating that, for an evenly sampled, strictly periodic signal y(t), the sum \eqref{eqn:dft_normal} will exactly yield the integral expression REF for the first \(\pm N/2\) coefficients.

It is pretty easy to derive the DFT strictly from a linear algebra perspective. Previously, we talked about a signal \(y(t)\). Now, we think of a vector \(y\in \mathbb{C}^N\). If you like, this is the sampled signal \(y(T_sk)\). We claim that the set of vectors

\begin{equation} \{v_k\}_{k=0}^{N-1} = \begin{bmatrix} e^{jk\w_oTs0}\\ e^{jk\w_oTs1}\\ \vdots\\ e^{jk\w_oTs(N-1)} \end{bmatrix} = \begin{bmatrix} e^{\frac{j2\pi k}{N}0}\\ e^{\frac{j2\pi k}{N}1}\\ \vdots\\ e^{\frac{j2\pi k}{N}(N-1)}\\ \end{bmatrix} \end{equation}

form a basis for \(\mathbb{C}^N\). In constrast to the infinite dimension case (the continous time Fourier Series) because we are now working in finite dimensional space, it is sufficient to show that \(v_k\perp v_{\ell},~k\neq\ell\) to show that the \(\{v_k\}\) form a basis. We proceed with direct calculation:

\begin{align} \langle v_k,~v_{\ell}\rangle &= \sum_{n=0}^{N-1} e^{\frac{-j2\pi k}{N}n} e^{\frac{j2\pi \ell}{N} n}\\ &= \sum_{n=0}^{N-1} e^{\frac{j2\pi (\ell - k)}{N}n}. \end{align}

It should be clear that if \(k=\ell\), then the exponential becomes one, so that the sum is equal \(N\). If \(k\neq\ell\), then we use the geometric sum to find that

\begin{align} \sum_{n=0}^{N-1} e^{\frac{j2\pi (\ell - k)}{N}n} = \frac{1 - e^{j2\pi (\ell - k)}} {1 - e^{\frac{j2\pi (\ell - k)}{N}}} = 0 \end{align}

Thus, the \(v_k\) form a basis for \(\mathbb{C}^N\). Given a vector \(y\in\mathbb{C}^N\), we can represent it in this basis as the linear combination of \(v_k\)

\begin{equation} y = \sum_{k=0}^{N-1} Y_kv_k, \end{equation}

where the coefficients \(Y_k\) are found the by projecting \(y\) onto each basis vector

\begin{equation} Y_k = \langle v_k,~ y\rangle. \end{equation}

In other words

\begin{equation} \begin{bmatrix} Y_0\\Y_1\\\vdots\\Y_{N-1} \end{bmatrix} = \begin{bmatrix} ---& v_0^*& ---\\ ---&v;_1^*&---\\ &\vdots&\\ ---&v;_{N-1}^*&--- \end{bmatrix} y \end{equation}

This is the "SFT" (Slow Fourier Transform) of \(y\). The FFT (Fast Fourier Transform) just does this same thing with high computational efficiency.

In general, when we hit our (supposedly) linear system \(g(\cdot)\) with an eigenvector \(\ejw\), we will not get back a perfect \(\Hjw\ejw\). Of the issues causing this are

1. Transients. The exponentials \(\ejw\) are only eigenvectors if they are applied for all time (starting infinitely long ago).
2. Noise. We will always collect a signal \(z(t) = y(t) + \eta(t)\), where \(eta(t)\) stems from some random process.
3. Non-linearity. Our system will never be perfectly linear. In the frequency response, we expect this to show up as non-trivial amounts of energy in specific higher harmonics, i.e., \(e^{j\omega n t}\).

We consider each item in turn.