Design of Systems with Prescribed Impulse Response. Based on Second-Order Cone Programming. Mladen Vucic and Goran Molnar. Faculty of Electrical ...
Design of Systems with Prescribed Impulse Response Based on Second-Order Cone Programming Mladen Vucic and Goran Molnar Faculty of Electrical Engineering and Computing University of Zagreb Unska 3, Zagreb, HR10000, Croatia
AbstractA method for the design of continuous-time systems approximating prescribed impulse response is presented. The design uses weighted least squares criterion, which is applied in the time domain. The method is developed based on iterative procedure solving a second-order cone programming (SOCP) subproblem in each iteration. It supports design criteria such as stability and minimum phase. Features of the proposed method are illustrated with a design of pulse-shaping filter for ultra-wideband (UWB) transmission systems.
I. INTRODUCTION Time-domain synthesis is used when a network is required which gives a prescribed time-domain response for a given excitation [1]. Such networks often appear in communication systems [2] and measuring equipment [3]. There are numerous methods developed for their design. A good overview of classic methods is given in [1]. In recent years, new methods have been developed for synthesis of networks approximating a particular impulse response [4], [5], [6], [7] and a particular step response shape [7], [8]. These methods have been used in the synthesis of pulse shapers for ultra-wideband (UWB) transmission, as illustrated in [4] and [5], and in the synthesis of high voltage pulse shapers, as shown in [7] and [8]. Here we present a new method for time-domain synthesis of continuous-time systems. The method is based on the second-order cone programming (SOCP), an optimization technique which has been popular in recent years. Second-order cone program is a convex optimization problem in which a linear function is minimized subject to a set of second-order cone constraints [9], [10], as in min x
subject to :
fTx
Ai ⋅ x + bi ≤ ciT x + d i , i = 1,..., m
(1) (2)
Many problems in signal processing can be expressed in a form of second-order cone program, [10]. On the other hand, there are problems that cannot be recognized as convex optimization problems at all. However, a method for solving these problems have been developed, which is based on sequence of linear updates, where each update is obtained by solving a SOCP subproblem. This method has been applied in the frequencydomain synthesis of recursive digital systems, as for example in [13] and [14]. The problem we are addressing here is the time-domain synthesis of continuous-time systems having prescribed impulse response. We will present a method for time-domain synthesis, which is based on iterative procedure containing a SOCP subproblem in each iteration. Furthermore, we will apply this method in a design of pulse-shaping filter for ultra-wideband transmission.
II. DESIGN METHOD A. Problem Formulation Our goal is to find a system whose impulse response minimizes the integral of the weighted squared error given by ∞
e wls ( x ) = w(t )[h(t , x ) − hd (t )]2 dt
∫
where h(t,x) is system's impulse response, x is the vector of system parameters, hd(t) is a desired impulse response and w(t) is a positive weighting function. We start our design by the choice of vector x, i.e. by the representation of the system being optimized. The transfer function of an Nth order system with M zeros is given by
is the optimization variable, and f ∈ ℜ n , Ai ∈ ℜ bi ∈ ℜ n i −1 , c i ∈ ℜ n , d i ∈ ℜ are the problem parameters. The norm appearing in (2) is the standard Euclidian norm, given by u = (u T u)1 2 . A constraint in (2) is called a second-order cone constraint because it is equivalent to
where
x ∈ ℜn ( n i −1) × n ,
Ai bi T x + ∈ Ci c d i i
(3)
where Ci is a second-order cone of dimension ni, given by u Ci = t
u ∈ ℜ ni −1 , t ∈ ℜ, u ≤ t
(4)
Efficient interior-point methods for solving problem in (1) and (2) have been described in literature, as for example in [11]. The appropriate software is also available, [12].
1-4244-0921-7/07 $25.00 © 2007 IEEE.
(5)
0
M
∏ (s − zi )
H ( s ) = H 0 i =1 N
∏
(6)
(s − pk )
k =1
where pk, zi and H0 denote system's poles, zeros and gain constant. Since the system is uniquely described by pk, zi and H0, these parameters are good candidates for the components of x. However, to form x as a real vector, we described complex pairs of the poles and the zeros by their real and imaginary parts. If the transfer function contains M1 real zeros, M2 complex zeros, N1 real poles and N2 complex poles, where M1+M2=M and N1+N2=N, x can be defined as
3311
x=
[H0,
z1 ,..., z M 1 , α1 , β1 ,..., α M 2 / 2 , β M 2 / 2 ,
]
p1 ,..., p N1 , σ 1 , ω1 ,..., σ N 2 / 2 , ω N 2 / 2 T
x
x
z M 1 + i = α i + jβ i , z ∗M 1 + i = α i − jβ i , i = 1,..., M 2 / 2
(8)
p N1 + k = σ k + jω k , p∗N1 + k = σ k − jωk , k = 1,..., N 2 / 2
(9)
Ts A ⋅ δ + b
min δ
By using η as a slack variable, problem in (17) and (18) can be written in a form min η
To solve the problem in (10), we formed an iterative procedure. It starts from an initial point x0. Based on the behavior of ewls(x) around x0, the appropriate correction δ0 is found and new point is obtained as x1=x0+δ0. This procedure is then repeated until the termination criterion is satisfied. In each iteration, the correction is obtained by solving an appropriate SOCP subproblem. Suppose that, after k iterations, we have arrived to the point xk and now we are searching for δk. If h(t,x) is a smooth function, in the vicinity of xk it can be approximated by a linear form, as in h(t , x k + δ) ≈ h(t , x k ) + g T (t , x k ) δ
(11)
where g(t,xk) is the gradient of h(t,x) with respect to x and evaluated at the point xk. Note that the approximation error is small within some close neighborhood, V, of the point xk, i.e. for δ∈V. Using (11), ewls(x) can be approximated by
}2
(t , x k ) δ + w(t ) [h(t , x k ) − hd (t )] dt (12)
0
In practice, h(t) and hd(t) can be considered as band-limited. Therefore, they can be represented by sets of discrete equidistant samples h(tq) and hd(tq), where tq=qTs, q=0,1,...,Q, and Ts is a sampling period which is small enough to ensure negligible aliasing error. Assuming h(qTs)≈0 for q>Q, integral in (12) can be approximated by a sum, and the error ewls(x) can be written as e wls ( x ) ≈ Ts
∑ { w(t q ) g T (t q , x k )δ +
Ts η
(20)
The above problem is recognized as a SOCP problem
B. Optimization Procedure
Q
(19)
A⋅ δ + b ≤ 1 δ ∈V
subject to :
(10)
x
T
(18)
η, δ
x opt = arg min[e wls ( x )]
∫ { w(t ) g
(16)
(17)
subject to : δ ∈ V
Using x in (7), the optimum system is found as
∞
]
δk can also be found by solving the problem
where αi, σk, βi and ωk are real and imaginary parts of complex zeros and poles, as given in
e wls ( x ) ≈
[
arg min[ewls ( x )] = arg min e wls ( x )
(7)
[
w(t q ) h(t q , x k ) − hd (t q )
]}
2
Expression (13) can be written in a matrix form as 2
u
subject to : Ai ⋅ u + bi ≤ c iT u + d i , i = 1,2
(22)
where the objective function, f, and the optimization variable, u, are given by f = [1 0 L 0] T ,
η u= δ
(23)
The first constraint in (22) contains the optimization problem itself. Its parameters are given by A1 = [0
[
A] , b1 = b , c1 = 1 Ts 0
L 0
] T , d1 = 0
(24)
The second constraint in (22) defines the neighborhood, V, within which the linear approximation is used. A possible choice for V is a ball, given by ||δ||≤δmax. Such a choice is appropriate for bounding movement of the poles and the zeros since they are grouped in the complex plane. On the other hand, the gain constant, H0, can take values that are several orders of magnitude higher and so cause a slow convergence rate. To avoid it, we formed the second constraint by using relative bound for H0, as in 0 1 / a A2 = 0 0 M M 0 0
0 1 M 0
K L O L
0 0 , , , M b2 = 0 c 2 = 0 d 2 = δ max 1
(25)
where a is the first component of xk. To form a SOCP subproblem, system's impulse response h(t,x) should be known. Using (8) and (9), x can be easily transformed into zeros and poles. If the poles are simple and M0. Formally, it can be expressed by a set of linear inequalities
In each iteration of the optimization procedure, a SOCP problem is solved, which is given by (19), (20) and (38). It is clear that the problem contains two second-order cone constraints, given in (20), and N1+N2/2 linear constraints, given in (38). Such a problem can be solved using an interior-point method for the second-order cone programming described for example in [10] and [11]. We have used the toolbox SeDuMi 1.05, which can be found in [12]. In the beginning of the optimization procedure, i.e. for small k, the magnitude of δk is usually determined by the second second-order cone constraint, ||A2u||≤δmax. When xk approaches to xopt, the magnitude starts with rapid decreasing. The iteration procedure finishes when ||A2u||
