# Transformation Matrix for Time Discretization Based on Tustin’s Method

Authors: Yiming Jiang; Xiaodong Hu; Sen Wu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper studies rules in transformation of transfer function through time discretization. A method of using transformation matrix to realize bilinear transform (also known as Tustin’s method) is presented. This method can be described as the conversion between the coefficients of transfer functions, which are expressed as transform by certain matrix. For a polynomial of degreen, the corresponding transformation matrix of order n exists and is unique. Furthermore, the transformation matrix can be decomposed into an upper triangular matrix multiplied with another lower triangular matrix. And both have obvious regularity. The proposed method can achieve rapid bilinear transform used in automatic design of digital filter. The result of numerical simulation verifies the correctness of the theoretical results. Moreover, it also can be extended to other similar problems. Example in the last throws light on this point.


---

## Body

## 1. Introduction

In this paper we proposed a new method of using transformation matrix to realize bilinear transform. We describe the principles and characteristics of the method in Section2. Also we have examinations and simulations on the method in Section 3.In Section2, Section 2.1 starts with exploration of how to express derivative of a function with discrete points under the method of bilinear transform. By studying the rules of that conversion, Section 2.2 reveals regularity in the arithmetic of bilinear transform operated on transfer function, so that it can be replaced by systematic mathematical derivations. Finally, the transformation matrix to realize Tustin transform is obtained by recursion. Also in the process, some inner rules of the Tustin transform can be revealed.In Section3, Section 3.1 provides a simple application to verify the mathematical correctness of the method and uses the obtained transformation matrix as an example. Section 3.2 explains how it works when the parameters of functions are uncertain. Section 3.3 shows an example about its expanded application as it can be used to construct transformation matrixes in other situations. The purpose is to enrich the method and to make it easier to understand. Of course they are all simplified.


### 1.1. About Tustin’s Method

Tustin’s method is also called the bilinear transform. As a convenient tool for time discretization, it is widely used in designing digital filters to obtain expected discrete time system with desired response characteristics. The transform, as a method to convertS domain function to Z domain function, is a first-order approximation of the S-to-Z mapping. Under the Laplace transform, it can be regarded as a discrete-time approximation [1]:
(1)
z
=
e
T
s
=
e
T
s
/
2
e
-
T
s
/
2
=
∑
n
=
0
∞
(
1
/
n
!
)
·
(
T
s
/
2
)
n
∑
n
=
0
∞
(
1
/
n
!
)
·
(
-
T
s
/
2
)
n
≈
1
+
T
s
/
2
1
-
T
s
/
2
.
Equivalently,
(2)
s
=
2
T
·
1
-
z
-
1
1
+
z
-
1
,
where T is the numerical integration step size of the trapezoidal rule used in the bilinear transform derivation.Compared withZ transform based on the pulse response method, the Tustin transform eliminates the aliasing effect generated by overlapping of frequency spectrum, relying on the linear transformation from S domain to Z domain as a one-to-one mapping. Consequently, it is suitable for the design of high-pass and band-stop filters. Furthermore, because of the simple algebraic relationship between S and Z, the transfer function of digital filter can be obtained conveniently by directing the S-to-Z mapping into the corresponding transfer function of analog filter [1–3].


### 1.2. Problems to Be Solved

Because the bilinear transform is a kind of fractional linear transformation, the amount of calculation increases sharply with the increase of the order of system. The problem of complex calculation appears when operating bilinear transform for high-order system. Particularly, in many situations the parameters of system are changeable; thus then traditional method fails to work. For these reasons, programs to realize automatic operation are required. An efficient way is to express the transform with matrix.Many scholars have proposed methods to obtain the matrix for bilinear transform, of which the literatures [4–6] were to derive coefficient matrix of different order under the discussion of relationships between the matrix elements, while [7, 8] used method of induction to derive coefficient matrix of order n from coefficient matrix of order zero. Other relative reference as [9–14] also studied quick algorithms for different forms of linear S-to-Z transform.The purpose of this paper is to introduce an easier way to obtain the transformation matrix. The method applies when the order of transfer function is large, using computer as a direct operation, or when the system function is uncertain, using the method to get the dynamic differential equation.


## 1.1. About Tustin’s Method

Tustin’s method is also called the bilinear transform. As a convenient tool for time discretization, it is widely used in designing digital filters to obtain expected discrete time system with desired response characteristics. The transform, as a method to convertS domain function to Z domain function, is a first-order approximation of the S-to-Z mapping. Under the Laplace transform, it can be regarded as a discrete-time approximation [1]:
(1)
z
=
e
T
s
=
e
T
s
/
2
e
-
T
s
/
2
=
∑
n
=
0
∞
(
1
/
n
!
)
·
(
T
s
/
2
)
n
∑
n
=
0
∞
(
1
/
n
!
)
·
(
-
T
s
/
2
)
n
≈
1
+
T
s
/
2
1
-
T
s
/
2
.
Equivalently,
(2)
s
=
2
T
·
1
-
z
-
1
1
+
z
-
1
,
where T is the numerical integration step size of the trapezoidal rule used in the bilinear transform derivation.Compared withZ transform based on the pulse response method, the Tustin transform eliminates the aliasing effect generated by overlapping of frequency spectrum, relying on the linear transformation from S domain to Z domain as a one-to-one mapping. Consequently, it is suitable for the design of high-pass and band-stop filters. Furthermore, because of the simple algebraic relationship between S and Z, the transfer function of digital filter can be obtained conveniently by directing the S-to-Z mapping into the corresponding transfer function of analog filter [1–3].


## 1.2. Problems to Be Solved

Because the bilinear transform is a kind of fractional linear transformation, the amount of calculation increases sharply with the increase of the order of system. The problem of complex calculation appears when operating bilinear transform for high-order system. Particularly, in many situations the parameters of system are changeable; thus then traditional method fails to work. For these reasons, programs to realize automatic operation are required. An efficient way is to express the transform with matrix.Many scholars have proposed methods to obtain the matrix for bilinear transform, of which the literatures [4–6] were to derive coefficient matrix of different order under the discussion of relationships between the matrix elements, while [7, 8] used method of induction to derive coefficient matrix of order n from coefficient matrix of order zero. Other relative reference as [9–14] also studied quick algorithms for different forms of linear S-to-Z transform.The purpose of this paper is to introduce an easier way to obtain the transformation matrix. The method applies when the order of transfer function is large, using computer as a direct operation, or when the system function is uncertain, using the method to get the dynamic differential equation.


## 2. Derivations

Bilinear transform is a bridge betweenS domain and Z domain. Laplace transform can change a time domain function witht into a frequency domain function with s, and it has a special quality that it can deform differential equations into algebraic equations to simplify mathematical operations. Z transform is similar to Laplace transform, but it is used to deal with discrete series of t. It can change a time domain function with t into Z domain. The result of every mapping from S domain to Z domain is replacement of continuous function by discrete points. Accordingly, when processing a continuous function, the result of mapping is expression of relations between numerical values, integration, and differentiation by discrete points. Because transfer functions can be converted to differential equations, expressions of differential forms by discrete points should be studied.


### 2.1. Definition for Derivatives with Discrete Points

In order to discuss distinctions between deformations to differential forms before and afterS-to-Z mapping by bilinear transform, such relations below will be necessary to deduct.Suppose{
x
n
,
x
n
-
1
,
x
n
-
2
,
…
,
x
n
-
m
} are (
m
+
1
) values at sampling points of discrete digital signal x
(
t
) (m is an integer less than or equal to n). If the sampling step size is T, then x
n
-
i
=
x
(
n
T
-
i
T
), i
=
0,1
,
2
,
…
,
m.If(
x
n
-
x
n
-
1
)
/
T represents first-order derivative with discrete values x
n, x
n
-
1 of two adjacent points, such recursion can be obtained as (
x
n
-
2
x
n
-
1
+
x
n
-
2
)
/
T
2 represents the second-order derivative with three, (
x
n
-
3
x
n
-
1
+
3
x
n
-
2
-
x
n
-
3
)
/
T
3 represents the third-order derivative with four, and so on. Here we define
(3)
d
i
x
→
n
=
∑
j
=
0
i
(
-
1
)
j
C
i
j
x
n
-
j
,
where d
i and x
→
n are (
m
+
1
)-dimensional row and column vectors and (
m
+
1
) is equal the number of samples. Let C
i
j
=
i
!
/
(
(
i
-
j
)
!
j
!
) (i
,
j are nonnegative integers, and i
≥
j):
(4)
d
i
=
[
C
i
0
-
C
i
1
⋯
(
-
1
)
i
C
i
i
0
⋯
0
]
,
x
→
n
=
[
x
n
x
n
-
1
⋯
x
n
-
m
]
T
.Thusd
i
x
→
n
/
T
i is the expression of the ith order derivative d
i
x
(
t
)
/
d
t
i for continuous time signal x with (
i
+
1
) values. In addition, each sampling point can be expressed by the linear combination of d
i
x
→
n
(
i
=
0,1
,
2
,
…
) as below:
(5)
x
n
=
x
n
,
x
n
-
1
=
x
n
-
(
x
n
-
x
n
-
1
)
,
x
n
-
2
=
x
n
-
2
(
x
n
-
x
n
-
1
)
+
(
x
n
-
2
x
n
-
1
+
x
n
-
2
)
,
x
n
-
3
=
x
n
-
3
(
x
n
-
x
n
-
1
)
+
3
(
x
n
-
2
x
n
-
1
+
x
n
-
2
)
000000
-
(
x
n
-
3
x
n
-
1
+
3
x
n
-
2
-
x
n
-
3
)
⋯
.
That is,
(6)
x
n
-
m
=
∑
i
=
0
m
(
-
1
)
i
C
m
i
(
∑
j
=
0
i
(
-
1
)
j
C
i
j
x
n
-
j
)
=
∑
i
=
0
m
(
-
1
)
i
C
m
i
(
d
i
x
→
n
)
.To facilitate the expression, let(7)
[
C
0
0
C
1
0
C
2
0
⋯
C
m
0
-
C
1
1
-
C
2
1
⋯
-
C
m
1
C
2
2
⋯
C
m
2
⋱
⋮
(
-
1
)
m
C
m
m
]
=
D
m
.D
m is an upper triangular matrix of (
m
+
1
)
×
(
m
+
1
) and satisfies the following properties:
(8)
D
m
=
(
D
m
)
-
1
=
[
d
0
T
d
1
T
⋯
d
m
T
]
,
(9)
x
→
n
T
=
[
d
0
x
→
n
d
1
x
→
n
⋯
d
m
x
→
n
]
D
m
.Extensionally, ifx
→
n
-
i is defined as an (
m
+
1
)-dimensional column vector [
x
n
-
i
x
n
-
i
-
1
⋯
x
n
-
m
+
1
x
n
-
m
0
⋯
0
]
T by moving forward the elements of x
→
n for i steps, we get another property (10) as
(10)
[
d
m
-
i
x
→
n
d
m
-
i
x
→
n
-
1
⋯
d
m
-
i
x
→
n
-
i
]
=
[
d
m
-
i
x
→
n
d
m
-
i
+
1
x
→
n
⋯
d
m
x
→
n
]
D
i
.For example, whenm
=
3, i
=
2,
(11)
D
3
=
[
1
1
1
1
-
1
-
2
-
3
1
3
-
1
]
,
D
2
=
[
1
1
1
-
1
-
2
1
]
,
x
→
n
T
=
[
x
n
x
n
-
1
x
n
-
2
x
n
-
3
]
=
[
d
0
x
→
n
T
d
1
x
→
n
T
d
2
x
→
n
T
d
3
x
→
n
T
]
[
1
1
1
1
-
1
-
2
-
3
1
3
-
1
]
,
[
d
1
x
→
n
T
d
1
x
→
n
-
1
T
d
1
x
→
n
-
2
T
]
=
[
(
x
n
-
x
n
-
1
)
(
x
n
-
1
-
x
n
-
2
)
(
x
n
-
2
-
x
n
-
3
)
]
=
[
d
1
x
→
n
T
d
2
x
→
n
T
d
3
x
→
n
T
]
[
1
1
1
-
1
-
2
1
]
,
where, d
0
=
[
1
0
0
0
], d
1
=
[
1
-
1
0
0
], d
2
=
[
1
-
2
1
0
], d
3
=
[
1
-
3
3
-
1
], and x
→
n
T
=
[
x
n
x
n
-
1
x
n
-
2
x
n
-
3
]
T, x
→
n
-
1
T
=
[
x
n
-
1
x
n
-
2
x
n
-
3
0
]
T, and x
→
n
-
2
T
=
[
x
n
-
2
x
n
-
3
0
0
]
T.


### 2.2. Principles and Steps to Obtain Transformation Matrix

With the definitions above, the transformation matrix for bilinear transform will be presented and the relations between coefficients before and after the transform will be deducted.Since the digital signal is not continuous,d
i
x
→
n
/
T
i is used as the ith order derivative of noncontinuous signal x. For a general transfer function,
(12)
H
(
s
)
=
Y
(
s
)
X
(
s
)
=
a
0
+
a
1
s
+
⋯
+
a
p
-
1
s
p
-
1
+
a
p
s
p
b
0
+
b
1
s
+
⋯
+
b
p
-
1
s
p
-
1
+
b
p
s
pis equal to the differential equation as follows:(13)
b
0
y
(
t
)
+
b
1
d
y
(
t
)
d
t
+
⋯
+
b
p
d
p
y
(
t
)
d
t
p
=
a
0
x
(
t
)
+
a
1
d
x
(
t
)
d
t
+
⋯
+
a
p
d
p
x
(
t
)
d
t
p
,
where, X
(
s
) and Y
(
s
) are polynomials of s, and b
p
≠
0. The general method of using bilinear transformation is substituting the S-to-Z mapping (2) into the transfer function. The Z domain transfer function H
(
z
) can be shown as follows:
(14)
H
(
z
)
=
α
0
+
α
1
z
-
1
+
⋯
+
α
p
-
1
z
-
(
p
-
1
)
+
α
p
z
-
p
β
0
+
β
1
z
-
1
+
⋯
+
β
p
-
1
z
-
(
p
-
1
)
+
β
p
z
-
p
.According to the definition ofZ transform, (14) is equal to the difference equation as follows:
(15)
β
0
y
n
+
β
1
y
n
-
1
+
⋯
+
β
p
y
n
-
p
=
α
0
x
n
+
α
1
x
n
-
1
+
⋯
+
α
p
x
n
-
p
.Here definea
→
=
[
a
0
a
1
⋯
a
p
]
T, b
→
=
[
b
0
b
1
⋯
b
p
]
T, α
→
=
[
α
0
α
1
⋯
α
p
]
T, and β
→
=
[
β
0
β
1
⋯
β
p
]
T. Combined with the relation between X and Y obtained by the transfer function H
(
s
), equation of X and Y can be obtained. The following are deduced on the relationship between a
→ and b
→, α
→, and β
→.Consider the simple case whena
→
=
[
1
0
]
T, b
→
=
[
0
1
]
T, and H
(
s
)
=
Y
(
s
)
/
X
(
s
)
=
1
/
s; then
(16)
s
Y
(
s
)
-
X
(
s
)
=
0
.According to the Laplace transform, (16) is equal to (17):
(17)
d
y
(
t
)
d
t
=
x
(
t
)
.Take input sequence{
x
n
,
x
n
-
1
} from signal x and output sequence {
y
n
,
y
n
-
1
} from signal, and the step size is T. Meanwhile, the average values are x
=
(
1
/
2
)
(
x
n
+
x
n
-
1
) and y
=
(
1
/
2
)
(
y
n
+
y
n
-
1
) weighed by {
C
1
0
/
2
,
C
1
1
/
2
}.Replace the first-order derivative ofy by d
1
y
→
n
/
T; according to (17), d
1
y
→
n
/
T equals x
-, as
(18)
y
n
-
y
n
-
1
T
=
x
=
1
2
(
x
n
+
x
n
-
1
)
which can be expressed as
(19)
d
1
y
→
n
T
=
1
2
(
d
0
x
→
n
+
d
0
x
→
n
-
1
)
=
1
2
∑
i
=
0
1
(
C
1
i
d
0
x
→
n
-
i
)
.Thus it can be concluded that when taking input sequence{
x
n
,
x
n
-
1
} and output sequence {
y
n
,
y
n
-
1
} into related operations, as a kind of approximation, (
1
/
2
)
(
d
0
x
→
n
+
d
0
x
→
n
-
1
) takes the position of x (d
0
x
→
n
=
x
n and d
0
x
→
n
-
1
=
x
n
-
1), and d
1
y
→
n
/
T takes the position of the first-order derivative ofy. The mapping can be shown as below
(20)
d
i
k
(
t
)
d
t
i
⟶
∑
j
=
0
1
-
i
(
C
1
-
i
j
d
i
k
→
n
-
j
)
2
1
-
i
T
i
,
(
i
=
0,1
)Take input sequence{
x
n
,
x
n
-
1
,
x
n
-
2
} and output sequence {
y
n
,
y
n
-
1
,
y
n
-
2
} sequentially. According to (17),
(21)
y
n
-
y
n
-
1
T
=
1
2
(
x
n
+
x
n
-
1
)
,
(22)
y
n
-
1
-
y
n
-
2
T
=
1
2
(
x
n
-
1
+
x
n
-
2
)
.Subtract (22) from (21) and then multiply by (
T
)
-
1; we get
(23)
y
n
-
2
y
n
-
1
+
y
n
-
2
T
2
=
1
2
(
x
n
-
x
n
-
2
)
T
=
1
2
(
x
n
-
x
n
-
1
)
+
(
x
n
-
1
-
x
n
-
2
)
T
as
(24)
d
2
y
→
n
T
2
=
1
2
(
d
1
x
→
n
+
d
1
x
→
n
-
1
)
T
=
1
2
T
∑
i
=
0
1
(
C
1
i
d
1
x
→
n
-
i
)
.Becaused
2
y
→
n
/
T
2 is the expression of second-order derivative, (
1
/
2
)
(
(
d
1
x
→
n
+
d
1
x
→
n
-
1
)
/
T
) is the average first-order derivative weighed by {
C
1
0
/
2
,
C
1
1
/
2
}. Similar to (17), (24) can be expressed by the differential equation as follows:
(25)
d
2
y
(
t
)
d
t
2
=
d
x
(
t
)
d
t
.After the Laplace transform, (26) can be obtained which can be deduced by (16):
(26)
s
2
Y
(
s
)
-
s
X
(
s
)
=
0
.Add (22) to (21) and we get
(27)
y
n
-
y
n
-
2
T
=
1
2
(
x
n
+
2
x
n
-
1
+
x
n
-
2
)
as
(28)
1
2
(
d
1
y
→
n
+
d
1
y
→
n
-
1
)
T
=
1
4
(
d
0
x
→
n
+
2
d
0
x
→
n
-
1
+
d
0
x
→
n
-
2
)
=
1
4
∑
i
=
0
2
(
C
2
i
d
0
x
→
n
-
i
)
.Because(
1
/
2
)
(
(
d
1
y
→
n
+
d
1
y
→
n
-
1
)
/
T
) is the expression of first-order derivative d
y
(
t
)
/
d
t, (
1
/
4
)
(
d
0
x
→
n
+
2
d
0
x
→
n
-
1
+
d
0
x
→
n
-
2
) is the average value of {
x
n
,
x
n
-
1
,
x
n
-
2
} weighed by {
C
2
0
/
4
,
C
2
1
/
4
,
C
2
2
/
4
}. Equation (28) can be expressed by the differential equation as (16).Thus it can be concluded that, when taking input sequence{
x
n
,
x
n
-
1
,
x
n
-
2
} and output sequence {
y
n
,
y
n
-
1
,
y
n
-
2
} into related operations, (
1
/
4
)
(
d
0
x
→
n
+
2
d
0
x
→
n
-
1
+
d
0
x
→
n
-
2
) takes the position of x, (
1
/
2
)
(
(
d
1
y
→
n
+
d
1
y
→
n
-
1
)
/
T
) takes the position of the first-order derivative of y, and d
2
y
→
n
/
T
2 takes the position of the second-order derivative of y. The mapping can be shown as follows:
(29)
d
i
k
(
t
)
d
t
i
⟶
∑
j
=
0
2
-
i
(
C
2
-
i
j
d
i
k
→
n
-
j
)
2
2
-
i
T
i
,
(
i
=
0,1
,
2
)
.Now consider whena
→
=
[
1
0
0
]
T, b
→
=
[
b
0
b
1
b
2
]
T, and H
(
s
)
=
1
/
(
b
0
+
b
1
s
+
b
2
s
2
). The corresponding differential equation is
(30)
b
0
y
(
t
)
+
b
1
d
y
(
t
)
d
t
+
b
2
d
2
y
(
t
)
d
t
2
=
x
(
t
)
.With conclusions above, taking input sequence{
x
n
,
x
n
-
1
,
x
n
-
2
} and output sequence {
y
n
,
y
n
-
1
,
y
n
-
2
}, (30) transforms into the following:
(31)
b
0
(
y
n
+
2
y
n
-
1
+
y
n
-
2
4
)
+
b
1
(
d
1
y
→
n
+
d
1
y
→
n
-
1
2
T
)
+
b
2
d
2
y
→
n
T
2
=
x
n
+
2
x
n
-
1
+
x
n
-
2
4
.The coefficients of{
x
n
,
x
n
-
1
,
x
n
-
2
} and {
y
n
,
y
n
-
1
,
y
n
-
2
} in (31) are the same as coefficients in the difference equation after bilinear transform.By this method, take input sequence{
x
n
,
x
n
-
1
,
…
,
x
n
-
m
} and output sequence {
y
n
,
y
n
-
1
,
…
,
y
n
-
m
} with (
m
+
1
) points, and the mapping can be expressed as
(32)
d
i
k
(
t
)
d
t
i
⟶
∑
j
=
0
m
-
i
(
C
m
-
i
j
d
i
k
→
n
-
j
)
2
m
-
i
T
i
,
(
i
=
0,1
,
2
,
…
,
m
)
.Such rules can be used to deal with a general transfer function (12), for (13) can be expressed as
(33)
b
0
2
i
∑
i
=
0
p
(
C
p
i
d
0
y
→
n
-
i
)
+
b
1
2
i
-
1
T
∑
i
=
0
p
-
1
(
C
p
-
1
i
d
1
y
→
n
-
i
)
+
⋯
+
b
p
T
i
d
p
y
→
n
=
a
0
2
i
∑
i
=
0
p
(
C
p
i
d
0
x
→
n
-
i
)
+
a
1
2
i
-
1
T
∑
i
=
0
p
-
1
(
C
p
-
1
i
d
1
x
→
n
-
i
)
+
⋯
+
a
p
T
i
d
p
x
→
n
.As a conclusion, for a discrete digital signalk whose interval is T, in order to get its mth order derivative, at least (
m
+
1
) values are needed to take into operations. If themth order derivative of k is expressed by (
m
+
l
) values, such mapping can be shown as follows:
(34)
d
m
k
(
t
)
d
t
m
⟶
T
-
m
·
2
-
(
l
-
1
)
∑
i
=
0
l
-
1
(
C
l
-
1
i
(
d
m
k
→
n
-
i
)
)
.Asl changes, the mapping shown by (34) satisfies the distribution rule of binomial coefficients [4, 5], shown in Figure 1.Figure 1
Using(
m
+
l
) values to express mth order derivative of k.The tree structure in Figure1 can be expressed by Figure 2, in which k
i
j is the element in ith row and jth column. Meanwhile, (34) can also be expressed as d
m
k
(
t
)
/
d
t
m
→
T
-
m
·
2
-
(
l
-
1
)
∑
j
=
1
l
(
k
l
j
(
d
m
k
→
n
-
j
+
1
)
).Figure 2
Relations between elements of Figure1 [4].For example, according to that theory, while taking{
k
n
,
k
n
-
1
,
…
,
k
n
-
m
} from signal k, according to (34), (
T
)
-
m
d
m
k
→
n represents the mth order derivative of k; 2
-
i
(
T
)
-
(
m
-
i
)
∑
j
=
0
i
(
C
i
j
d
m
-
i
k
→
n
-
j
) represents the (
m
-
i
)th order derivative of k (i
≤
m). And according to (10), such equation exists as follows:
(35)
∑
j
=
0
i
(
C
i
j
d
m
-
i
k
→
n
-
j
)
=
[
d
m
-
i
k
→
n
d
m
-
i
+
1
k
→
n
⋯
d
m
k
→
n
]
D
i
[
C
i
0
C
i
1
⋮
C
i
i
]
.For example, whenm
=
3, i
=
2,
(36)
∑
j
=
0
2
(
C
2
j
d
1
k
→
n
-
j
)
=
[
d
1
k
→
n
d
1
k
→
n
-
1
d
1
k
→
n
-
2
]
[
1
2
1
]
=
[
d
1
k
→
n
d
2
k
→
n
d
3
k
→
n
]
D
2
[
1
2
1
]
=
[
d
1
k
→
n
d
2
k
→
n
d
3
k
→
n
]
[
4
-
4
1
]
.Back to the deduction on the relationship betweena
→ and b
→, and α
→ and β
→, first multiply both sides by 2
p; referring to the conclusions above, the right side of (12) can be transformed by such mapping:
(37)
2
p
(
a
0
x
(
t
)
+
a
1
d
x
(
t
)
d
t
+
⋯
+
a
p
d
p
x
(
t
)
d
x
p
)
⟶
∑
i
=
0
p
(
a
i
(
2
T
)
(
p
-
i
)
∑
j
=
0
i
(
C
i
j
d
p
-
i
x
→
n
-
j
)
)
.According to (34), the equation can be shown as follows:
(38)
∑
i
=
0
p
(
a
i
(
2
T
)
(
p
-
i
)
∑
j
=
0
i
(
C
i
j
d
p
-
i
x
→
n
-
j
)
)
=
[
d
0
x
→
n
d
1
x
→
n
⋯
d
p
x
→
n
]
H
p
a
→
,
where H
p is lower triangular matrix of (
p
+
1
)
×
(
p
+
1
); H
p is defined as
(39)
H
p
=
[
h
→
p
(
2
T
)
h
→
p
-
1
⋯
(
2
T
)
(
p
-
i
)
h
→
p
-
i
⋯
(
2
T
)
p
h
→
0
]
.
h
→
0
,
h
→
1
,
…
,
h
→
p are (
p
+
1
)-dimensional column vectors, and
(40)
h
i
=
[
O
D
i
]
[
C
i
0
C
i
1
⋯
C
i
i
]
T
,
where i
=
0,1
,
2
,
…
,
p, [
O
D
i
] is partitioned matrix of (
p
+
1
)
×
(
i
+
1
) and O is zero matrix of (
p
-
i
)
×
(
i
+
1
).According to (9), the right side of (15) transforms as follows:
(41)
α
0
x
n
+
α
1
x
n
-
1
+
⋯
+
α
p
x
n
-
p
=
x
→
n
T
α
→
=
[
d
0
x
→
n
d
1
x
→
n
⋯
d
p
x
→
n
]
D
p
α
→
.ObviouslyD
p
α
→ is the coefficient vector of [
d
0
x
→
n
d
1
x
→
n
⋯
d
p
x
→
n
]. Since (13) equals (15) and (38) equals (41), consequently,
(42)
D
p
α
→
=
H
p
a
→
.In the same way,(43)
D
p
β
→
=
H
p
b
→
.According to (8), (
D
p
)
-
1
=
D
p, the coefficient of (15) can be expressed as follows:
(44)
α
→
=
D
p
H
p
a
→
,
β
→
=
D
p
H
p
b
→
,
where a
→ and b
→ express the coefficients of the numerator and denominator polynomials of H
(
s
). The transformation matrix for bilinear transformation is obtained by multiplying D
p and H
p, D
p, and H
p exist and are unique when p is set to a fixed value.


## 2.1. Definition for Derivatives with Discrete Points

In order to discuss distinctions between deformations to differential forms before and afterS-to-Z mapping by bilinear transform, such relations below will be necessary to deduct.Suppose{
x
n
,
x
n
-
1
,
x
n
-
2
,
…
,
x
n
-
m
} are (
m
+
1
) values at sampling points of discrete digital signal x
(
t
) (m is an integer less than or equal to n). If the sampling step size is T, then x
n
-
i
=
x
(
n
T
-
i
T
), i
=
0,1
,
2
,
…
,
m.If(
x
n
-
x
n
-
1
)
/
T represents first-order derivative with discrete values x
n, x
n
-
1 of two adjacent points, such recursion can be obtained as (
x
n
-
2
x
n
-
1
+
x
n
-
2
)
/
T
2 represents the second-order derivative with three, (
x
n
-
3
x
n
-
1
+
3
x
n
-
2
-
x
n
-
3
)
/
T
3 represents the third-order derivative with four, and so on. Here we define
(3)
d
i
x
→
n
=
∑
j
=
0
i
(
-
1
)
j
C
i
j
x
n
-
j
,
where d
i and x
→
n are (
m
+
1
)-dimensional row and column vectors and (
m
+
1
) is equal the number of samples. Let C
i
j
=
i
!
/
(
(
i
-
j
)
!
j
!
) (i
,
j are nonnegative integers, and i
≥
j):
(4)
d
i
=
[
C
i
0
-
C
i
1
⋯
(
-
1
)
i
C
i
i
0
⋯
0
]
,
x
→
n
=
[
x
n
x
n
-
1
⋯
x
n
-
m
]
T
.Thusd
i
x
→
n
/
T
i is the expression of the ith order derivative d
i
x
(
t
)
/
d
t
i for continuous time signal x with (
i
+
1
) values. In addition, each sampling point can be expressed by the linear combination of d
i
x
→
n
(
i
=
0,1
,
2
,
…
) as below:
(5)
x
n
=
x
n
,
x
n
-
1
=
x
n
-
(
x
n
-
x
n
-
1
)
,
x
n
-
2
=
x
n
-
2
(
x
n
-
x
n
-
1
)
+
(
x
n
-
2
x
n
-
1
+
x
n
-
2
)
,
x
n
-
3
=
x
n
-
3
(
x
n
-
x
n
-
1
)
+
3
(
x
n
-
2
x
n
-
1
+
x
n
-
2
)
000000
-
(
x
n
-
3
x
n
-
1
+
3
x
n
-
2
-
x
n
-
3
)
⋯
.
That is,
(6)
x
n
-
m
=
∑
i
=
0
m
(
-
1
)
i
C
m
i
(
∑
j
=
0
i
(
-
1
)
j
C
i
j
x
n
-
j
)
=
∑
i
=
0
m
(
-
1
)
i
C
m
i
(
d
i
x
→
n
)
.To facilitate the expression, let(7)
[
C
0
0
C
1
0
C
2
0
⋯
C
m
0
-
C
1
1
-
C
2
1
⋯
-
C
m
1
C
2
2
⋯
C
m
2
⋱
⋮
(
-
1
)
m
C
m
m
]
=
D
m
.D
m is an upper triangular matrix of (
m
+
1
)
×
(
m
+
1
) and satisfies the following properties:
(8)
D
m
=
(
D
m
)
-
1
=
[
d
0
T
d
1
T
⋯
d
m
T
]
,
(9)
x
→
n
T
=
[
d
0
x
→
n
d
1
x
→
n
⋯
d
m
x
→
n
]
D
m
.Extensionally, ifx
→
n
-
i is defined as an (
m
+
1
)-dimensional column vector [
x
n
-
i
x
n
-
i
-
1
⋯
x
n
-
m
+
1
x
n
-
m
0
⋯
0
]
T by moving forward the elements of x
→
n for i steps, we get another property (10) as
(10)
[
d
m
-
i
x
→
n
d
m
-
i
x
→
n
-
1
⋯
d
m
-
i
x
→
n
-
i
]
=
[
d
m
-
i
x
→
n
d
m
-
i
+
1
x
→
n
⋯
d
m
x
→
n
]
D
i
.For example, whenm
=
3, i
=
2,
(11)
D
3
=
[
1
1
1
1
-
1
-
2
-
3
1
3
-
1
]
,
D
2
=
[
1
1
1
-
1
-
2
1
]
,
x
→
n
T
=
[
x
n
x
n
-
1
x
n
-
2
x
n
-
3
]
=
[
d
0
x
→
n
T
d
1
x
→
n
T
d
2
x
→
n
T
d
3
x
→
n
T
]
[
1
1
1
1
-
1
-
2
-
3
1
3
-
1
]
,
[
d
1
x
→
n
T
d
1
x
→
n
-
1
T
d
1
x
→
n
-
2
T
]
=
[
(
x
n
-
x
n
-
1
)
(
x
n
-
1
-
x
n
-
2
)
(
x
n
-
2
-
x
n
-
3
)
]
=
[
d
1
x
→
n
T
d
2
x
→
n
T
d
3
x
→
n
T
]
[
1
1
1
-
1
-
2
1
]
,
where, d
0
=
[
1
0
0
0
], d
1
=
[
1
-
1
0
0
], d
2
=
[
1
-
2
1
0
], d
3
=
[
1
-
3
3
-
1
], and x
→
n
T
=
[
x
n
x
n
-
1
x
n
-
2
x
n
-
3
]
T, x
→
n
-
1
T
=
[
x
n
-
1
x
n
-
2
x
n
-
3
0
]
T, and x
→
n
-
2
T
=
[
x
n
-
2
x
n
-
3
0
0
]
T.


## 2.2. Principles and Steps to Obtain Transformation Matrix

With the definitions above, the transformation matrix for bilinear transform will be presented and the relations between coefficients before and after the transform will be deducted.Since the digital signal is not continuous,d
i
x
→
n
/
T
i is used as the ith order derivative of noncontinuous signal x. For a general transfer function,
(12)
H
(
s
)
=
Y
(
s
)
X
(
s
)
=
a
0
+
a
1
s
+
⋯
+
a
p
-
1
s
p
-
1
+
a
p
s
p
b
0
+
b
1
s
+
⋯
+
b
p
-
1
s
p
-
1
+
b
p
s
pis equal to the differential equation as follows:(13)
b
0
y
(
t
)
+
b
1
d
y
(
t
)
d
t
+
⋯
+
b
p
d
p
y
(
t
)
d
t
p
=
a
0
x
(
t
)
+
a
1
d
x
(
t
)
d
t
+
⋯
+
a
p
d
p
x
(
t
)
d
t
p
,
where, X
(
s
) and Y
(
s
) are polynomials of s, and b
p
≠
0. The general method of using bilinear transformation is substituting the S-to-Z mapping (2) into the transfer function. The Z domain transfer function H
(
z
) can be shown as follows:
(14)
H
(
z
)
=
α
0
+
α
1
z
-
1
+
⋯
+
α
p
-
1
z
-
(
p
-
1
)
+
α
p
z
-
p
β
0
+
β
1
z
-
1
+
⋯
+
β
p
-
1
z
-
(
p
-
1
)
+
β
p
z
-
p
.According to the definition ofZ transform, (14) is equal to the difference equation as follows:
(15)
β
0
y
n
+
β
1
y
n
-
1
+
⋯
+
β
p
y
n
-
p
=
α
0
x
n
+
α
1
x
n
-
1
+
⋯
+
α
p
x
n
-
p
.Here definea
→
=
[
a
0
a
1
⋯
a
p
]
T, b
→
=
[
b
0
b
1
⋯
b
p
]
T, α
→
=
[
α
0
α
1
⋯
α
p
]
T, and β
→
=
[
β
0
β
1
⋯
β
p
]
T. Combined with the relation between X and Y obtained by the transfer function H
(
s
), equation of X and Y can be obtained. The following are deduced on the relationship between a
→ and b
→, α
→, and β
→.Consider the simple case whena
→
=
[
1
0
]
T, b
→
=
[
0
1
]
T, and H
(
s
)
=
Y
(
s
)
/
X
(
s
)
=
1
/
s; then
(16)
s
Y
(
s
)
-
X
(
s
)
=
0
.According to the Laplace transform, (16) is equal to (17):
(17)
d
y
(
t
)
d
t
=
x
(
t
)
.Take input sequence{
x
n
,
x
n
-
1
} from signal x and output sequence {
y
n
,
y
n
-
1
} from signal, and the step size is T. Meanwhile, the average values are x
=
(
1
/
2
)
(
x
n
+
x
n
-
1
) and y
=
(
1
/
2
)
(
y
n
+
y
n
-
1
) weighed by {
C
1
0
/
2
,
C
1
1
/
2
}.Replace the first-order derivative ofy by d
1
y
→
n
/
T; according to (17), d
1
y
→
n
/
T equals x
-, as
(18)
y
n
-
y
n
-
1
T
=
x
=
1
2
(
x
n
+
x
n
-
1
)
which can be expressed as
(19)
d
1
y
→
n
T
=
1
2
(
d
0
x
→
n
+
d
0
x
→
n
-
1
)
=
1
2
∑
i
=
0
1
(
C
1
i
d
0
x
→
n
-
i
)
.Thus it can be concluded that when taking input sequence{
x
n
,
x
n
-
1
} and output sequence {
y
n
,
y
n
-
1
} into related operations, as a kind of approximation, (
1
/
2
)
(
d
0
x
→
n
+
d
0
x
→
n
-
1
) takes the position of x (d
0
x
→
n
=
x
n and d
0
x
→
n
-
1
=
x
n
-
1), and d
1
y
→
n
/
T takes the position of the first-order derivative ofy. The mapping can be shown as below
(20)
d
i
k
(
t
)
d
t
i
⟶
∑
j
=
0
1
-
i
(
C
1
-
i
j
d
i
k
→
n
-
j
)
2
1
-
i
T
i
,
(
i
=
0,1
)Take input sequence{
x
n
,
x
n
-
1
,
x
n
-
2
} and output sequence {
y
n
,
y
n
-
1
,
y
n
-
2
} sequentially. According to (17),
(21)
y
n
-
y
n
-
1
T
=
1
2
(
x
n
+
x
n
-
1
)
,
(22)
y
n
-
1
-
y
n
-
2
T
=
1
2
(
x
n
-
1
+
x
n
-
2
)
.Subtract (22) from (21) and then multiply by (
T
)
-
1; we get
(23)
y
n
-
2
y
n
-
1
+
y
n
-
2
T
2
=
1
2
(
x
n
-
x
n
-
2
)
T
=
1
2
(
x
n
-
x
n
-
1
)
+
(
x
n
-
1
-
x
n
-
2
)
T
as
(24)
d
2
y
→
n
T
2
=
1
2
(
d
1
x
→
n
+
d
1
x
→
n
-
1
)
T
=
1
2
T
∑
i
=
0
1
(
C
1
i
d
1
x
→
n
-
i
)
.Becaused
2
y
→
n
/
T
2 is the expression of second-order derivative, (
1
/
2
)
(
(
d
1
x
→
n
+
d
1
x
→
n
-
1
)
/
T
) is the average first-order derivative weighed by {
C
1
0
/
2
,
C
1
1
/
2
}. Similar to (17), (24) can be expressed by the differential equation as follows:
(25)
d
2
y
(
t
)
d
t
2
=
d
x
(
t
)
d
t
.After the Laplace transform, (26) can be obtained which can be deduced by (16):
(26)
s
2
Y
(
s
)
-
s
X
(
s
)
=
0
.Add (22) to (21) and we get
(27)
y
n
-
y
n
-
2
T
=
1
2
(
x
n
+
2
x
n
-
1
+
x
n
-
2
)
as
(28)
1
2
(
d
1
y
→
n
+
d
1
y
→
n
-
1
)
T
=
1
4
(
d
0
x
→
n
+
2
d
0
x
→
n
-
1
+
d
0
x
→
n
-
2
)
=
1
4
∑
i
=
0
2
(
C
2
i
d
0
x
→
n
-
i
)
.Because(
1
/
2
)
(
(
d
1
y
→
n
+
d
1
y
→
n
-
1
)
/
T
) is the expression of first-order derivative d
y
(
t
)
/
d
t, (
1
/
4
)
(
d
0
x
→
n
+
2
d
0
x
→
n
-
1
+
d
0
x
→
n
-
2
) is the average value of {
x
n
,
x
n
-
1
,
x
n
-
2
} weighed by {
C
2
0
/
4
,
C
2
1
/
4
,
C
2
2
/
4
}. Equation (28) can be expressed by the differential equation as (16).Thus it can be concluded that, when taking input sequence{
x
n
,
x
n
-
1
,
x
n
-
2
} and output sequence {
y
n
,
y
n
-
1
,
y
n
-
2
} into related operations, (
1
/
4
)
(
d
0
x
→
n
+
2
d
0
x
→
n
-
1
+
d
0
x
→
n
-
2
) takes the position of x, (
1
/
2
)
(
(
d
1
y
→
n
+
d
1
y
→
n
-
1
)
/
T
) takes the position of the first-order derivative of y, and d
2
y
→
n
/
T
2 takes the position of the second-order derivative of y. The mapping can be shown as follows:
(29)
d
i
k
(
t
)
d
t
i
⟶
∑
j
=
0
2
-
i
(
C
2
-
i
j
d
i
k
→
n
-
j
)
2
2
-
i
T
i
,
(
i
=
0,1
,
2
)
.Now consider whena
→
=
[
1
0
0
]
T, b
→
=
[
b
0
b
1
b
2
]
T, and H
(
s
)
=
1
/
(
b
0
+
b
1
s
+
b
2
s
2
). The corresponding differential equation is
(30)
b
0
y
(
t
)
+
b
1
d
y
(
t
)
d
t
+
b
2
d
2
y
(
t
)
d
t
2
=
x
(
t
)
.With conclusions above, taking input sequence{
x
n
,
x
n
-
1
,
x
n
-
2
} and output sequence {
y
n
,
y
n
-
1
,
y
n
-
2
}, (30) transforms into the following:
(31)
b
0
(
y
n
+
2
y
n
-
1
+
y
n
-
2
4
)
+
b
1
(
d
1
y
→
n
+
d
1
y
→
n
-
1
2
T
)
+
b
2
d
2
y
→
n
T
2
=
x
n
+
2
x
n
-
1
+
x
n
-
2
4
.The coefficients of{
x
n
,
x
n
-
1
,
x
n
-
2
} and {
y
n
,
y
n
-
1
,
y
n
-
2
} in (31) are the same as coefficients in the difference equation after bilinear transform.By this method, take input sequence{
x
n
,
x
n
-
1
,
…
,
x
n
-
m
} and output sequence {
y
n
,
y
n
-
1
,
…
,
y
n
-
m
} with (
m
+
1
) points, and the mapping can be expressed as
(32)
d
i
k
(
t
)
d
t
i
⟶
∑
j
=
0
m
-
i
(
C
m
-
i
j
d
i
k
→
n
-
j
)
2
m
-
i
T
i
,
(
i
=
0,1
,
2
,
…
,
m
)
.Such rules can be used to deal with a general transfer function (12), for (13) can be expressed as
(33)
b
0
2
i
∑
i
=
0
p
(
C
p
i
d
0
y
→
n
-
i
)
+
b
1
2
i
-
1
T
∑
i
=
0
p
-
1
(
C
p
-
1
i
d
1
y
→
n
-
i
)
+
⋯
+
b
p
T
i
d
p
y
→
n
=
a
0
2
i
∑
i
=
0
p
(
C
p
i
d
0
x
→
n
-
i
)
+
a
1
2
i
-
1
T
∑
i
=
0
p
-
1
(
C
p
-
1
i
d
1
x
→
n
-
i
)
+
⋯
+
a
p
T
i
d
p
x
→
n
.As a conclusion, for a discrete digital signalk whose interval is T, in order to get its mth order derivative, at least (
m
+
1
) values are needed to take into operations. If themth order derivative of k is expressed by (
m
+
l
) values, such mapping can be shown as follows:
(34)
d
m
k
(
t
)
d
t
m
⟶
T
-
m
·
2
-
(
l
-
1
)
∑
i
=
0
l
-
1
(
C
l
-
1
i
(
d
m
k
→
n
-
i
)
)
.Asl changes, the mapping shown by (34) satisfies the distribution rule of binomial coefficients [4, 5], shown in Figure 1.Figure 1
Using(
m
+
l
) values to express mth order derivative of k.The tree structure in Figure1 can be expressed by Figure 2, in which k
i
j is the element in ith row and jth column. Meanwhile, (34) can also be expressed as d
m
k
(
t
)
/
d
t
m
→
T
-
m
·
2
-
(
l
-
1
)
∑
j
=
1
l
(
k
l
j
(
d
m
k
→
n
-
j
+
1
)
).Figure 2
Relations between elements of Figure1 [4].For example, according to that theory, while taking{
k
n
,
k
n
-
1
,
…
,
k
n
-
m
} from signal k, according to (34), (
T
)
-
m
d
m
k
→
n represents the mth order derivative of k; 2
-
i
(
T
)
-
(
m
-
i
)
∑
j
=
0
i
(
C
i
j
d
m
-
i
k
→
n
-
j
) represents the (
m
-
i
)th order derivative of k (i
≤
m). And according to (10), such equation exists as follows:
(35)
∑
j
=
0
i
(
C
i
j
d
m
-
i
k
→
n
-
j
)
=
[
d
m
-
i
k
→
n
d
m
-
i
+
1
k
→
n
⋯
d
m
k
→
n
]
D
i
[
C
i
0
C
i
1
⋮
C
i
i
]
.For example, whenm
=
3, i
=
2,
(36)
∑
j
=
0
2
(
C
2
j
d
1
k
→
n
-
j
)
=
[
d
1
k
→
n
d
1
k
→
n
-
1
d
1
k
→
n
-
2
]
[
1
2
1
]
=
[
d
1
k
→
n
d
2
k
→
n
d
3
k
→
n
]
D
2
[
1
2
1
]
=
[
d
1
k
→
n
d
2
k
→
n
d
3
k
→
n
]
[
4
-
4
1
]
.Back to the deduction on the relationship betweena
→ and b
→, and α
→ and β
→, first multiply both sides by 2
p; referring to the conclusions above, the right side of (12) can be transformed by such mapping:
(37)
2
p
(
a
0
x
(
t
)
+
a
1
d
x
(
t
)
d
t
+
⋯
+
a
p
d
p
x
(
t
)
d
x
p
)
⟶
∑
i
=
0
p
(
a
i
(
2
T
)
(
p
-
i
)
∑
j
=
0
i
(
C
i
j
d
p
-
i
x
→
n
-
j
)
)
.According to (34), the equation can be shown as follows:
(38)
∑
i
=
0
p
(
a
i
(
2
T
)
(
p
-
i
)
∑
j
=
0
i
(
C
i
j
d
p
-
i
x
→
n
-
j
)
)
=
[
d
0
x
→
n
d
1
x
→
n
⋯
d
p
x
→
n
]
H
p
a
→
,
where H
p is lower triangular matrix of (
p
+
1
)
×
(
p
+
1
); H
p is defined as
(39)
H
p
=
[
h
→
p
(
2
T
)
h
→
p
-
1
⋯
(
2
T
)
(
p
-
i
)
h
→
p
-
i
⋯
(
2
T
)
p
h
→
0
]
.
h
→
0
,
h
→
1
,
…
,
h
→
p are (
p
+
1
)-dimensional column vectors, and
(40)
h
i
=
[
O
D
i
]
[
C
i
0
C
i
1
⋯
C
i
i
]
T
,
where i
=
0,1
,
2
,
…
,
p, [
O
D
i
] is partitioned matrix of (
p
+
1
)
×
(
i
+
1
) and O is zero matrix of (
p
-
i
)
×
(
i
+
1
).According to (9), the right side of (15) transforms as follows:
(41)
α
0
x
n
+
α
1
x
n
-
1
+
⋯
+
α
p
x
n
-
p
=
x
→
n
T
α
→
=
[
d
0
x
→
n
d
1
x
→
n
⋯
d
p
x
→
n
]
D
p
α
→
.ObviouslyD
p
α
→ is the coefficient vector of [
d
0
x
→
n
d
1
x
→
n
⋯
d
p
x
→
n
]. Since (13) equals (15) and (38) equals (41), consequently,
(42)
D
p
α
→
=
H
p
a
→
.In the same way,(43)
D
p
β
→
=
H
p
b
→
.According to (8), (
D
p
)
-
1
=
D
p, the coefficient of (15) can be expressed as follows:
(44)
α
→
=
D
p
H
p
a
→
,
β
→
=
D
p
H
p
b
→
,
where a
→ and b
→ express the coefficients of the numerator and denominator polynomials of H
(
s
). The transformation matrix for bilinear transformation is obtained by multiplying D
p and H
p, D
p, and H
p exist and are unique when p is set to a fixed value.


## 3. Examples

### 3.1. An Example of Digital Filter Design

For example, the design of a 4-order Butterworth low-pass filter with cut-off frequencyω
c
=
0.5
π.SetT
=
1 s; then the cut-off frequency of the corresponding analog filter Ω
c can be obtained [4, 5]:
(45)
Ω
c
=
2
T
tan
ω
c
T
2
=
2
rad
/
s
.According to 4-order Butterworth normalized system function [2],
(46)
H
n
(
u
)
=
1
1
+
2.613
u
+
3.414
u
2
+
2.613
u
3
+
u
4
.Takeu
=
s
/
Ω
c into H
n
(
u
); thus
(47)
H
(
s
)
=
16
16
+
20.904
s
+
13.656
s
2
+
5.226
s
3
+
s
4
=
[
1
s
s
2
s
3
s
4
]
a
→
[
1
s
s
2
s
3
s
4
]
b
→
.
So a
→
=
[
16
0
0
0
0
]
T and b
→
=
[
16
20.904
13.656
5.226
1
]
T; according to (41),
(48)
α
→
=
D
4
H
4
a
→
,
β
→
=
D
4
H
4
b
→
,
where, according to (7), (39), and T
=
1,
(49)
D
4
=
[
1
1
1
1
1
-
1
-
2
-
3
-
4
1
3
6
-
1
-
4
1
]
,
H
4
=
[
16
-
32
16
24
-
24
16
-
8
12
-
16
16
1
-
2
4
-
8
16
]
.
For H
4, each column can be calculated by such rules:
(50)
[
16
-
32
24
-
8
1
]
=
D
4
[
1
4
6
4
1
]
,
[
16
-
24
12
-
2
]
=
2
·
D
3
[
1
3
3
1
]
,
[
16
-
16
4
]
=
2
2
·
D
2
[
1
2
1
]
,
[
16
-
8
]
=
2
3
·
D
1
[
1
1
]
,
16
=
2
4
·
D
0
,
where
(51)
D
0
=
1
,
D
1
=
[
1
1
-
1
]
,
D
2
=
[
1
1
1
-
1
-
2
1
]
,
D
3
=
[
1
1
1
1
-
1
-
2
-
3
1
3
-
1
]
.
Take D
4 and H
4, a
→ and b
→ into (47); α
→ and β
→ can be obtained as
(52)
α
→
=
16
[
1
4
6
4
1
]
T
=
[
16
64
96
64
16
]
T
,
β
→
=
[
170.24
0
82.752
0
3.008
]
T
.
Thus,
(53)
H
(
z
)
=
[
1
z
-
1
z
-
2
z
-
3
z
-
4
]
α
→
[
1
z
-
1
z
-
2
z
-
3
z
-
4
]
β
→
=
16
+
64
z
-
1
+
96
z
-
2
+
64
z
-
3
+
16
z
-
4
170.24
+
82.752
z
-
2
+
3.008
z
-
4
=
0.0940
+
0.376
z
-
1
+
0.564
z
-
2
+
0.376
z
-
3
+
0.0940
z
-
4
1
+
0.486
z
-
2
+
0.0177
z
-
4
.H
(
z
) could be expressed by difference equation:
(54)
y
n
+
0.486
y
n
-
2
+
0.0177
y
n
-
4
=
0.0940
x
n
+
0.376
x
n
-
1
+
0.564
x
n
-
2
+
0.376
x
n
-
3
+
0.0940
x
n
-
4
.The amplitude-frequency response of Butterworth digital low-pass filter obtained by Tustin transform is shown in Figure3.Simulation of amplitude frequency response.
(a)
The amplitude-frequency response of normalized analog low-pass filter
(b)
The amplitude-frequency response of digital low-pass filter


### 3.2. An Example of ComputingZ-Domain Transfer Function

Reference [15] mentioned design of switch capacitor circuit by Tustin transform. The S-domain transfer function of first-order and second-order are as follows:
(55)
H
1
(
s
)
=
(
±
)
A
0
+
A
1
s
B
0
+
s
,
H
2
(
s
)
=
(
±
)
A
0
+
A
1
s
+
A
2
s
2
B
0
+
B
1
s
+
s
2
,
where A
0, A
1, A
2, B
0, and B
1 are real.In order to getZ-domain transfer function, operate H
1
(
s
) by Tustin transformation as
(56)
α
→
1
=
D
1
H
1
a
→
1
=
[
1
1
0
-
1
]
[
2
0
-
1
2
T
]
[
A
0
A
1
]

=

[
A
0
+
2
A
1
T
A
0
-
2
A
1
T
]
,
β
→
1
=
D
1
H
1
b
→
1
=
[
1
1
0
-
1
]
[
2
0
-
1
2
T
]
[
B
0
1
]
=
[
B
0
+
2
T
B
0
-
2
T
]
,
(57)
H
1
(
z
)
=
(
±
)
A
0
+
2
A
1
/
T
+
(
A
0
-
2
A
1
/
T
)
z
-
1
B
0
+
2
/
T
+
(
B
0
-
2
/
T
)
z
-
1
;then operateH
2
(
s
) by Tustin transformation as
(58)
α
→
2
=
D
2
H
2
a
→
2
=
[
1
1
1
0
-
1
-
2
0
0
1
]
[
4
0
0
-
4
4
T
0
1
-
2
T
4
T
2
]
[
A
0
A
1
A
2
]
=
[
A
0
+
2
A
1
T
+
4
A
2
T
2
2
A
0
-
8
A
2
T
2
A
0
-
2
A
1
T
+
4
A
2
T
2
]
,
β
→
2
=
D
2
H
2
b
→
2
=
[
1
1
1
0
-
1
-
2
0
0
1
]
[
4
0
0
-
4
4
T
0
1
-
2
T
4
T
2
]
[
B
0
B
1
1
]
=
[
B
0
+
2
B
1
T
+
4
T
2
2
B
0
-
8
T
2
B
0
-
2
B
1
T
+
4
T
2
]
,
(59)
H
2
(
z
)
=
(
±
)
(
A
0
+
2
A
1
T
+
4
A
2
T
2
+
(
2
A
0
-
8
A
2
T
2
)
z
-
1
00000000
+
(
A
0
-
2
A
1
T
+
4
A
2
T
2
)
z
-
2
)
×
(
B
0
+
2
B
1
T
+
4
T
2
+
(
2
B
0
-
8
T
2
)
z
-
1
00000
+
(
B
0
-
2
B
1
T
+
4
T
2
)
z
-
2
)
-
1
.Thus,H
1
(
z
) and H
2
(
z
) are the Z-domain transfer functions.


### 3.3. An Example of Extensional Application

Reference [9] presented a matrix method for a biquadratic transformation. By the algorithm recommended in this paper, we get the matrix with little deformation.For annth order discrete-time polynomial
(60)
D
(
z
)
=
∑
i
=
0
n
d
i
z
i
=
[
1
z
⋯
z
n
]
d
→
,
where d
→ is an n-dimensional column vector, d
→
=
[
d
0
d
1
⋯
d
n
]
T.Take the transformation form (61) into D
(
z
):
(61)
z
=
s
2
+
s
+
1
s
2
-
s
+
1
;
then D
(
z
) transformed into a continuous-time polynomial as
(62)
C
(
s
)
=
∑
i
=
0
2
n
c
i
s
i
=
[
1
s
⋯
s
2
n
]
c
→
,
where c
→ is an 2n-dimensional column vector, c
→
=
[
c
0
c
1
⋯
c
2
n
]
T.The relationship between coefficients in (60) and (62) can be expressed as
(63)
c
→
=
Q
n
d
→
,
where Q
n is the (
2
n
+
1
)
×
(
n
+
1
) transformation matrix. Q
n can be obtained by this method.Equal to (61), the transformation form can be deformed as
(64)
z
=
s
2
+
s
+
1
s
2
-
s
+
1
=
1
+
(
s
/
(
s
2
+
1
)
)
1
-
(
s
/
(
s
2
+
1
)
)
→
u
=
-
(
s
/
(
s
2
+
1
)
)
z
=
1
-
u
1
+
u
.Take (64) into D
(
z
), then D
(
z
) transforms into F
(
u
) as follows:
(65)
D
(
z
)
=
[
1
z
⋯
z
n
]
d
→
→
z
=
(
1
-
u
)
/
(
1
+
u
)
F
(
u
)
=
[
1
u
⋯
u
n
]
f
→
.According to (44), the relation between d
→ and f
→ can be expressed as (66); here consider the step as T
=
2:
(66)
f
→
=
D
n
H
n
d
→
.Last, take the equation betweenu and s into C
(
s
), as follows:
(67)
F
(
u
)
=
[
1
u
⋯
u
n
]
f
→
→
u
=
-
(
s
/
(
s
2
+
1
)
)
C
(
s
)
=
[
1
s
⋯
s
2
n
]
c
→
.Similar to the deduction of (44), the relation between f
→ and c
→ can be expressed as
(68)
c
→
=
G
n
f
→
,
where G
n is a (
2
n
+
1
)
×
(
n
+
1
) matrix:
(69)
G
n
=
[
C
n
0
0
0
⋯
0
0
-
C
n
-
1
0
0
⋯
⋮
C
n
1
0
C
n
-
2
0
⋯
0
⋮
⋮
⋮
⋯
(
-
1
)
n
C
0
0
C
n
n
-
1
0
C
n
-
2
n
-
2
⋯
0
0
-
C
n
-
1
n
-
1
0
⋯
⋮
C
n
n
0
0
⋯
0
]
=
[
g
→
n
-
g
→
n
-
1
⋯
(
-
1
)
i
g
→
n
-
i
⋯
(
-
1
)
n
g
→
0
]
g
→
2
k
=
[
⋯
︸
n
-
3
0
C
2
k
k
-
1
0
C
2
k
k
0
C
2
k
k
+
1
0
⋯
︸
n
-
3
]
T
g
→
2
k
+
1
=
[
⋯
︸
n
-
3
C
2
k
+
1
k
-
1
0
C
2
k
+
1
k
0
C
2
k
+
1
k
+
1
0
C
2
k
+
1
k
+
2
⋯
︸
n
-
3
]
T
00000000000000000000000000000000000
k
=
0,1
,
2
,
…
.Finally,Q
n can be expressed as
(70)
Q
n
=
G
n
D
n
H
nwhilen
=
2:
(71)
G
2
=
[
1
0
0
0
-
1
0
2
0
1
0
-
1
0
1
0
0
]
,
D
2
=
[
1
1
1
0
-
1
-
2
0
0
1
]
,
H
2
=
[
4
0
0
-
4
2
0
1
-
1
1
]
(
T
=
2
)
;
thus Q
2 can be obtained by (70), which is the same as shown in [8]:
(72)
Q
2
=
G
2
D
2
H
2
=
[
1
1
1
-
2
0
2
3
1
3
-
2
0
2
1
1
1
]
.


## 3.1. An Example of Digital Filter Design

For example, the design of a 4-order Butterworth low-pass filter with cut-off frequencyω
c
=
0.5
π.SetT
=
1 s; then the cut-off frequency of the corresponding analog filter Ω
c can be obtained [4, 5]:
(45)
Ω
c
=
2
T
tan
ω
c
T
2
=
2
rad
/
s
.According to 4-order Butterworth normalized system function [2],
(46)
H
n
(
u
)
=
1
1
+
2.613
u
+
3.414
u
2
+
2.613
u
3
+
u
4
.Takeu
=
s
/
Ω
c into H
n
(
u
); thus
(47)
H
(
s
)
=
16
16
+
20.904
s
+
13.656
s
2
+
5.226
s
3
+
s
4
=
[
1
s
s
2
s
3
s
4
]
a
→
[
1
s
s
2
s
3
s
4
]
b
→
.
So a
→
=
[
16
0
0
0
0
]
T and b
→
=
[
16
20.904
13.656
5.226
1
]
T; according to (41),
(48)
α
→
=
D
4
H
4
a
→
,
β
→
=
D
4
H
4
b
→
,
where, according to (7), (39), and T
=
1,
(49)
D
4
=
[
1
1
1
1
1
-
1
-
2
-
3
-
4
1
3
6
-
1
-
4
1
]
,
H
4
=
[
16
-
32
16
24
-
24
16
-
8
12
-
16
16
1
-
2
4
-
8
16
]
.
For H
4, each column can be calculated by such rules:
(50)
[
16
-
32
24
-
8
1
]
=
D
4
[
1
4
6
4
1
]
,
[
16
-
24
12
-
2
]
=
2
·
D
3
[
1
3
3
1
]
,
[
16
-
16
4
]
=
2
2
·
D
2
[
1
2
1
]
,
[
16
-
8
]
=
2
3
·
D
1
[
1
1
]
,
16
=
2
4
·
D
0
,
where
(51)
D
0
=
1
,
D
1
=
[
1
1
-
1
]
,
D
2
=
[
1
1
1
-
1
-
2
1
]
,
D
3
=
[
1
1
1
1
-
1
-
2
-
3
1
3
-
1
]
.
Take D
4 and H
4, a
→ and b
→ into (47); α
→ and β
→ can be obtained as
(52)
α
→
=
16
[
1
4
6
4
1
]
T
=
[
16
64
96
64
16
]
T
,
β
→
=
[
170.24
0
82.752
0
3.008
]
T
.
Thus,
(53)
H
(
z
)
=
[
1
z
-
1
z
-
2
z
-
3
z
-
4
]
α
→
[
1
z
-
1
z
-
2
z
-
3
z
-
4
]
β
→
=
16
+
64
z
-
1
+
96
z
-
2
+
64
z
-
3
+
16
z
-
4
170.24
+
82.752
z
-
2
+
3.008
z
-
4
=
0.0940
+
0.376
z
-
1
+
0.564
z
-
2
+
0.376
z
-
3
+
0.0940
z
-
4
1
+
0.486
z
-
2
+
0.0177
z
-
4
.H
(
z
) could be expressed by difference equation:
(54)
y
n
+
0.486
y
n
-
2
+
0.0177
y
n
-
4
=
0.0940
x
n
+
0.376
x
n
-
1
+
0.564
x
n
-
2
+
0.376
x
n
-
3
+
0.0940
x
n
-
4
.The amplitude-frequency response of Butterworth digital low-pass filter obtained by Tustin transform is shown in Figure3.Simulation of amplitude frequency response.
(a)
The amplitude-frequency response of normalized analog low-pass filter
(b)
The amplitude-frequency response of digital low-pass filter


## 3.2. An Example of ComputingZ-Domain Transfer Function

Reference [15] mentioned design of switch capacitor circuit by Tustin transform. The S-domain transfer function of first-order and second-order are as follows:
(55)
H
1
(
s
)
=
(
±
)
A
0
+
A
1
s
B
0
+
s
,
H
2
(
s
)
=
(
±
)
A
0
+
A
1
s
+
A
2
s
2
B
0
+
B
1
s
+
s
2
,
where A
0, A
1, A
2, B
0, and B
1 are real.In order to getZ-domain transfer function, operate H
1
(
s
) by Tustin transformation as
(56)
α
→
1
=
D
1
H
1
a
→
1
=
[
1
1
0
-
1
]
[
2
0
-
1
2
T
]
[
A
0
A
1
]

=

[
A
0
+
2
A
1
T
A
0
-
2
A
1
T
]
,
β
→
1
=
D
1
H
1
b
→
1
=
[
1
1
0
-
1
]
[
2
0
-
1
2
T
]
[
B
0
1
]
=
[
B
0
+
2
T
B
0
-
2
T
]
,
(57)
H
1
(
z
)
=
(
±
)
A
0
+
2
A
1
/
T
+
(
A
0
-
2
A
1
/
T
)
z
-
1
B
0
+
2
/
T
+
(
B
0
-
2
/
T
)
z
-
1
;then operateH
2
(
s
) by Tustin transformation as
(58)
α
→
2
=
D
2
H
2
a
→
2
=
[
1
1
1
0
-
1
-
2
0
0
1
]
[
4
0
0
-
4
4
T
0
1
-
2
T
4
T
2
]
[
A
0
A
1
A
2
]
=
[
A
0
+
2
A
1
T
+
4
A
2
T
2
2
A
0
-
8
A
2
T
2
A
0
-
2
A
1
T
+
4
A
2
T
2
]
,
β
→
2
=
D
2
H
2
b
→
2
=
[
1
1
1
0
-
1
-
2
0
0
1
]
[
4
0
0
-
4
4
T
0
1
-
2
T
4
T
2
]
[
B
0
B
1
1
]
=
[
B
0
+
2
B
1
T
+
4
T
2
2
B
0
-
8
T
2
B
0
-
2
B
1
T
+
4
T
2
]
,
(59)
H
2
(
z
)
=
(
±
)
(
A
0
+
2
A
1
T
+
4
A
2
T
2
+
(
2
A
0
-
8
A
2
T
2
)
z
-
1
00000000
+
(
A
0
-
2
A
1
T
+
4
A
2
T
2
)
z
-
2
)
×
(
B
0
+
2
B
1
T
+
4
T
2
+
(
2
B
0
-
8
T
2
)
z
-
1
00000
+
(
B
0
-
2
B
1
T
+
4
T
2
)
z
-
2
)
-
1
.Thus,H
1
(
z
) and H
2
(
z
) are the Z-domain transfer functions.


## 3.3. An Example of Extensional Application

Reference [9] presented a matrix method for a biquadratic transformation. By the algorithm recommended in this paper, we get the matrix with little deformation.For annth order discrete-time polynomial
(60)
D
(
z
)
=
∑
i
=
0
n
d
i
z
i
=
[
1
z
⋯
z
n
]
d
→
,
where d
→ is an n-dimensional column vector, d
→
=
[
d
0
d
1
⋯
d
n
]
T.Take the transformation form (61) into D
(
z
):
(61)
z
=
s
2
+
s
+
1
s
2
-
s
+
1
;
then D
(
z
) transformed into a continuous-time polynomial as
(62)
C
(
s
)
=
∑
i
=
0
2
n
c
i
s
i
=
[
1
s
⋯
s
2
n
]
c
→
,
where c
→ is an 2n-dimensional column vector, c
→
=
[
c
0
c
1
⋯
c
2
n
]
T.The relationship between coefficients in (60) and (62) can be expressed as
(63)
c
→
=
Q
n
d
→
,
where Q
n is the (
2
n
+
1
)
×
(
n
+
1
) transformation matrix. Q
n can be obtained by this method.Equal to (61), the transformation form can be deformed as
(64)
z
=
s
2
+
s
+
1
s
2
-
s
+
1
=
1
+
(
s
/
(
s
2
+
1
)
)
1
-
(
s
/
(
s
2
+
1
)
)
→
u
=
-
(
s
/
(
s
2
+
1
)
)
z
=
1
-
u
1
+
u
.Take (64) into D
(
z
), then D
(
z
) transforms into F
(
u
) as follows:
(65)
D
(
z
)
=
[
1
z
⋯
z
n
]
d
→
→
z
=
(
1
-
u
)
/
(
1
+
u
)
F
(
u
)
=
[
1
u
⋯
u
n
]
f
→
.According to (44), the relation between d
→ and f
→ can be expressed as (66); here consider the step as T
=
2:
(66)
f
→
=
D
n
H
n
d
→
.Last, take the equation betweenu and s into C
(
s
), as follows:
(67)
F
(
u
)
=
[
1
u
⋯
u
n
]
f
→
→
u
=
-
(
s
/
(
s
2
+
1
)
)
C
(
s
)
=
[
1
s
⋯
s
2
n
]
c
→
.Similar to the deduction of (44), the relation between f
→ and c
→ can be expressed as
(68)
c
→
=
G
n
f
→
,
where G
n is a (
2
n
+
1
)
×
(
n
+
1
) matrix:
(69)
G
n
=
[
C
n
0
0
0
⋯
0
0
-
C
n
-
1
0
0
⋯
⋮
C
n
1
0
C
n
-
2
0
⋯
0
⋮
⋮
⋮
⋯
(
-
1
)
n
C
0
0
C
n
n
-
1
0
C
n
-
2
n
-
2
⋯
0
0
-
C
n
-
1
n
-
1
0
⋯
⋮
C
n
n
0
0
⋯
0
]
=
[
g
→
n
-
g
→
n
-
1
⋯
(
-
1
)
i
g
→
n
-
i
⋯
(
-
1
)
n
g
→
0
]
g
→
2
k
=
[
⋯
︸
n
-
3
0
C
2
k
k
-
1
0
C
2
k
k
0
C
2
k
k
+
1
0
⋯
︸
n
-
3
]
T
g
→
2
k
+
1
=
[
⋯
︸
n
-
3
C
2
k
+
1
k
-
1
0
C
2
k
+
1
k
0
C
2
k
+
1
k
+
1
0
C
2
k
+
1
k
+
2
⋯
︸
n
-
3
]
T
00000000000000000000000000000000000
k
=
0,1
,
2
,
…
.Finally,Q
n can be expressed as
(70)
Q
n
=
G
n
D
n
H
nwhilen
=
2:
(71)
G
2
=
[
1
0
0
0
-
1
0
2
0
1
0
-
1
0
1
0
0
]
,
D
2
=
[
1
1
1
0
-
1
-
2
0
0
1
]
,
H
2
=
[
4
0
0
-
4
2
0
1
-
1
1
]
(
T
=
2
)
;
thus Q
2 can be obtained by (70), which is the same as shown in [8]:
(72)
Q
2
=
G
2
D
2
H
2
=
[
1
1
1
-
2
0
2
3
1
3
-
2
0
2
1
1
1
]
.


## 4. Conclusions

The method of using corresponding matrix for Tustin transform representation reflects the inherent properties; on the other hand, it allows computer programs to deal with the transfer functions of ordern automatically. In this paper, the transformation matrix is obtained by multiplying D
n and H
n, which dramatically reduces the complexity of the program without any tabulation. Furthermore, this method can also be modified properly to suit similar problems about coefficients after linear transform.



---

*Source: 905791-2014-08-18.xml*