Beyond Nyquist in Frequency Response Function Identification: Applied to Slow-Sampled Systems

Fast-sampled models are essential for control design, e.g., to address intersample behavior. The aim of this paper is to develop a non-parametric identification technique for fast-sampled models of systems that have relevant dynamics and actuation above the Nyquist frequency of the sensor, such as vision-in-the-loop systems. The developed method assumes smoothness of the frequency response function, which allows to disentangle aliased components through local models over multiple frequency bands. The method identifies fast-sampled models of slowly-sampled systems accurately in a single identification experiment. Finally, an experimental example demonstrates the effectiveness of the technique.


I. INTRODUCTION
Systems that have actuation and dynamics above the Nyquist frequency of the sensor, known as slow-sampled systems, are becoming increasingly common in for example vision-in-the-loop systems.As a consequence of the Nyquist-Shannon sampling theorem [1], slow-sampled systems are typically identified up to the Nyquist frequency of the slowsampled sensor.In sharp contrast, fast-sampled models of systems are typically required for control design and performance evaluation, e.g., for the use in evaluating intersample performance [2].
Non-parametric frequency-domain representations are often used for performance evaluation and controller design of linear time-invariant (LTI) systems.An example is manual loop-shaping [3] and parametric system identification [4].A common method for frequency-domain representation is through Frequency Response Functions (FRFs).FRFs can directly be identified from input-output data and are fast, accurate, and inexpensive to obtain [5], [6].Finally, FRFs allow for direct evaluation of stability, performance, and robustness [7].
The identification of fast-sampled models for slowsampled systems is challenging, since the maximum achiev-able identification frequency of traditional FRF identification for LTI systems is limited by the Nyquist frequency of the slow-sampled sensor.The key reason is that fast-sampled outputs are aliased when sampled by a slow-sensor, resulting in indistinguishable contributions in the output, and hence, a fast-sampled model cannot be uniquely recovered [8].As a result, techniques for identifying fast-sampled models for slow-sampled systems, that are required for control design and performance evaluation, are necessary.
Important developments have been made in identification techniques for slow-sampled systems, primarily in continuous-time and multirate parametric system identification.First, continuous-time system identification aims to identify a continuous-time parametric model using inputoutput data, as outlined in [9].Typically, these methods require intersample assumptions on the input signal, e.g., zero-order hold or bandlimited signals [10].Second, parametric identification of slow-sampled systems are developed and include methods for impulse response [11] and outputerror [12] model estimation.Lifting techniques, such as using subspace [13], frequency-domain [14] or hierarchical identification techniques [15], are also developed.These methods focus on parametric identification, require intersample assumptions on the input signal and do not exploit fastsampled inputs, and consequently, do not disentangle aliased components.
Although methods for identification beyond the Nyquist frequency of slow-sampled systems have been developed, an efficient and systematic methodology for single-experiment FRF identification of fast-sampled models, that disentangles aliased components with arbitrary input signals, is currently lacking.In this paper, slow-sampled systems are identified with excitation signals that cover the full frequency spectrum, where aliased components are disentangled from each other through exploiting the assumption of smooth behavior in the frequency response of a system.Generally, this assumption is at the basis of modern FRF identification, as seen in techniques for LTI single-rate systems, such as Local Polynomial Modeling (LPM) [5] or local rational modeling [16].In fact, LPM for LTI single-rate systems is recovered as a special case of the developed framework.The key contributions of this paper include the following.C1 Identification of non-parametric fast-sampled models for slow-sampled systems, by exciting the full frequency spectrum and aliased components are disentangled from each other by assuming smooth behavior in the frequency domain (Section IV).C2 Validation of the framework for identification of slow- sampled systems in an experimental setup (Section V).
Notation: Fast-sampled signals are denoted by subscript h and slow-sampled signals by subscript l.The N -points and M -points Discrete Fourier Transform (DFT) for respectively finite-time fast-sampled and slow-sampled signals are given by with frequency bin k, sampling times T h and T l , discretetime indices for fast-sampled signals with integers Z and N, M the amount of data points of the fast and slow sampled signals, and generalized frequency variable The sampling times of the slow-sampled and fast-sampled signals relate as T l = F T h , with downsampling factor F ∈ Z >0 .The complex conjugate of A is denoted as A and the complex conjugate transpose as A H .The complement of sets is given by A \ B. The expected value of a random variable X is given by E {X}.

II. PROBLEM FORMULATION
In this section, the problem that is considered in this paper is presented.First, the identification setting is presented.Finally, the problem addressed in this paper is defined.

A. Identification Setting
The goal is to identify a fast-sampled non-parametric model Ĝ with sampling rate f h = 1 T h , using the slowsampled output y l with sampling rate , where the open-loop control structure considered is visualized in Fig. 1.The system G and the noise model H are LTI Single-Input Single-Output (SISO) systems.The input-output behavior of LTI SISO systems in the frequency domain is given by with system transient term T G (Ω k ) [5].This reveals that a single frequency of Y l is influenced by a single frequency of U h , also called the frequency-separation principle.The measured slow-sampled output is a downsampled version of the fast-sampled output as shown in Fig. 1, i.e., , where E(k) is filtered zeromean white noise, and is assumed to be independent and identically distributed.In time-domain, the downsampling operation in (4) equates to S d y h (n) = y h (nF ).By applying the downsampling operation in the frequency-domain in (4), the DFT of the slow-sampled output is given by [17] Y By substituting the input-output behavior of the fast-sampled system G(Ω k ) in ( 3), the slow-sampled output is given by

B. Problem Definition and Approach
The DFT Y l in ( 6), for a single frequency bin k, is influenced by F frequencies of G and U h .This is caused by aliasing due to the downsampling operation.Hence, the fastrate system G(Ω k ) can in general not be uniquely identified with the slow-sampled output Y l for arbitrary inputs U h .
The problem considered in this paper is as follows.Given fast-sampled input data u h and slow-sampled output signal y l with DFTs U h and Y l shown in (4), identify a fast-sampled model of G(Ω k ) in the frequency-domain for bins k ∈ Z [0,N −1] , i.e., up until the fast-sampled sampling frequency f h , for the identification setup seen in Fig. 1.The approach is developed in two steps.
1) Development of an intuitive idea in Section III for identifying slow-sampled systems using a dedicated input signal for a sparse frequency grid.2) Development of the full approach in Section IV using arbitrary input signals and full frequency grids, leading to contribution C1.

III. INTUITIVE IDEA: IDENTIFICATION WITH SPARSE
FREQUENCY SPECTRUM The first step in Section II-B is developed, where aliasing is precisely traced for each input signal such that the slowsampled output Y l (k) is only influenced by a single fastsampled input U h (k).This step is intended for conveying the intuitive idea, leading to the full approach, i.e., the second step in Section II-B, in Section IV.
The transient contribution T G = 0 in this section, meaning that the transient is neglected to facilitate the development of the intuitive idea, in that case ( 6) is equal to Illustration of problem (top) and intuitive idea (bottom) for identifying a fast-sampled model of G(Ω k ) ( ) with slow-sampled output Y l .The Nyquist frequency and multiples when sampling the output a factor F = 3 slower is shown as ( ) .Top: Non-zero input U (k + iM ) ∀i ∈ {0, 1, 2} and associated gain |G(Ω k+iM )| ( ) show that the slow-sampled output Y l (k) ( ) is a summation as in (7), resulting in F = 3 unknowns, but only 1 equation.Bottom: Non-zero input U (k + i(M + 1)) ∀i ∈ {0, 1, 2} and associated gain |G(Ω k+i(M +1) )| ( ) result in a single contribution of the summation in (7) influencing Y l (k + i) ( ) by deliberately not exciting specific bins ( ).Hence, the fast-sampled system G(Ω k+i(M +1) ) can be uniquely recovered at frequency bins k + i(M + 1).
showing that the slow-sampled output is a summation of the baseband response for f = 0, and aliased components for . Due to the summation of the baseband and aliased contributions, the fast-sampled system can in general not be uniquely recovered from slow-sampled output Y l .An example is seen in the top of Fig. 2.
The key idea is that by designing the input signal U h such that the output at frequency bin k is only influenced by a single contribution of the summation G (Ω k+M f ) U h (k + M f ) in (7), the fast-sampled system can be recovered for a subset of all frequency bins.This means that input signals of the form result in the input-output behavior which shows that the summation disappears, and hence, the fast-sampled system G can be uniquely recovered for the frequency bins k +iM .The sparse set of excited frequencies is from now on denoted by S, i.e., the signal in (8) can be represented by U h (k) ̸ = 0 ∀k ∈ S and U h (k) = 0 ∀k / ∈ S. The general concept of identifying slow-sampled systems with sparse multisines is shown in the bottom of Fig. 2. By noting that due to the M -periodicity of the DFT, ( 9) is rewritten as An estimate of the system G is given by Ĝ given that the excitation signal is designed such that the output Y l (k) is only influenced by a single frequency of U h (k) according to (8).An example of a sparse multisine is shown in Example 1.
Example 1 An example of a sparse multisine, that achieves broadband excitation, is given by Essentially, the sparse multisines in Example 1 avoid interference of different frequency bands, and by appropriately selecting the inputs a system estimate Ĝ is obtained at a sparse set of frequencies in each frequency band, including the bands beyond the Nyquist frequency.
The resulting estimate Ĝ(Ω k ) is only obtained at the sparse set of excited bins k ∈ S. By ranging i over multiple sets of experiments in {0, 1, . . .F − 1}, the system Ĝ(Ω k ) can be uniquely recovered for the full frequency spectrum.As a result, F experiments are necessary to identify the system, which leads to a time-intensive procedure.In the next section, a time-efficient single-experiment identification approach is developed.

IV. IDENTIFICATION WITH FULL EXCITATION SPECTRUM
In this section, the approach is developed to identify a fast-sampled model of G(Ω k ) for all frequency bins k ∈ Z [0,N −1] , given slow-sampled outputs, therewith constituting contribution C1.This is realized by exciting the full frequency spectrum, where aliased contributions are disentangled by exploiting a smoothness condition on G and the transient.In contrast to Section III, where for each frequency bin k the F unknowns in (7) were reduced to a single unknown as seen in (10), in this section the amount of equations are increased to F for each frequency bin k.First, the method is presented.Second, a covariance analysis is provided.Finally, the developed approach is summarized in a procedure.

A. Identification of Slow-Sampled Systems with Full Excitation Spectrum
The frequency response of the system G(Ω k ) and transient T G (Ω k ) are assumed to be smooth, as is formalized in Assumption 1 and Assumption 2. The smoothness assumption enables disentangling aliased components, and consequently, the need for a sparse excitation signal in Section III is relaxed.

Assumption 1
The frequency response of the fast-sampled system G(Ω k ) can be approximated in a local window r ∈ Z [−nw,nw] , with 2n w + 1 the window size, as an R th order polynomial as Assumption 2 The summation of transients, seen in (6), is assumed to be smooth in the local window r ∈ Z [−nw,nw] , i.e., The assumption of a locally smooth system and transient is commonly imposed and at the basis of modern FRF identification, and is valid since G(Ω k ) and T G (Ω m ) are functions with continuous derivatives up to any order [5], [16].Hence, the assumption of a smooth summation of transients is equally reasonable, since the individual transient contributions are smooth.By substituting the polynomial models for G(Ω k ) in ( 13) and for the transient in ( 14), the slow-sampled output in ( 6) is rewritten with parameter vector Θ and data vector K for the local window r as The parameter vector Θ(k) ∈ C 1×(R+1)(F +1) is given by and data vector K(k + r) ∈ C (R+1)(F +1)×1 is given by with input vector where where Y l,nw ∈ C 1×2nw+1 , K nw ∈ C (R+1)(F +1)×2nw+1 and V nw ∈ C 1×2nw+1 are constructed as The fast-sampled system G(Ω k ) and transient T (Ω k ) are uniquely identifiable for all k ∈ Z [0,N −1] , in the presence of aliasing and with a full excitation spectrum, as in Theorem 1.
designed such that K nw is of full row rank, an estimate of the parameter vector Θ in the least-squares sense of (20) is uniquely determined as and the fast-sampled model of system 14) can be obtained.
Proof: The matrix inverse K nw K H nw −1 uniquely exists if K nw K H nw has full rank, which is achieved if the rank of K nw is equal to the row rank of K nw and K nw is full row rank.The interpretation of Theorem 1 is as follows.First, sufficient data 2n w + 1 should be available, such that (22) leads to a unique solution of the (R + 1)(F + 1) parameters, hence 2n w + 1 ≥ (F + 1)(R + 1), which is satisfied by design of wide matrix K nw .In other words, by applying the smoothness assumption, the estimation of the (R+1)(F +1) parameters Θ in (22) utilizes 2n w + 1 outputs in Y l (k + r).This explains how the smoothness assumption allows to disentangle F aliased contributions at a frequency bin k.Additionally, no overlapping between windows k + r + iM ∀i ∈ Z [0,F −1] is allowed, otherwise K nw is not full row rank, and hence, M + 1 > 2n w + 1.Second, the system in (22) can be solved uniquely if K nw is full row rank.As a consequence, aliased and transient contributions can be disentangled if all inputs in the local window and at the aliased windows are sufficiently 'rough', that is formalized for a single local window in [18].For Theorem 1, this means that the spectral difference . This condition is fulfilled by, e.g., random-phase multisines [5].The developed framework is illustrated in Fig. 3.
Remark 1 Note that traditional LPM for single-rate LTI systems is recovered as a special case of the developed framework by setting F = 1.

Remark 2
The identified FRF by sparse multisines from Section III can be interpolated at the non-excited frequency bins k / ∈ S, similar to [19], by seeing them as a special case of the framework in this section.The condition on U (k + r) in Theorem 1 is inherently satisfied by sparse excitation.

B. Variance Estimate of the FRF
The developed framework enables estimation of the variance of the identified FRF.The variance of Ĝ, that is estimated using (22), is given by Theorem 2.

Theorem 2
The estimated variance of the FRF Ĝ that is estimated using (23) is given by that is an estimate of the true variance of the identified FRF with C V the variance of the noise and an estimate based on measurements ĈV , noise interpolation error O int H [5], and Proof: The proof extends [5,Appendix 7.E] to F frequency bands.In particular, by combining (20) and ( 22) into (Y nw − ΘK nw ) S = V nw S, Θ 1 0 the factor F appears in the difference between the true and estimated system G and Ĝ.The variance of the noise is equal to An estimate of the noise variance is calculated by taking an average over the local window, see [5,Appendix 7.B] for technical details, i.e., with Vnw = Y l,nw − Θ(k)K nw .

C. Developed Procedure
The developed approach is summarized in Procedure 1, that links the main results in this paper.

V. VALIDATION
In this section, the developed method is validated on an experimental setup, leading to contribution C2.
Procedure 1 (Identifying slow-sampled systems with full frequency spectrum) 1) Construct u h such that it satisfies the requirements of U given in Theorem 1. 2) Apply input u h to system and record the output y l .
3) Take the DFT of input u h and output y l using (1).

A. Measurement Setup
The experimental setup is shown in Fig. 4. The setup consists of two rotating masses connected via a rubber band.The rotating masses are each actuated by a DC motor.The input u h and output y l are respectively the torque and rotation of the first mass.The second mass is virtually suspended to the fixed world by a feedback controller.The excitation signal for the developed approach that excites the full frequency spectrum is a random-phase multisine and has an root-mean-square value of 8.3•10 −3 Nm.For comparison purposes, the intuitive idea from Section III uses sparse multisines with a root-mean-square value of 9.6 • 10 −3 Nm and are designed as in Example 1.A photograph and a schematic overview of the test setup are seen in Fig. 4. The settings used during identification are shown in Table I.

B. Experimental Results
The fast-sampled system is identified using the developed approach from Section IV.For comparison purposes, the sparse multisine approach from Section III, a traditional approach and using the fast-sampled output y, that is not available for the other approaches, are used to identify the fast-sampled system.The results for the developed, sparse multisine and traditional approach are seen in Fig. 5, Fig. 6 and Fig. 7.The following observations are made 1) From Fig. 5 and Fig. 6, it is observed that the developed approach and the sparse multisines are capable of identifying the dynamics above the Nyquist frequency of the sensor.The developed approach of Section IV that assumes smooth behavior in the frequency-domain has significantly lower variance and a factor F higher frequency resolution.2) From Fig. 7, it is observed that exciting the system with the full frequency spectrum and performing FRF identification by Ĝ(Ω k ) = Y l (k)/U h (k) cannot accurately identify the fast-sampled system G, due to aliasing.Additionally, the estimated FRF is dominated by a periodic behavior, which is explained because the slow-sampled DFT Y l (k), that is M -periodic, is used for all frequencies and aliasing is not accounted for.3) From Fig. 5, the variance of the developed approach that excites the full frequency spectrum repeats every f l .For example, the increase of variance around 23 Hz repeats at 53 Hz.This is explained because the variance from ( 24) is determined with the estimated noise variance Ĉv from (29), that uses the slow-sampled data Y l , and hence, is periodic in f l .Additionally, the mirroring effect that is observed, e.g., the increase of variance at 23 Hz is mirrored to 7 and 37 Hz, is caused because the DFT of Ĉv is symmetrical in 1 2 f l .

VI. CONCLUSIONS
The results in this paper enable identifying FRFs of slowsampled systems where aliasing occurs.The key step is assuming smooth behavior of the system FRF, which allows to appropriately disentangle aliased contributions when exciting the full frequency spectrum.Furthermore, covariance estimates of the FRF are provided.Finally, the framework is validated through experimental results.The dual case, where outputs are fast-sampled and inputs slow-sampled, is trivial by assuming appropriate interpolator behavior.The developed approach is a key enabler for closed-loop, multivariable and parametric system identification and control design for slow-sampled systems, such as vision-in-the-loop systems.

Fig. 4 .
Fig. 4. Left: Picture of experimental setup used.Right: Schematic overview of experimental setup.

Fig. 5 .Fig. 6 .
Fig.5.Identified FRF Ĝ(Ω k ) for excitation by random-phase multisines covering the full frequency spectrum from Section IV ( ) with covariance estimate from (24) ( ) (right axis).The identified FRF based on fastsampled data is shown as ( ) and multiples of the Nyquist frequency of the slow sensor as ( ) .

Fig. 7 .
Fig.7.Identified FRF Ĝ(Ω k for excitation by full frequency spectrum random-phase multisines by performing Ĝ(Ω k ) = Y l (k)/U h (k) using the same dataset as the full approach ( ) .The identified FRF based on fastsampled data is shown as ( ) and multiples of the Nyquist frequency of the slow sensor as ( ) .
(15)notes the Kronecker product.Collecting the column vectors from(15)in matrices for the window r ∈ Z [−nw,nw] gives

TABLE I EXPERIMENTAL
SETTINGS.