# Second Order Discrete-Time System Demonstration

Discrete-time systems are remarkable: the time response can be computed from mere difference equations, and the coefficients a_{i}, b_{i} of these equations are also the coefficients of H(z). Here, I try to illustrate this remarkableness by converting a continuous-time second-order system to an approximately equivalent discrete-time system. With a discrete-time model, we can then easily compute the time response to any input. But note that the goal here is as much to understand the discrete-time model as it is to find the response.

Consider a second-order continuous-time system [1,2] with the following transfer function:

$$H(s)=\frac{\omega_n^2}{s^2+2\zeta\omega_ns+\omega_n^2} \qquad(1)$$

ζ (zeta) is called the damping ratio and ω_{n} is called the natural frequency. These two parameters completely define H(s). The system has two finite poles and two zeros at ω = ∞.

The poles of H(s) are $-\zeta\omega_n \pm \omega_n\sqrt{\zeta^2-1}$ . For 0 < ζ < 1, the system is said to be underdamped, and the poles are complex-conjugate. For ζ > 1, the system is overdamped, and the poles are real. For ζ = 1, the system is critically damped, and there are two real poles, both at -ω_{n}. Figure 1a shows the s-plane of an underdamped system with complex-conjugate poles, and Figure 1b shows an overdamped system with real poles. When we convert these systems to discrete-time ones, we’ll find that, for sample frequency f_{s} = 100 Hz, the z-plane poles are as shown in Figures 1c and 1d.

In Section 1, we’ll develop a Matlab function that converts H(s) to a discrete-time system H(z), and plots the system’s poles, impulse response, and frequency response. Then we’ll show several example system responses. We’ll see that having a discrete-time model makes it easy to find the time response to an arbitrary input. In section 2, we’ll compare the time and frequency responses of the discrete-time system to those of the continuous-time system.

Figure 1.Second Order System Poles, ω_{n} = 2π*12 rad/s.

a. s-domain, ζ = 0.5. b. s-domain, ζ = 1.5

c. z-domain, ζ = 0.5, f_{s}= 100 Hz. d. z-domain, ζ = 1.5, f_{s} = 100 Hz.

## 1. Finding the discrete-time system transfer function

A second-order discrete-time system’s transfer function can be written as [3]:

$$H(z)= K\frac{b_0+b_1z^{-1}+b_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}} \qquad(2)$$

To find the coefficients of H(z), we’ll start by finding the pole locations in the z-plane. The continuous time system’s poles are:

$$p_{a1}=-\zeta\omega_n+\omega_n\sqrt{\zeta^2-1}$$

$$p_{a2}=-\zeta\omega_n-\omega_n\sqrt{\zeta^2-1}$$

As already mentioned, for ζ < 1, the poles are complex-conjugate. Using the definition of z: z= exp(sT_{s}), the z-plane poles are:

$$p_1= exp(p_{a1}T_s) \qquad \quad$$

$$p_2= exp(p_{a2}T_s) \qquad(3)$$

The denominator of H(z) is then:

$$ den=(1-p_1z^{-1})(1-p_2z^{-1}) $$

$$ =1-(p_1+p_2)z^{-1}+p_1p_2z^{-2} $$

Now consider the numerator of H(z). As mentioned, the continuous time system has two zeros at ω = ∞. For the discrete-time system, we have the option of using two, one, or no zeros. Unfortunately, we cannot place zeros at ω = ∞. The location of the zeros will affect how similar the discrete time system’s time and frequency responses are to those of the continuous time system (We’ll compare the results for two different zero locations in Section 2). If we place a zero at z = -1 + j0, this falls on the unit circle at e^{j}^{π}. You may recall that this point corresponds to f = f_{s}/2: exp(jωT_{s}) = exp(j2πf_{s}/2 * Ts) = e^{j}^{π}. For our initial examples, we’ll place two zeros at z = -1 + j0. The numerator of Equation 2 becomes:

$$ num=(1+z^{-1})(1+z^{-1}) $$

$$ = 1+2z^{-1}+z^{-2} $$

So now we can list the coefficients of Equation 2:

[b_{0} b_{1} b_{2}] = [1 2 1]

[a_{0} a_{1} a_{2}] = [1 -(p_{1} + p_{2}) p_{1}p_{2}],

where p_{1} and p_{2} are the z-plane poles from Equation 3. Finally, we note that H(jω) has gain of 1 at ω = 0. So we will choose K in Equation 2 for a gain of 1 at f = 0. At f = 0, z = e^{j0} = 1.Thus:

H(z=1) = 1 = K*sum(b)/sum(a), and

K = sum(a)/sum(b)

We now have all the equations needed to compute the coefficients of H(z), given ζ, ω_{n}, and f_{s}. The Matlab function so_demo.m is listed in the appendix of the pdf document. Besides computing the coefficients of H(z), it plots the poles of H(s), poles and zeros of H(z), and the impulse and frequency responses of H(z).

Here is a summary of the steps we used to convert H(s) to H(z), where H(z) has the form shown in Equation 2:

- Compute the poles p
_{a1}and p_{a2}of H(s). - Convert these poles to z-plane poles p
_{1}and p_{2}. - Find H(z) denominator coefficients
*a*from p_{1}and p_{2}. - Assign two z-plane zeros at z = -1 + j0, giving H(z) numerator coefficients b = [1 2 1].
- Compute K to provide gain of 1 at f = 0.

.

### Example 1 – underdamped system

We’ll use so_demo.m to create a system with ζ = 0.2, f_{n} = 7 Hz, and f_{s} = 100 Hz. The Matlab code is as follows:

zeta= 0.2; % damping ratio fn= 7; % Hz natural frequency fs= 100; % Hz sample frequency nz= 2; % number of zeros [K,b,a]= so_demo(zeta,fn,fs,nz)

The function calculates K, b, and a as:

K = 0.0436 b = 1 2 1 a = 1.0000 -1.6641 0.8387

The graphic results are plotted by so_demo.m, as shown in Figure 2. As math dictates, the poles in the left half of the s-plane are transformed to fall inside the unit circle of the z-plane. There are nz = 2 z-domain zeros at z = -1 + j0. The impulse response is oscillatory with decaying amplitude – an underdamped response. Note that if nz were equal to 1, so_demo.m would assign a single zero at z = 0 +j0 (See Section 2).

Figure 2. Second-order system, ζ = 0.2, f_{n} = 7 Hz, and f_{s} = 100 Hz.

a. s-domain poles. b. z-domain poles and zeros

c. Impulse response. d. Frequency response magnitude.

*Step response*

We can find the response to any input simply by using the system’s difference equations, computed by the Matlab function *filter* [4]. The Matlab code to compute the step response is:

N= 64; x= ones(1,N); % unit step y= filter(K*b,a,x); % step response

The step response is shown in Figure 3.

Figure 3. Second-order system step response, ζ = 0.2, f_{n} = 7 Hz, and f_{s} = 100 Hz.

*Response to a slew-rate limited step*

Here, the input a signal is a slew-rate limited step that reaches its final value after 8 samples. Here is the code to generate the input and compute the time response:

N= 64; u= ones(1,N); v= ones(1,8)/8; x= conv(u,v); % slew-rate limited pulse y= filter(K*b,a,x); % response

The input signal and time response are shown in Figure 4.

Figure 4. Top: Slew-rate limited step input.

Bottom: Second-order system output, ζ = 0.2, f_{n} = 7 Hz, and f_{s} = 100 Hz.

### Example 2 – overdamped system

This example again uses f_{n} = 7 Hz, and f_{s} = 100 Hz, but ζ = 1.5, resulting in an overdamped system with real poles. See Figure 5. The impulse response has no oscillations.

Figure 5. Second-order system, ζ = 1.5, f_{n} = 7 Hz, and f_{s} = 100 Hz.

a. s-domain poles. b. z-domain poles and zeros

c. Impulse response. d. Frequency response magnitude.

### Example 3 – sinusoidal impulse response

Here we let f_{n} = 7 Hz, and f_{s} = 100 Hz, with ζ = 0. This is an unstable system with s-plane poles on the jω axis and z-plane poles on the unit circle. Its impulse response is a sinewave at f_{n} (Figure 6).

Figure 6. Second-order system, ζ = 0, f_{n} = 7 Hz, and f_{s} = 100 Hz.

a. s-domain poles. b. z-domain poles and zeros.

c. Impulse response. d. Frequency response.

### Example 4 – Butterworth lowpass filter

Here, f_{n} = 10 Hz, f_{s} = 100 Hz, and ζ = 1/sqrt(2) ~= 0.707.This results in s-plane poles at ω_{n}*(0.707 +/- 0.707). These are the poles of a second-order Butterworth lowpass filter with -3 dB frequency of f_{n} Hz. The discrete-time version has a -3 dB point that is slightly less than 10 Hz due to the zeros at f_{s}/2 (Figure 7). A practical approach for designing Butterworth filters of any order uses the bilinear transform [5].

Figure 7. Second-order system, ζ = 0.707, f_{n} = 10 Hz, and f_{s} = 100 Hz (Butterworth lowpass).

a. s-domain poles. b. z-domain poles and zeros.

c. Impulse response. d. Frequency response magnitude.

## 2. Comparing discrete-time and continuous-time responses

Let’s consider an underdamped (ζ < 1) second-order system with H(s) given by Equation 1. The impulse response is [2]:

$$h(t)=\frac{\omega_n}{\sqrt{1-\zeta^2}}exp(-\zeta\omega_nt) \cdot sin(\omega_n\sqrt{1-\zeta^2}\,t)\,u(t)\qquad(\zeta<1)$$

We’ll compare the impulse responses and frequency responses of the discrete and continuous-time systems for two different choices for zeros of H(z): two zeros at z = -1 + j0, or a single zero at z = 0 + j0. For the single-zero case, zero location and K were computed using the impulse-invariance method (see Appendix). Reference [6] provides a second-order system example of the impulse-invariance method. Note that for both cases that follow, we had to multiply the discrete-time impulse response by f_{s} to match the amplitude of the continuous-time impulse response.

*Two zeros at z = -1 + j0*

Figure 8 compares the impulse responses and frequency responses of the discrete-time and continuous-time systems for ζ = 0.707 (Butterworth response), f_{n}= 10 Hz, and f_{s} = 100 Hz. You can see some mismatch of the impulse responses. The frequency responses diverge at higher frequencies due to the discrete system’s zeros at f_{s}/2 = 50 Hz. Accuracy of the discrete-time responses can be improved by increasing the sample rate: here, for example, we could double f_{s} to 200 Hz.

Figure 8. Discrete-time (red) and continuous-time (blue) responses, two zeros at z = -1 + j0,

ζ = 0.707, f_{n}= 10 Hz, and f_{s} = 100 Hz. a. Impulse response. b. Freq. response magnitude

*One zero at z = 0 + j0 (Impulse-Invariance Method)*

This example uses the same continuous-time system as the previous example. Using impulse-invariance, the discrete-time impulse response matches the continuous-time impulse response exactly at the sample instants (Figure 9a). The same cannot be said for the frequency response: the discrete response has a slight loss a f = 0 of about 0.3 dB, and the insertion loss is less than that of the continuous-time system as frequency approaches f_{s}/2. Accuracy of the frequency response can be improved by increasing the sample rate.

Placing a zero at z = 0 + j0 gives numerator coefficients b = [0 1]. This has no effect on the impulse response amplitude, but serves to align the timing with that of the continuous-time response.

Figure 9. Discrete-time (red) and continuous-time (blue) responses, impulse-invariant,

ζ = 0.707, f_{n}= 10 Hz, and f_{s} = 100 Hz. a. Impulse response. b. Freq. response magnitude.

## Appendix Matlab Function so_demo.m

The function so_demo.m converts the second-order system H(s) described in Equation 1 to a discrete-time system H(z), and plots the system’s poles, impulse response, and frequency response. H(z) is shown in Equation 2. The function call is:

[K,b,a]= so_demo(zeta,fn,fs,nz);

where

zeta = damping ratio. zeta > 0 for a stable system.

fn = natural frequency, Hz

fs = sample frequency, Hz

nz = number of zeros (1 or 2)

K = numerator scale factor of H(z)

b = vector of numerator coefficients of H(z)

a = vector of denominator coefficients of H(z)

nz must be 1 or 2:

- For nz = 1, the impulse-invariance method [6] gives a zero at z = 0 + j0. It turns out that K =
T
_{s}*h(T_{s}), where h(T_{s}) is the continuous-time system’s impulse response evaluated at t = T_{s}. The formula for h(T_{s}) varies, depending on whether ζ < 1, ζ = 1, or ζ > 1. Note that if we didn’t care about exact amplitude scaling, we could have instead set K = sum(a)/sum(b) = sum(a). - For nz = 2, the zeros are at z = -1 + j0.

Table 1 shows how the function handles the two cases of nz.

For a listing of so_demo.m, see the attached pdf file.

## References

1. Boulet, Benoit, *Fundamentals of Signals and Systems*, Thomson Learning, 2006, p 300 – 305.

https://mlichouri.files.wordpress.com/2013/10/fundamentals-of-signals-and-systems.pdf__
__

2. Sachs, Jason, “Second-Order Systems, Part I:Boing!”, Embeddedrelated website, 2014,

https://www.embeddedrelated.com/showarticle/671.php

3. *Signal Processing Toolbox User’s Guide*, MathWorks, 2018, p 3.

4.* Signal Processing Toolbox User’s Guide*, p 6.

5. Robertson, Neil, “Design IIR Butterworth Filters Using 12 Lines of Code”, DSPrelated website, 2018,

https://www.dsprelated.com/showarticle/1119.php

6. Lyons, Richard G., *Understanding Digital Signal Processing*, 3^{rd} Ed., section 6.10.

__
__

__
__

__
__

Neil Robertson April 1, 2020

**Previous post by Neil Robertson:**

A Simplified Matlab Function for Power Spectral Density

- Comments
- Write a Comment Select to add a comment

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.

Registering will allow you to participate to the forums on ALL the related sites and give you access to all pdf downloads.