By Scott Hirsch, MathWorks and John Glass, MathWorks

The System Identification Toolbox enables you to use measured input-output data to build and evaluate linear models of dynamic systems. You can use this toolbox to develop models of a black-box system without having to fully characterize the mathematics governing the system behavior.

Version 6.0 of the System Identification Toolbox introduced important new options for developing models. The toolbox now supports estimation and validation using frequency domain data. This release also added support for the creation of low-order continuous-time process models. This article explores the practical application of these new features, using a simple linear circuit as the device under test. We also demonstrate data collection using the Data Acquisition Toolbox and Instrument Control Toolbox.

For the purpose of this article, we will use a known device under test, in this case a low-pass filter. The first-order filter is implemented as an RC circuit, as shown in Figure 1. The low-pass filter is one of the simplest linear time-invariant systems. It can be characterized mathematically in both the frequency and time domains. This makes it an ideal demonstration platform for understanding the system identification procedure, as we can validate our results with theoretical predictions.

The transfer function of this circuit, with approximate resistance \(R = 1 k\Omega\) and capacitance \(C=0.1 \mu F\), is:

\[T(j\omega) = \frac{1}{1+j\omega C R}\]

The cutoff frequency (the frequency at which the input signal is attenuated by 3 dB) is

\[\omega_o = \frac{1}{CR} = 10,000 \text{ Rad/s} = 1,592 \text{ Hz}\]

Our goal is to use the System Identification Toolbox to develop a linear time-invariant model of the circuit.

The most crucial step in system identification is to collect an appropriate data set for estimation. Using the Data Acquisition Toolbox and the Instrument Control Toolbox, we can collect this data directly into MATLAB^{®}. This significantly reduces the time for iterating as we experiment with different data sets (eliminating the need to work with multiple tools for data collection and modeling) and also lets us build automated test applications that incorporate our design (modeling) tasks.

This test uses an Agilent 33120A Function Generator to excite the device under test. The Instrument Control Toolbox controls this function generator via GPIB (General Purpose Interface Bus). The system excitation and response are captured with the line input of a Windows sound card. The Data Acquisition Toolbox is used to read (stream) the data into MATLAB so that we can configure and acquire data from the sound card.

This section details the procedure for using the Instrument Control Toolbox and the Data Acquisition Toolbox to collect data for system identification. We start by configuring the sound card to acquire 1.0 seconds of data on two channels at 96 kHz.

We start by creating an analog input object (using our Windows sound card), and adding two channels to the object.

```
ai = analoginput('winsound');
addchannel(ai,1:2);
```

By default, the sound card uses only a small set of specific sample rates (8.000 kHz, 11.025 kHz, 22.050 kHz, and 44.100 kHz). We turn off StandardSampleRates to gain access to 96 kHz.

set(ai,'StandardSampleRates','Off') Fs = 96000; set(ai,'SampleRate', Fs)

We define a manual trigger to let us start the card running before we are ready to acquire data and to eliminate start-up transients from our acquired data. We configure the trigger command to acquire 96,000 samples.

set(ai,'TriggerType','manual'); set(ai,'SamplesPerTrigger', Fs)

We start the sound card:

start(ai)

We now configure the function generator to create our system excitation. Using the Test & Measurement Tool, a new graphical user interface provided with the Instrument Control Toolbox, we can quickly connect to and configure our hardware.

```
tmtool % Open TMTOOL
```

Selecting **Scan For All Hardware** from the **Tools** menu finds all of the instruments connected to our computer. Navigating through the hardware node, we find the Agilent 33120A function generator connected to the GPIB board. We select this device to display a communication panel on the right pane of TMTOOL, and then click **Connect** to connect MATLAB to the function generator. This creates an instrument interface object, which is a MATLAB variable that handles the communication between MATLAB and this particular GPIB device. The current state of TMTOOL is shown in Figure 2.

Version 2.0 of the Instrument Control Toolbox introduced support for industry standard instrument drivers via MATLAB device objects, which let us communicate with our instrument without knowing the specific underlying commands. When you right-click on the Instrument Drivers node and choose **Scan for Instrument Drivers**, you find MATLAB instrument drivers and IVI and VXIplug&play drivers on the computer. The MATLAB instrument driver has already been created for the 33120A, so it appears under the MATLAB Instrument Drivers node. Right-click this driver and select **Create Device Object Using Driver**. This action creates a device object to work with the 33120A in MATLAB, utilizing the interface object we created above. Select this device object in the Instrument Objects node to gain point-and-click access to all of the instrument’s functionality, which has been programmed into the driver.

We use the **Configure** tab to set the function generator to output 2V RMS random noise (Figure 3).

The Session Log shows all of the MATLAB code used for this instrument control session (Figure 4). We can copy commands from this window or click **Save Session** to generate an M-file to use again later.

Now that the analog input object is ready, and the instrument is configured appropriately, we can collect our data by triggering the acquisition and bringing the data into MATLAB.

trigger(ai) % Trigger the acquisition [data,time] = getdata(ai); % Bring the data into MATLAB

We collected the first data set with white noise system excitation. Because white noise input simultaneously excites information over all frequencies below the data collection Nyquist rate, it is particularly suited to frequency domain modeling. There are two benefits to using frequency domain data for modeling linear systems such as this filter:

- The frequency response function is computed before any model is estimated. This gives a good indication of where the break frequency is and serves as a check when you initially examine the experimental data.
- The overall size of the data set used during system identification is reduced, as the frequency response function is typically estimated by averaging smaller chunks of a large dataset.

However, frequency response function estimation requires an extra step in the model-generation process and can heavily depend on the identification and data collection parameters.

There are frequency response function estimation algorithms available in both the System Identification Toolbox and the Signal Processing Toolbox. For this article, we used the frequency response estimation functions in the System Identification Toolbox.

We applied two simple guidelines when designing the input signal and acquisition system to estimate the frequency response function:

- Choose the input amplitude so that the output response is a minimum of 10x larger then the noise in the system.
- Choose an acquisition sample time on the order of 10x the system bandwidth. This provides a good balance between model fidelity and minimizing the size of the identification data set.

Figure 5 shows the measured system input and response. The mean of the data has been removed by eliminating any offset. As expected for white noise passed through a low-pass filter, the response is of lower amplitude than the input.

Now that the data has been collected we are ready to estimate a model. First, we store the data in a data object to be used throughout the identification process.

Create an `iddata`

system identification data object.

z = iddata(Output,Input,Ts);

The data can now be detrended by removing the mean values or linear trends from the signals:

zd = detrend(z);

For frequency domain identification, the next step is to compute the frequency response function for this data. We use `spafdr`

from the System Identification Toolbox, which estimates frequency response and spectrum by spectral analysis with frequency dependent resolution. The last input argument specifies that the frequency response should be computed over 100 logarithmically spaced frequencies from 50 to 60,000 Hz.

```
zd_fdr = spafdr(zd,[],{50,60000,100});
bode(zd_fdr); % Plot the bode response
```

We chose a range where the frequency response matches the first-order behavior that we expected. Upon examination of the computed frequency response, it appears that our estimate of the frequency response function is valid only over the range 620 and 60000 rad/s. Outside of this range it does not appear that the data has converged to an accurate frequency response estimate.

At this point we can perform an estimation of the data. First, we select the range with which we plan on working. We use the System Identification Toolbox `idfilt`

function to filter the frequency response data with an ideal frequency domain filter:

zd_fdr_f = idfilt(zd_fdr,[620,60000]);

The System Identification Toolbox offers many model structures to identify discrete and continuous models. We want to identify a first-order model, so we use the `idproc`

model structure. This model structure represents simple continuous-time process models, characterized by static gain, possible dead time, and dominating time constant(s). The structure allows the specification of the order of the desired model. Examples of this type of model structure include the following transfer functions:

\[P(S) = \frac{K}{\left(T_{p_1 s} + 1\right)\left(T_{p_2 s} + 1\right)}, \quad P(S) = \frac{Ke^{-sT_d}}{T_{p_1 s} + 1}, \quad P(S) = \frac{K}{(T_{ws})^2 + 2 \zeta T_{ws} + 1}\]

For the identification of the filter we used a model of the form

\[P(S) = \frac{K}{T_{p_1 s} + 1}\]

We identify this model with the following code:

P1 = idproc('P1'); % Create a first-order idproc object P1.Tp1.min = 1e-8; % Set the lower bound for the pole location P1_wn = pem(zd_fdr_f,P1); % Identify the model using the frequency

Information about the model is provided in the display summary of `P1_wn`

:

P1_wn

` Process model with transfer function`

*with K = 0.95488
Tp1 = 8.8552e-005
Estimated using PEM from data set zd_fdr_f
Loss function 0.000828687 and FPE 0.000845428*

Relating the transfer function to the equations governing our system, we see that the cutoff frequency of the model is simply the inverse of `Tp1`

:

omega_c_wn = 1/P1_wn.Tp1.value; fc_wn = omega_c_wn/(2*pi)

` fc_wn = 1.7973e+003`

The percent error in cutoff frequency between the model and the measured data is:

fc_err = abs(((fc_wn - fc_mw)/fc_mw)*100)

` fc_err = 11.2025`

We find that this model has an error of greater than 10% in cutoff frequency. This is a first approximation, where the limiting factor in this case is the limited frequency range that was identified. There are two approaches to improving the frequency response estimation, (1) Collect more white noise data or (2) use a different input signal. In many applications it may not be reasonable to expect that a larger data set can be acquired. Therefore, it is worth looking at a different input signal.

Instead of white noise, we use a chirp signal (swept sine). The chirp signal is particularly useful because the frequencies that are put through the system are directly specified, while using a limited amount of data.

We return to tmtool to define a chirp signal that sweeps logarithmically from 50 Hz to 10,000 Hz in 1.5 seconds as our system excitation.

Figure 7 shows the measured system input and response. We clearly see the effect of the low-pass filter as the frequency of the input signal increases.

Following the same procedure as before, we compute and display the frequency response of the system.

zsw = iddata(Output,Input,Ts); zswd = detrend(zsw); zf_sw = spafdr(zswd,[],{50,60000,100}); bode(zf_sw);

Notice that the frequency response is much cleaner for this data set, so we do not need to apply a filter.

Creating a model with this new data set yields much better results:

```
P1 = idproc('P1');
P1.Tp1.min = 1e-8;
P1sw = pem(zf_sw,P1)
```

` Process model with transfer function`

IMAGE

with K = 0.99483

Tp1 = 9.905e-005

Estimated using PEM from data set zf_sw

Loss function 0.000851581 and FPE 0.000868785

omega_c_ws = 1/P1sw.Tp1.value; fc_ws = omega_c_ws / (2*pi)

` fc_ws = 1.6068e+003`

fc_err = abs(((fc_ws - fc_mw)/fc_mw)*100)

` fc_err = 0.5834`

We find that the chirp system excitation yields a more precise model. We can also compare the time-domain response of this model with the measured data (Figure 9):

compare(zswd,P1sw)

This article discussed applied frequency domain system identification for developing a low-order continuous-time model of a low-pass filter. We found that for this particular system, exciting the system with a chirp signal yielded the most accurate model with a small set of collected data. We also used the Data Acquisition Toolbox and Instrument Control Toolbox to easily capture data directly into our design environment. Combining test and measurement capabilities with our system identification tools accelerates the iterative process of developing system models with measured data.

Published 2004