Simulation of Wave Propagation by Multidimensional Digital Filters
Simulation of Wave Propagation by Multidimensional Digital Filters H. Krau, R. Rabenstein, M. Gerken Lehrstuhl fur Nachrichtentechnik, Universitat Erlangen{Nurnberg, Cauerstr. 7, D-91058 Erlangen, Germany Tel.: +9131/85-7101 Fax: +9131/303840 Email: [email protected]
Abbreviated title: Simulation of Wave Propagation Abstract
The application of multidimensional digital ltering theory to the simulation of continuous systems with distributed parameters has been treated in various recent contributions. Here, two di erent concepts will be compared, using a lossy transmission line of nite extent as an example. One approach is based on multidimensional wave digital lters, the other on the derivation of a simulating structure by functional transformations for each independent variable. The comparison includes a review of the underlying theoretical concepts as well as the presentation of the design process for the digital system and simulation results for various examples.
Keywords:
Computer Simulation, Wave Propagation, Digital Signal Processing, Digital Filters, Wave Digital Filters
Escola Politecnica da USP, S~ao Paulo, Brazil, worked with a post-doc fellowship at the University of Erlangen{Nurnberg, Germany.
1
1 Introduction The simulation of lumped parameter systems by onedimensional digital signal processing methods is well known. It has been used also for the simulation of distributed parameter systems like transmission lines by applying a spatial discretization to the describing partial di erential equation. This procedure is known as the method of lines and gives a sequence of interconnected lumped parameter systems. Each of these can be replaced by an appropriate onedimensional digital lter for simulation purposes. Multidimensional digital signal processing methods have been con ned to image and seismic processing for a long time. Recently, the application of the theory of multidimensional digital lters to the numerical solution of partial di erential equations and to the simulation of continuous systems with distributed parameters has been treated in the literature [4, 5, 12, 17, 18]. These methods allow the replacement of the distributed parameter system directly by a multidimensional digital lter, without the need for an initial conversion of the continuous system into a series of lumped parameter systems. Two of these methods will be compared here with respect to the simulation of a lossy transmission line. The transmission line problem has been chosen as an example of a hyperbolic partial di erential equation since both methods can be applied in a more general way to a variety of wave propagation phenomena. What is more, the method of section 3 which is based on a functional transformation approach is applicable to partial di erential equations of the parabolic type. Thus it it is quali ed e.g. for the simulation of heat and mass transfer problems [17, 18, 20]. The method based on multidimensional wave digital lter principles can be extended to propagation problems in more than one spatial dimension [4, 5, 6] and to more elaborate boundary value problems [11]. Section 2 gives a mathematical description of the transmission line problem. The digital ltering concepts to be compared are presented in sections 3 and 4. The examples for the numerical experiments and the results of the comparison are given in section 5.
2
2 Problem Statement We consider a transmission line of length A described by the two partial di erential equations @ u(x; t) + l @ i(x; t) + ri(x; t) = 0 @x @t (1) @ i(x; t) + c @ u(x; t) + gu(x; t) = 0 @x @t for voltage u(x; t) and current i(x; t) along the line. Both voltage and current depend on the space coordinate x with 0 x A and on the time coordinate t with 0 t. The coecients are given in terms of the following physical quantities: l : serial inductance in H=m, c : shunt capacitance in F=m, r : serial resistance in =m, g : shunt conductance in S=m. The two partial di erential equations (1) of rst order may be transformed into the second-order partial di erential equation @ y(x; t) = lc @ y(x; t) + (lg + rc) @ y(x; t) + rgy(x; t) ; (2) @x @t @t the so-called telegraph equation, where y(x; t) denotes voltage u(x; t) or current i(x; t). Unless otherwise stated y(x; t) in the sequel stands for the voltage u(x; t). The line shall be initially discharged, that is y(x; 0) = 0 ; (3) @ y(x; t) = 0: @t t At the ends of the line, we prescribe boundary conditions of the rst kind y(0; t) = yb (t) ; (4) y(A; t) = yb (t) or of the second kind @ y(x; t) = yb (t) ; @x x (5) @ y(x; t) = yb (t) : @x x A The transmission line is said to be free of distortions if r =g (6) l c holds. The following sections describe two multidimensional digital systems for the simulation of this continuous system. 2
2
2
2
=0
1
2
1
=0
2
=
3
3 Functional Transformation Approach (FTA) This section describes the specialization of a general method for the discrete simulation of linear multidimensional continuous systems to the problem stated previously. A full account of the method and an application to the heat ow equation is given in [16, 17, 18]. The application to wave propagation problems is given in [19]. The key steps are the representation of the continuous system by a transfer function, the derivation of the transfer function of a corresponding discrete system and nally the speci cation of a realizing structure, which can be easily cast into a computer program. This procedure is well known from the analysis of onedimensional systems. For time dependent problems like e.g. electrical networks, the Laplace{transformation [14, 21] provides an elegant tool to represent an ordinary di erential equation along with its initial condition(s) by appropriate transfer functions. This method works well because the integration range of the Laplace{transformation matches the region on which the inital value problem is de ned (0 t < 1 in both cases). An extension of this method to the space{ and time{dependent initial{boundary value problem described in section 2 calls for a further functional transformation for the space variable with 0 x A. Due to reasons explained in [18], it is of advantage to use an odd (even) continuation of the spatial functions into the interval ?A < x < 0 for boundary conditions of the rst (second) kind and a periodic repetition with period 2A outside of ?A x A. Then the nite Fourier{transformation with an integration range of ?A x A will include boundary conditions with respect to x in the same way into a transfer function as the Laplace{transformation does for the initial conditions with respect to t. Mixed type boundary conditions, i.e. boundary conditions of the rst kind on the one end and of the second kind at the other end, may be accomplished using a periodic repetition with period 4A. See [18, 21] for a discussion of the nite Fourier{ transformation and its properties. A de nition will be given later in this section. The procedure will now be demonstrated for the telegraph equation (2).
3.1 Continuous System
In order to obtain a transfer function from the partial di erential equation representation (2) of the transmission line, we at rst apply the one{sided Laplace{transformation with respect to the time coordinate t 1
Z
Lfy(x; t)g = Y (x; s) = y(x; t)e?stdt ; s 2 C :
(7)
I
0
Since the in uence of the boundary conditions is a linear superposition of yb (t) and yb (t), we can assume yb (t) = 0 for simplicity and apply the Laplace{transformation only to the non{vanishing boundary value at x = 0. The result is @ Y (x; s) = (s) Y (x; s) ; Y (s) = Lfy (t)g ; (8) b b @x with the propagation function 1
2
2
2
2
1
2
q
(s) = (sl + r)(sc + g) : 4
1
(9)
For later convenience, we introduce the phase velocity v and the abbreviations and by (10) v = p1 ; = ? rl ; = ? gc : lc The squared propagation function is thus (11)
(s) = v1 (s ? )(s ? ) : The next step is the application of the nite Fourier{transformation ( n. FT) 0
1
2
1
0
2
2
1
2 0
FAfY (x; s)g = Y (; s) =
A
Z
2
Y (x; s)e?j =A x dx ; 2 ZZ (
?A
)
(12)
1 1 FA = Y (x; s) = 2A Y (; s)ej =A x (13) ?1 in space direction. For boundary conditions of the rst kind, we have to assume an odd continuation and periodic repetition in the sense stated above. With the derivation theorem of the nite Fourier{transformation [18, 21], we obtain ? Y (; s) = (s) Y (; s) + 2j Yb (s) (14) with (15) = A : This algebraic equation can be solved for Y (; s) to obtain Y (; s) = G b (; s) Yb (s) (16) with G b (; s) = ? (2sj : (17) ) + Note, that equations (16,17) are an alternate description of the continuous system given by (2,3,4). It is not based on the space and time coordinates x and t, but on the corresponding transformed variables and s, respectively. We can envisage G b (; s) in (17) as the twodimensional transfer function between the input Yb (s) and the output Y (; s). For a given transfer function, (16) is a concise description of the problem previously given by the partial di erential equation (2) together with initial (3) and boundary conditions (4). Using (10,11,15), we can write the denominator of the transfer function (17) as (18)
(s) + = v1 (s ? )(s ? ) + ! where ! = A v = 2T : (19) The quantity p (20) T = 2vA = 2A lc ?1 fY (; s)g
X
(
)
=
2
2
1
1
1
1
2
2
1
2
1
1
2
h
2
2 0
0
0
0
5
0
2
2
i
can be regarded as the time it takes for a wave with phase velocity v to travel across an (extended) line of length 2A. Equation (18) is a polynomial in s with the zeros 0
s
s1 = ( ) = +2 ?2 ? ! ; (21) which constitute the poles of the transfer function (17) v G b (; s) = ?j (s ? s (2! : (22) 1 ))(s ? s1 ( )) We can describe the in uence of a boundary value yb (t) in the time domain by 1
1 2
2
1
2
2
2
2
0
1
1
2
1
y(; t) = gb (; t) t yb (t) (23) with y(; t) = FAfy(x; t)g = L? fY (; s)g and gb (; t) = L? fG b (; s)g = ?j s (2)!? vs ( ) [exp(s1 ( )t) ? exp(s1 ( )t)] ? (t) : 1 1 (24) (? (t) is the unit step function.) By setting yb (t) = (t), we see that gb (; t) describes the reaction of the transmission line to an impulse at the boundary. Equation (24) speci es this reaction as an in nite set of time functions, each valid for a di erent discrete spatial frequency with ?1 < < 1 (see (12,13)). The corresponding description in the space{time domain can be obtained by inverse nite Fourier{transformation (13) as 1 gb (x; t) = 21A gb (; t)ej x : (25) ?1 1
1
1
1
1
0
1
1
1
2
2
1
1
1
0
1
X
1
1
=
3.2 Discrete System
The equations (2,3,4,25), (16,17), and (23,24) give di erent representations of the transmission line in the (x; t), (; s) and (; t) domain, respectively. We will now specify a procedure to obtain a multidimensional digital lter for the simulation of the continuous transfer function G b (; s). Its purpose is to calculate a close approximation to the output y(x; t) of the continuous system at the grid points (x; t) = (nh; kT ), where h and T are the spatial and temporal step sizes respectively, and n and k are the corresponding discrete spatial and temporal variables. The design of the twodimensional digital lter will proceed in two steps. At rst we will discretize the time variable, based on equations (23,24) separately for each spatial frequency . The second step is the discretization of the space variable, based on the inverse nite Fourier{transformation (13). 1
3.2.1 Discretization of the Time Variable
The discretization of the time variable is based on the well known impulse{invariant transformation [9, 14, 22]. To this end, we consider gb (; t) as the impulse response of an 1
6
analog lter with a discrete parameter . Then the impulse response gbh (; k) of a hybrid system that is discrete in time, but continuous in space is 1
gbh (; k) = T gb (; kT ) :
(26)
1
1
From (24), we obtain
gbh (; k) = ?j s (2)!? vs ( ) z1k ( ) ? z1k ( ) ? (k) 1 1 h
0
1
1
i
1
2
2
(27)
1
with
z1i( ) = exp(s1i( )T ) ; i = 1; 2 : (28) Z{transformation [8] gives the transfer function with z G hb (; z) = z + cb(b()z)+ (29) c ( ) ; where bb ( ) = ?2j!v z1 ( ) ? z1 ( ) ; (30) s1 ( ) ? s1 ( ) c ( ) = ?(z1 ( ) + z1 ( )) ; (31) c ( ) = z1 ( ) z1 ( ) : (32) For further evaluation of the coecients bb ( ), c ( ), c ( ), we have to distinguish, whether the zeros of (18), as given by (21) are both real and distinct, equal to each other or a complex conjugate pair. In terms of (18), the cases j j >< 2?! (33) are possible. The results for all three cases can be summarized as bb ( ) = ?2j!v T exp( T ) sin ! T ; (34) ! T c ( ) = ?2 exp( T ) cos ! T ; (35) c ( ) = c = exp(2 T ) : (36) with = 21 ( + ) = ? 21 rl + gc (37) ! = ! ? 2?! : (38) 1
1
2
1
1
0
0
1
1
2
1
2
1
0
2
1
2
1
=
1
0
1
0
0
0
2
1
1
0
2
0
0
0
1
2
s
2
1
2
2
According to the di erent cases given by (33), ! may be real, zero, or imaginary. The special case of a transmission line without distortions (6) is characterized by = = and ! = ! . Using the transfer function G hb (; z) from (29), we can describe the hybrid system by Y h (; z) = G hb (; z)Ybd (z) (39) 1
1
1
7
1
2
0
with the input Ybd (z) = Z fyb (kT )g. Insertion of (29), reordering of terms and inverse z{transformation gives yh(; k) = ?c ( )yh (; k ? 1) ? c yh(; k ? 2) + bb ( )ybd (k ? 1) : (40) 1
1
1
0
1
1
By inverse n. FT, we could obtain
yh(x; k) = ?c (x) x yh(x; k ? 1) ? c yh(x; k ? 2) + bb (x)ybd (k ? 1) ; 1
0
1
1
(41)
i. e. the output of a continuous space, discrete time system with
yh(x; k) y(x; kT ) :
(42)
Since we are not interested in hybrid simulation, we use the hybrid system (39,40,41) as an intermediate stage to a fully discrete system. The remaining step | the discretization of the space variable | will be covered next.
3.2.2 Discretization of the Space Variable
A rst attempt to a fully discrete system is the evaluation of (41) at discrete points in space x = nh. This poses no problem for the second and third term on the right hand side. However, the convolution term cannot be solved exactly on a discrete grid for general c (x) and yh(x; k ? 1). We will therefore focus our attention on a discrete approximation of the term 1
f h (x; k) x
= c (x) x yh(x; k) nh h
1
=
i
n
= FA? c ( )yh(; k) nh 1
x=
1
o
x=nh
:
(43)
An exact solution would always be possible, if yh(x; k) was a bandlimited function in space. This prerequisite is not met, since bb (x) is not bandlimited (see its n. FT in (34)). We may, however, apply the rectangular rule of integration to intervals of length (44) h = NA ; N 2 IN for the convolution (43). This yields the approximate solution 1
1
f d (n; k)
=
f h(x = nh; k) =
h
N X =?N +1
A
Z
?A
c ()yh (nh ? ; k)d 1
c (h)yh((n ? )h; k) = 1
N
X
=?N +1
cd()yd(n ? ; k) ; 1
(45)
where cd(n) and yd(n; k) denote the samples of h c (x) and yh(x; k) at locations x = nh. Choosing the xed relationship h=v T (46) 1
1
0
In case of a transmission line without distortions (see (6)) equation (45) even provides the exact solution. 1
8
between the spatial step size h and the temporal sampling interval T it is shown in [19] that cd(n) and yd(n; k) in (45) may be calculated by means of an inverse Discrete Fourier Transform (DFT) of length 2N . To this aim the Discrete Fourier Transform is applied in the form 1
n
N
o
y^d(; k) = DFT N yd(n; k) = 2
X
n=?N +1 N 1 X
yd(n; k) exp(?j N n ) ;
y^d(; k) exp(j N n ) : yd(n; k) = DFT?N y^d(; k) = 2N ?N Then one gets yd(n; k) = DFT?N y^d(; k) with the DFT coecients y^d(; k) = h1 yh(; k) ; ?N < N : On the other hand cd (n) may be obtained via 2
1
n
o
=
2
+1
n
1
o
(47) (48) (49) (50)
1
cd(n) = DFT?N 1
2
1
8 < X
1
:
i=?1
c ( ? i 2N ) 1
9 =
(51)
;
up to any desired accuracy. Concerning the convergence and the practical evaluation of the in nite series in (51) we refer to [19]. Thus we can write (45) in the form
f d(n; k) =
cd(n)
2N yd(n; k)
1
(52)
with 2N denoting cyclic convolution of length 2N . In summary, the continuous convolution term f h(x; k) from (43) can be evaluated approximately at a spatial grid point (x = nh = nv T ) via a discrete cyclic convolution between the voltage distribution yd(n; k) and a coecient sequence cd(n). The latter can eciently be computed with fast algorithms for the DFT using the spectral components according to (35,50,51). Note, that also the third term on the right hand side of (41) is obtained with the DFT up to any desired accuracy as 0
1
bdb (n) = DFT?N 1
by
2
1
8 < :
1 hi
1
X
?1
=
9 =
bb ( ? i 2N ) : 1
(53)
;
All together, the fully discrete system with the grid points (x; t) = (nh; kT ) is given
yd(n; k) = ? cd(n) 1
2N yd(n; k ? 1) ? c yd(n; k ? 2) + bdb (n) ybd (k ? 1) : 0
1
1
(54)
A realizing structure is shown in g. 1. It is a second order recursive system with respect to the discrete time variable k. In each time step a vector of 2N spatial samples is 9
Figure 1: Structure of the discrete system (see (54)) processed in parallel. The boundary value ybd (k) acts as the input of the system. Its impact on the spatial voltage distribution yd(n; k) is modelled by the sequence bdb (n). Finally we remark that the structure according to g. 1 is not based on a replacement of di erential operators in (2) by di erence operators, as is done in many other numerical methods. Consequently the CFL-condition (cf. section 4.2) governing the stability does notpapply here. This means that the structure is also stable for step sizes T > h=v ; v = 1= lc. A stability criterion for recursive systems of this type can be found in [15]. 1
1
0
10
0
4 Multidimensional Wave Digital Filter (MWDF) In this section a di erent approach based on multidimensional wave digital lter (MWDF) theory will be outlined in brevity. This method has been presented in various publications by Fettweis and Nitsche, e.g. applied to the problem of two conducting plates in [4] and [5]. For a full account of the method we therefore refer to [4, 5, 6], for a survey of MWDF applications to [7]. The key steps of this approach are the derivation of a multidimensional Kirchho circuit from the underlying system of partial di erential equations and its simulation by a multidimensional wave digital lter. Thus a passive continuous system will be simulated by a passive discrete system.
4.1 Derivation of the Kirchho Circuit
In contrast to section 3 now the two partial di erential equations in (1) for voltage u(x; t) and current i(x; t) across a transmission line, having constant parameters r; l; g and c, serve as the basis for the investigations. In order to deal with same quantities a normalization according to i (x; t) := i(x; t) and i (x; t) := u(rx; t) (55) is used where r is an arbitrary positive resistance. One then gets the two equations @ i (x; t) + l @ i (x; t) + r i (x; t) = 0 r @x (56) @t @ i (x; t) + r c @ i (x; t) + r g i (x; t) = 0 (57) @x @t 1
2
2
2
2
2
1
1
2
1
2
2
2
We now go on describing the frequency domain version of the method which is completely equivalent to a direct space-time approach in case of constant parameters r; l; g and c as was shown in [4, 5]. Assuming steady state behaviour one can set i (x; t) = I esT t ; = 1; 2 (58) with complex amplitudes I and vectors sT = [sx; st] ; t = [x; t]T : (59) Inserting (58) and (59) into (56,57) gives the algebraic equations (stl + r)I + sxr I = 0 (60) sxr I + (stcr + gr )I = 0 (61) where the last equation has been multiplied by the factor r . One can easily verify that equations (60) and (61) may be interpreted as mesh equations of a twodimensional Kirchho circuit as depicted in g. 2a. Furthermore g. 2b shows a modi cation of the circuit where an arbitrary positive constant has been introduced. Transformation of the T-con guration of g. 2b into a symmetric lattice structure gives the circuit of g. 2c where the impedances Z ; Z ; Z und Z are given by Z = st(l ? ) ; Z = st(cr ? ) (62) Z = st ? s x r ; Z = st + s x r : (63) 1
2 2
2 2
2 1
2 2
2
2
1
1
3
2
3
2
2
4
11
4
2 2
2
Figure 2: a) Equivalent Kirchho circuit; b) modi ed Kirchho circuit; c) equivalent symmetric lattice structure In order to guarantee passivity of the circuit the following inequalities l ; cr are supposed to hold. With the choice 2 2
(64)
s
= l ; r = cl the impedances Z and Z degenerate to Z = Z = 0: 2
1
(65)
2
1
2
4.2 Derivation of the MWDF
(66)
Similar to the 1D-case the transition from voltage and current quantities u and i to incident and re ected wave quantities a and b, respectively, plays an essential role in the derivation of the 2D WDF. Furthermore the bilinear transformation is applied, this time, however, in a generalized way as the frequency domain equivalent to a generalized trapezoidal rule of integration [4, 5]. 12
In accordance with [4, 5] incident and re ected waves in the frequency domain are introduced as voltage waves
A = U + R I ; B = U ? R I ;
(67)
where the proper orientation of voltage/current and wave quantities can be seen in g. 3. The factor R , called the port resistance, has to be chosen in a suitable way.
Figure 3: Voltage/current and wave quantity orientations of an electrical 1-port The full derivation of the MWDF for the simulation of the Kirchho circuit of g. 2c is beyond the scope of this paper. Here we just give an outline of major steps of this derivation and refer the reader to more profound WDF literature, especially [3] and [4, 5]. As pointed out in [3] the lattice twoport of g. 4a possesses the wave equivalent shown in g. 4b. Note that the twoport of g. 4a constitutes a part of the Kirchho circuit of g. 2c. In order to get the structure of g. 4b, the wave description (67) with orientation
Figure 4: a) Symmetric lattice twoport; b) equivalent wave representation as in g. 3 must be applied to each of the two ports of the lattice. By choosing the two port resistances to be equal to the same positive value R , which is not necessary but convenient for symmetry reasons, the scattering quantities S and S in g. 4b take the 0
1
13
2
form
R ; (68) S = ZZ ? +R R : (69) S = ZZ ? +R Application of the generalized bilinear transformation [4, 5] in this context leads to the approximate substitutions (70) Z = st ? sxr ?! 2Tr tanh 21 (stT ? sxh) Z = st + sxr ?! 2Tr tanh 21 (stT + sxh) (71) 1
2
3
3
0
3
0
4
0
4
0
2
2
4
2
2
where the relationship
2 = 2r =: r (72) T h between spatial and temporal step sizes h and T , respectively, holds. By insertion of the right hand side terms of (70,71) into (68,69) and the appropriate choice R = r the scattering quantities S and S are transformed in the following way 2
0
0
0
1
2
S ?! ?e?stT esx h =: ?zt? zx ; (73) S ?! ?e?stT e?sx h =: ?zt? zx? : (74) As a consequence the lattice twoport of g. 4a possesses the WDF equivalent of g. 4b with S and S being replaced by the right hand side terms of (73,74). Here zt and zx denote temporal delay and spatial shift operators of step sizes T and h, respectively. Thus after discretization via the generalized bilinear transformation, the following relationships in the space-time domain between the corresponding discrete-time/discrete-space wave quantities b0(n; k) = a0(n + 1; k ? 1) (75) 00 00 b (n; k) = a (n ? 1; k ? 1) (76) in g. 4b result. Wave quantities a0; a00; b0, and b00 here are the space-time equivalents to A0; A00; B 0, and B 0 of g. 4b. Applying WDF-concepts, namely system description by wave quantities and bilinear transformation, in a similar way to the rest of the circuit of g. 2c leads to the twodimensional WDF depicted in g. 5. 3-port series adaptors have been used for the representation of port interconnections. The port resistances were chosen as follows: R = r ; R = gr (77) R = 2(l T? ) ; R = 2(crT? ) (78) (79) R = R = r = 2hr = 2T : 1
1
1
2
1
2
1
3
2
5
4
6
2
0
14
2 2
2 2
1
Figure 5: Multidimensional WDF corresponding to the Kirchho circuit of g. 2c The minus sign of the factors ? in the gure in contrast to g. 4b is due to equations (73, 74). Choosing = l; r = l=c (cf. (65)) we get port resistances 1 2
q
2
R =R =0 2
(80)
4
which means that the two 3-port adaptors in g. 5 degenerate to 2-port adaptors. Thus the mere time delay elements in the MWDF structure have no e ect anymore. This special case corresponds to an integration along the characteristics of (1) for a transmission line without distortions by using the rectangular rule of integration. Note that the generalized bilinear transformation corresponds to a generalized rectangular rule of integration in the space-time domain [4, 5]. Of essential importance is how the desired voltage and current at every grid point (x; t) = (nh; kT ) may be obtained from the wave quantities of g. 5. Recalling (67) this can be accomplished in various ways by calculating e.g. b = a ?b = ? b i = i = a 2? (81) R 2R 2R ? b = r a ? b = ?r b u = r i = r a 2R (82) 2R 2R at every space location of interest. Concerning the operation of the structure in g. 5 let us give some additional remarks. In order to calculate voltage and current along the transmission line the MWDF-structure must be applied for every interesting spatial sampling point x = nh; h = NA separately. So for N + 1 locations x 2 [0; A] the structure must be used N + 1 times. A closer look at the structure furthermore shows that in order to calculate wave quantities at (x; t) = (nh; kT ) the appropriate quantities at grid points ((n ? 1)h; (k ? 1)T ), (nh; (k ? 1)T ) 1
2 2
5
5
2
2
0
2
1
2
6
6
0
15
2
1
4
4
4
2
3
3
and ((n + 1)h; (k ? 1)T ) must be known. This means that wave quantities of the former timestep located at the same place as well as at adjacent points contribute to grid point (x; t) = (nh; kT ). This dependency is depicted in g. 6.
Figure 6: Dependence of wave quantities between subsequent time steps Furthermore g. 6 reveals that care has to be taken in the calculation of wave quantities at boundary points, i.e. at x = 0 and x = A, because the structure requires knowledge of quantities outside the actual transmission line. How this problem may be circumvented in an appropriate manner will be shown in section 4.3. A crucial point in simulating partial di erential equations is the ratio of step sizes such that stability is guaranteed. Taking into account inequalities (64) and the coupling between spatial and temporal step sizes by (72) one easily veri es that the upper bound for the temporal step size p T h lc = vh (83) 0
must not be exceeded. This is the well-known Courant-Friedrichs-Levy (CFL) bound which is the stability criterion for many approaches to the numerical solution of the wave equation, i.e. for the situation of a wave travelling freely across a lossless transmission p line (r = g = 0) with phase velocity v = 1= lc. As has been pointed out already the same MWDF structure could have been obtained by a direct space-time approach based on the two partial di erential equations (1). For more information towards this approach we refer to [4, 5]. 0
4.3 Boundary conditions
In the context of the presented WDF treatment the boundary conditions (4) have to be formulated as
u(0; t) = yb (t) u(A; t) = yb (t) 1
2
x = 0; x = A:
(84) (85)
Let us start with the realization of the left boundary condition (x = 0). The crucial point in combining voltage and current calculation at x = 0 with the realization of (84) is the fact that the structure in g. 5 requires knowledge of the state x (?h; kT ) in order to calculate wave quantities a and a for (x = 0; t = (k + 1)T ). Our aim is to nd an expression for x (0; (k +1)T ) incorporating boundary condition (84) without being obliged 2
5
6
2
16
to know x (?h; kT ). To this end we make use of equations (55, 82) and the 3-port series adaptor equation b = a ? (a + a ) (86) with the adaptor coecient of the right series adaptor. This coecient has to be calculated from the appropriate port resistances as it is usual in WDF theory. Equation (86) may also be written as a ? b = (a + a ) : (87) Thus with identity (55) the boundary condition (84) turns into r (a + a ) = y : (88) b 2r Solving for wave quantity a gives a = r2r yb ? a : (89) Moreover a closer look at the MWDF structure in g. 5 reveals the coupling 1 ?1 x a 2 2 (90) = 1 1 a ?2 ?2 x and in inverse form a 1 ?1 x : (91) = ?1 ?1 a x From summation of the two equations (91) 2
6
6
6
4
6
6
6
6
2 6
6
4
4
6
6
1
0
6
0
6
3
2
5
6 4
7 5
6
3
2
1
6 4
1
2 6
2
3
6 6 6 6 4
2 7 76 74 7 5
3
1
76 54
6 4
7 5
2
3
32
2
7 5
4
5
7 5
6
2
x = ?(x + 2a ) 2
1
6
results. Substituting a here by (89) then gives x = ?x ? r4r yb + 2a (92) = ? x + r4r yb + 2x ; where we used the identity a = ?x . Equation (92) constitutes the result we were looking for. It enables us to perform the calculation of voltage and current at x = 0 plus satisfying boundary condition (84). So in every timestep we calculate x ; x and x for x = 0 with the required operations of g. 5 (remember, this poses no problem at all) and then determine x for x = 0 with (92). It is quite obvious that we arrive at the similar expression (93) x = ? x + r4r yb + 2x 6
2
0
1
1
2 6
4
!
0
1
4
2 6
1
4
4
1
3
4
2
!
1
2
0
2 6
17
2
4
in order to satisfy computability and validity of boundary condition (85) at x = A. The outlined boundary treatment can easily be modi ed to accomplish a given current at the boundary. To this end on simply has to use (81) instead of (82) and the 3-port series adaptor relation of the left series adaptor relating b and a and subsequently perform similar steps as described right before. Because of the strict locality of the MWDF algorithm mixed boundary conditions with voltage given at the one side and and current at the other side can be realized. Other types of boundary conditions, such as voltage sources with internal resistances or load resistances can be incorporated into the MWDF algorithm [11]. Terminations with more general impedances are under investigation. Finally we remark that it is equally possible to solve initial-value problems with the MWDF structure of g. 5. Prescribing initial values for voltage and current across the line, i.e. u(x; t = 0) and i(x; t = 0) ; 0 x A ; (94) these values may be transformed into initial values of the individual state variables x to x. 5
5
1
4
18
5 Simulation results As a rst example we consider a transmission line without distortions. The parameters of the transmission line are: A = 2 m, r = 10 /m, l = 10 mH/m, g = 4 mS/m and c = 4 F/m. We prescribe the boundary value 0 t 16 T otherwise
(
V ; yb (t) = 01 V ; 1
(95)
at the beginning of the line and yb (t) = 0; 8 t. The simulation methods described in sections 3 and 4 have been used for the simulation of the voltage distribution across the line. With N = 32 we get (96) h = 0:0625 m and T = vh = 12:5 s as the spatial and temporal step sizes. The result is depicted in g. 7 for times t 2 [ 0; 1:5 ms ]. 2
0
y(x,t) in V
1
0.5
0
−0.5 2
0 1.5 0.5
1 1 t in ms
0.5 1.5
0
Figure 7: Simulation result according to section 3 for the rst example 19
x in m
Figure 8: Boundary value yb (t) and simulation result at x = A=2 according to section 3 for the rst example 1
Furthermore g. 8 shows the the excitation at the line beginning and the simulation result at location x = A=2. In this case where the transmission line causes no distortions, also an analytical solution is possible. It has been used as a reference, against which the simulation results have been compared. For the method of section 3 the simulation error turned out to be less than 10? V, since, as has already been mentioned, equ. (45) provides the exact solution to (43). Fig. 8 clearly demonstrates the absence of numerical dispersion in the calculation using the FTA approach for the chosen spatial and temporal step sizes. On the other hand the maximumdeviation for the MWDF approach is depicted in g. 9 for step sizes T 0 = 2Tm ; h0 = 2hm ; m = 0(1)6 (97) versus the required computation time. All computations were run under MATLAB on a 25 MHz PC-80386. It can be seen that the simultaneous reduction of h and T by a factor of 2 leads to a decrease in the simulation error by a factor of 4. Also shown in g. 9 are the results of a computation method originating from a discretization of di erential operators in (2) in such a way that an explicit di erence equation results. This discretization of (2) was performed using central di erences with respect to grid point (x; t) = (nh; kT ). As a second example we consider the coaxial cable presented in [13] at constant temperature. Its characteristic parameters are: A = 10 km, r = 46:80 m /m, l = 5:60 H/m, g = 0 and c = 120 pF/m. This time we take the raised cosine impulse 12
yb (t) = 1
8 < :
sin
2
t 24 T
V ; 0V ;
0 t 23 T otherwise
(98)
as the excitation of the line. Again the simulation methods of sections 3 and 4 were applied with N = 32 and (99) h = 312:5 m and T = vh 8:10 s : 0
20
Figure 9: Maximum simulation error versus computation time for the rst example. EXPL: explicit di erence equation, MWDF: multidimensional WDF approach The voltage distribution obtained by simulation is shown in g. 10. Additionally g. 11 shows the boundary value yb (t) as well as the simulation result at the selected point x = A=2. Also for this case of a transmission line with distortions a reference solution has been calculated using a technique for numerical inverse Laplace{transformation. It was found to be accurate up to 10? V. Again the maximum simulation error versus required computation time is to be seen in g. 12 for the methods of sections 3 and 4 and the aforementioned explicit numerical method. Starting with step sizes h and T as in (99), the subsequent grid re nement according to (97) for both methods results in a reduction of the simulation error by the factor 4. In terms of the required computation time the MWDF method appears to be superior whereas the method of section 3 produces results that are four times more accurate for the same spatial and temporal sampling grid than the other methods. Moreover, restricting the in nite series in (51) and (53) to the baseband values of the nite Fourier transforms from ?N < N produces still more accurate results for the special boundary value yb (t) of (98). 1
11
1
21
y(x,t) in V
1 0.8 0.6 0.4 0.2 0 -0.2 0
10 8
0.2 6
0.4 4
0.6 2
0.8 1
t in ms
0
x in km
Figure 10: Simulation result according to section 3 for the second example
Figure 11: Boundary value yb (t) and simulation result at x = A=2 according to section 3 for the second example 1
22
Figure 12: Maximum simulation error versus computation time for the second example. FTA: functional transformation approach, EXPL: explicit di erence equation, MWDF: multidimensional WDF approach
23
6 Summary We have presented two methods for the discrete simulation of continuous systems applied to the problem of wave propagation across a uniform transmission line. The method of section 3 applies functional transformations to a single higher order partial di erential equation and subsequently develops a discrete system to match the continuous system's behaviour by discretizing the input/output characteristic which consists of a 2D transfer function. Thus, in the presented form, its use is restricted to partial di erential equations with constant coecients, but as the simulation results show, it produces better results than the MWDF approach of section 4 for a given spatial and temporal sampling grid. Applying functional transformations in a more general way this approach also allows to treat nonlinear partial di erential equations. On the other hand the MWDF method starts from the original set of partial di erential equations and makes use of generalizations to 1D WDF theory in order to obtain a simulating discrete system. As has been shown in [5] the case of non{constant coecients as well as more enhanced problems such as the 3D Maxwell equations [6] may be treated with this approach, too. In case of the problem under investigation the method produces acceptable results at relatively low computational expense. Summarizing, it was shown that both methods are well{suited for the numerical solution of the given problem, each of them exhibiting speci c advantages.
Acknowledgement: The authors wish to thank P. Ste en for many useful discussions and T. Eckart and M. Schetelig for programming work.
24
References [1] R. Bracewell, The Fourier{Transformation and its Applications (McGraw{Hill, New York, 1978). [2] G. Doetsch, Handbuch der Laplace{Transformation, Band III (Birkhauser Verlag, Basel, 1956). [3] A. Fettweis, Wave digital lters: theory and practice, Proc. IEEE 74 (2) (1986) 270-327. [4] A. Fettweis and G. Nitsche, Transformation approach to numerically integrating PDEs by means of WDF principles, Multidimensional Systems and Signal Processing (2) (1991) 127-159. [5] A. Fettweis and G. Nitsche, Numerical integration of partial di erential equations using principles of multidimensional wave digital lters, Journal of VLSI Signal Processing (3) (1991) 7-24. [6] A. Fettweis, Multidimensional wave digital lters for discrete-time modelling of Maxwell's equations, Int. Journal for Numerical Modelling: Electronic Networks, Devices and Fields 5 (1992) 183-201. [7] A. Fettweis, Multidimensional wave-digital principles: from ltering to numerical integration, in: Proceedings Int. Conf. Acoustics, Speech, and Signal Processing , Adelaide, Australia (1994) VI-173{VI-181. [8] E.I. Jury, Theory and Application of the z{Transform{Method (John Wiley, New York, 1964). [9] Z. Kowalczuk, Discrete approximation of continuous-time systems: a survey, IEE Proceedings-G , 140 (4) (1993) 264-278. [10] H. Krau and R. Rabenstein, Numerical solution of the telegraph equation by multidimensional digital lters, in: Proc. Int. Conf. DSP & CAES , Nicosia, Cyprus (1993) 170-175. [11] H. Krau and R. Rabenstein: Application of multidimensional wave digital lters to boundary value problems, to be published in IEEE Signal Processing Letters 2 (7) (1995). [12] C.-C.J Kuo and B.C. Levy, Discretisation and solution of elliptic PDEs { A digital signal processing approach, Proc. IEEE 78 (12) (1990) 1808-1842. [13] H. Mokhtari, A. Nyeck, C. Tosser{Roussey, A. Tosser{Roussey: Finite di erence method and Pspice simulation applied to the coaxial cable in a linear temperature gradient, IEE Proceedings{A 139 (1) (1992) 39-41. [14] A.V. Oppenheim and A.S. Willsky: Signals and Systems (Prentice{Hall Inc., Englewood Cli s, USA, 1983). 25
[15] R. Rabenstein and P. Ste en, Stability of image sequence processing systems, Archiv fur Elektronik und U bertragungstechnik (AEU ) 37 (7/8) (1983) 261-266. [16] R. Rabenstein, Simulation of distributed linear systems, in: Proc. of the EUROSIM'92 Simulation Congress, Capri, Italy (1992) 41-46. [17] R. Rabenstein, Simulation of linear continuous systems with distributed parameters, Simulation Practice and Theory 1 (1993) 93-107. [18] R. Rabenstein, Discrete simulation of linear multidimensional continuous systems, Multidimensional Systems and Signal Processing 5 (1994) 7-40. [19] R. Rabenstein and H. Krau, Discrete simulation of uniform transmission lines by multidimensional digital lters, to be published in: Int. Journal for Numerical Modelling: Electronic Networks, Devices and Fields . [20] R. Rabenstein, Discrete simulation of dynamical boundary value problems, to be published in: Proc. of the EUROSIM'95 Simulation Congress, Vienna, Austria (1995) [21] R. Sauer and I. Szabo, Mathematische Hilfsmittel des Ingenieurs, Teil I (Springer{ Verlag, Berlin, 1967). [22] H.W. Schuler, A signalprocessing approach to simulation, Frequenz 35 (1981) 174184.
26
Abbreviated title: Simulation of Wave Propagation Abstract
The application of multidimensional digital ltering theory to the simulation of continuous systems with distributed parameters has been treated in various recent contributions. Here, two di erent concepts will be compared, using a lossy transmission line of nite extent as an example. One approach is based on multidimensional wave digital lters, the other on the derivation of a simulating structure by functional transformations for each independent variable. The comparison includes a review of the underlying theoretical concepts as well as the presentation of the design process for the digital system and simulation results for various examples.
Keywords:
Computer Simulation, Wave Propagation, Digital Signal Processing, Digital Filters, Wave Digital Filters
Escola Politecnica da USP, S~ao Paulo, Brazil, worked with a post-doc fellowship at the University of Erlangen{Nurnberg, Germany.
1
1 Introduction The simulation of lumped parameter systems by onedimensional digital signal processing methods is well known. It has been used also for the simulation of distributed parameter systems like transmission lines by applying a spatial discretization to the describing partial di erential equation. This procedure is known as the method of lines and gives a sequence of interconnected lumped parameter systems. Each of these can be replaced by an appropriate onedimensional digital lter for simulation purposes. Multidimensional digital signal processing methods have been con ned to image and seismic processing for a long time. Recently, the application of the theory of multidimensional digital lters to the numerical solution of partial di erential equations and to the simulation of continuous systems with distributed parameters has been treated in the literature [4, 5, 12, 17, 18]. These methods allow the replacement of the distributed parameter system directly by a multidimensional digital lter, without the need for an initial conversion of the continuous system into a series of lumped parameter systems. Two of these methods will be compared here with respect to the simulation of a lossy transmission line. The transmission line problem has been chosen as an example of a hyperbolic partial di erential equation since both methods can be applied in a more general way to a variety of wave propagation phenomena. What is more, the method of section 3 which is based on a functional transformation approach is applicable to partial di erential equations of the parabolic type. Thus it it is quali ed e.g. for the simulation of heat and mass transfer problems [17, 18, 20]. The method based on multidimensional wave digital lter principles can be extended to propagation problems in more than one spatial dimension [4, 5, 6] and to more elaborate boundary value problems [11]. Section 2 gives a mathematical description of the transmission line problem. The digital ltering concepts to be compared are presented in sections 3 and 4. The examples for the numerical experiments and the results of the comparison are given in section 5.
2
2 Problem Statement We consider a transmission line of length A described by the two partial di erential equations @ u(x; t) + l @ i(x; t) + ri(x; t) = 0 @x @t (1) @ i(x; t) + c @ u(x; t) + gu(x; t) = 0 @x @t for voltage u(x; t) and current i(x; t) along the line. Both voltage and current depend on the space coordinate x with 0 x A and on the time coordinate t with 0 t. The coecients are given in terms of the following physical quantities: l : serial inductance in H=m, c : shunt capacitance in F=m, r : serial resistance in =m, g : shunt conductance in S=m. The two partial di erential equations (1) of rst order may be transformed into the second-order partial di erential equation @ y(x; t) = lc @ y(x; t) + (lg + rc) @ y(x; t) + rgy(x; t) ; (2) @x @t @t the so-called telegraph equation, where y(x; t) denotes voltage u(x; t) or current i(x; t). Unless otherwise stated y(x; t) in the sequel stands for the voltage u(x; t). The line shall be initially discharged, that is y(x; 0) = 0 ; (3) @ y(x; t) = 0: @t t At the ends of the line, we prescribe boundary conditions of the rst kind y(0; t) = yb (t) ; (4) y(A; t) = yb (t) or of the second kind @ y(x; t) = yb (t) ; @x x (5) @ y(x; t) = yb (t) : @x x A The transmission line is said to be free of distortions if r =g (6) l c holds. The following sections describe two multidimensional digital systems for the simulation of this continuous system. 2
2
2
2
=0
1
2
1
=0
2
=
3
3 Functional Transformation Approach (FTA) This section describes the specialization of a general method for the discrete simulation of linear multidimensional continuous systems to the problem stated previously. A full account of the method and an application to the heat ow equation is given in [16, 17, 18]. The application to wave propagation problems is given in [19]. The key steps are the representation of the continuous system by a transfer function, the derivation of the transfer function of a corresponding discrete system and nally the speci cation of a realizing structure, which can be easily cast into a computer program. This procedure is well known from the analysis of onedimensional systems. For time dependent problems like e.g. electrical networks, the Laplace{transformation [14, 21] provides an elegant tool to represent an ordinary di erential equation along with its initial condition(s) by appropriate transfer functions. This method works well because the integration range of the Laplace{transformation matches the region on which the inital value problem is de ned (0 t < 1 in both cases). An extension of this method to the space{ and time{dependent initial{boundary value problem described in section 2 calls for a further functional transformation for the space variable with 0 x A. Due to reasons explained in [18], it is of advantage to use an odd (even) continuation of the spatial functions into the interval ?A < x < 0 for boundary conditions of the rst (second) kind and a periodic repetition with period 2A outside of ?A x A. Then the nite Fourier{transformation with an integration range of ?A x A will include boundary conditions with respect to x in the same way into a transfer function as the Laplace{transformation does for the initial conditions with respect to t. Mixed type boundary conditions, i.e. boundary conditions of the rst kind on the one end and of the second kind at the other end, may be accomplished using a periodic repetition with period 4A. See [18, 21] for a discussion of the nite Fourier{ transformation and its properties. A de nition will be given later in this section. The procedure will now be demonstrated for the telegraph equation (2).
3.1 Continuous System
In order to obtain a transfer function from the partial di erential equation representation (2) of the transmission line, we at rst apply the one{sided Laplace{transformation with respect to the time coordinate t 1
Z
Lfy(x; t)g = Y (x; s) = y(x; t)e?stdt ; s 2 C :
(7)
I
0
Since the in uence of the boundary conditions is a linear superposition of yb (t) and yb (t), we can assume yb (t) = 0 for simplicity and apply the Laplace{transformation only to the non{vanishing boundary value at x = 0. The result is @ Y (x; s) = (s) Y (x; s) ; Y (s) = Lfy (t)g ; (8) b b @x with the propagation function 1
2
2
2
2
1
2
q
(s) = (sl + r)(sc + g) : 4
1
(9)
For later convenience, we introduce the phase velocity v and the abbreviations and by (10) v = p1 ; = ? rl ; = ? gc : lc The squared propagation function is thus (11)
(s) = v1 (s ? )(s ? ) : The next step is the application of the nite Fourier{transformation ( n. FT) 0
1
2
1
0
2
2
1
2 0
FAfY (x; s)g = Y (; s) =
A
Z
2
Y (x; s)e?j =A x dx ; 2 ZZ (
?A
)
(12)
1 1 FA = Y (x; s) = 2A Y (; s)ej =A x (13) ?1 in space direction. For boundary conditions of the rst kind, we have to assume an odd continuation and periodic repetition in the sense stated above. With the derivation theorem of the nite Fourier{transformation [18, 21], we obtain ? Y (; s) = (s) Y (; s) + 2j Yb (s) (14) with (15) = A : This algebraic equation can be solved for Y (; s) to obtain Y (; s) = G b (; s) Yb (s) (16) with G b (; s) = ? (2sj : (17) ) + Note, that equations (16,17) are an alternate description of the continuous system given by (2,3,4). It is not based on the space and time coordinates x and t, but on the corresponding transformed variables and s, respectively. We can envisage G b (; s) in (17) as the twodimensional transfer function between the input Yb (s) and the output Y (; s). For a given transfer function, (16) is a concise description of the problem previously given by the partial di erential equation (2) together with initial (3) and boundary conditions (4). Using (10,11,15), we can write the denominator of the transfer function (17) as (18)
(s) + = v1 (s ? )(s ? ) + ! where ! = A v = 2T : (19) The quantity p (20) T = 2vA = 2A lc ?1 fY (; s)g
X
(
)
=
2
2
1
1
1
1
2
2
1
2
1
1
2
h
2
2 0
0
0
0
5
0
2
2
i
can be regarded as the time it takes for a wave with phase velocity v to travel across an (extended) line of length 2A. Equation (18) is a polynomial in s with the zeros 0
s
s1 = ( ) = +2 ?2 ? ! ; (21) which constitute the poles of the transfer function (17) v G b (; s) = ?j (s ? s (2! : (22) 1 ))(s ? s1 ( )) We can describe the in uence of a boundary value yb (t) in the time domain by 1
1 2
2
1
2
2
2
2
0
1
1
2
1
y(; t) = gb (; t) t yb (t) (23) with y(; t) = FAfy(x; t)g = L? fY (; s)g and gb (; t) = L? fG b (; s)g = ?j s (2)!? vs ( ) [exp(s1 ( )t) ? exp(s1 ( )t)] ? (t) : 1 1 (24) (? (t) is the unit step function.) By setting yb (t) = (t), we see that gb (; t) describes the reaction of the transmission line to an impulse at the boundary. Equation (24) speci es this reaction as an in nite set of time functions, each valid for a di erent discrete spatial frequency with ?1 < < 1 (see (12,13)). The corresponding description in the space{time domain can be obtained by inverse nite Fourier{transformation (13) as 1 gb (x; t) = 21A gb (; t)ej x : (25) ?1 1
1
1
1
1
0
1
1
1
2
2
1
1
1
0
1
X
1
1
=
3.2 Discrete System
The equations (2,3,4,25), (16,17), and (23,24) give di erent representations of the transmission line in the (x; t), (; s) and (; t) domain, respectively. We will now specify a procedure to obtain a multidimensional digital lter for the simulation of the continuous transfer function G b (; s). Its purpose is to calculate a close approximation to the output y(x; t) of the continuous system at the grid points (x; t) = (nh; kT ), where h and T are the spatial and temporal step sizes respectively, and n and k are the corresponding discrete spatial and temporal variables. The design of the twodimensional digital lter will proceed in two steps. At rst we will discretize the time variable, based on equations (23,24) separately for each spatial frequency . The second step is the discretization of the space variable, based on the inverse nite Fourier{transformation (13). 1
3.2.1 Discretization of the Time Variable
The discretization of the time variable is based on the well known impulse{invariant transformation [9, 14, 22]. To this end, we consider gb (; t) as the impulse response of an 1
6
analog lter with a discrete parameter . Then the impulse response gbh (; k) of a hybrid system that is discrete in time, but continuous in space is 1
gbh (; k) = T gb (; kT ) :
(26)
1
1
From (24), we obtain
gbh (; k) = ?j s (2)!? vs ( ) z1k ( ) ? z1k ( ) ? (k) 1 1 h
0
1
1
i
1
2
2
(27)
1
with
z1i( ) = exp(s1i( )T ) ; i = 1; 2 : (28) Z{transformation [8] gives the transfer function with z G hb (; z) = z + cb(b()z)+ (29) c ( ) ; where bb ( ) = ?2j!v z1 ( ) ? z1 ( ) ; (30) s1 ( ) ? s1 ( ) c ( ) = ?(z1 ( ) + z1 ( )) ; (31) c ( ) = z1 ( ) z1 ( ) : (32) For further evaluation of the coecients bb ( ), c ( ), c ( ), we have to distinguish, whether the zeros of (18), as given by (21) are both real and distinct, equal to each other or a complex conjugate pair. In terms of (18), the cases j j >< 2?! (33) are possible. The results for all three cases can be summarized as bb ( ) = ?2j!v T exp( T ) sin ! T ; (34) ! T c ( ) = ?2 exp( T ) cos ! T ; (35) c ( ) = c = exp(2 T ) : (36) with = 21 ( + ) = ? 21 rl + gc (37) ! = ! ? 2?! : (38) 1
1
2
1
1
0
0
1
1
2
1
2
1
0
2
1
2
1
=
1
0
1
0
0
0
2
1
1
0
2
0
0
0
1
2
s
2
1
2
2
According to the di erent cases given by (33), ! may be real, zero, or imaginary. The special case of a transmission line without distortions (6) is characterized by = = and ! = ! . Using the transfer function G hb (; z) from (29), we can describe the hybrid system by Y h (; z) = G hb (; z)Ybd (z) (39) 1
1
1
7
1
2
0
with the input Ybd (z) = Z fyb (kT )g. Insertion of (29), reordering of terms and inverse z{transformation gives yh(; k) = ?c ( )yh (; k ? 1) ? c yh(; k ? 2) + bb ( )ybd (k ? 1) : (40) 1
1
1
0
1
1
By inverse n. FT, we could obtain
yh(x; k) = ?c (x) x yh(x; k ? 1) ? c yh(x; k ? 2) + bb (x)ybd (k ? 1) ; 1
0
1
1
(41)
i. e. the output of a continuous space, discrete time system with
yh(x; k) y(x; kT ) :
(42)
Since we are not interested in hybrid simulation, we use the hybrid system (39,40,41) as an intermediate stage to a fully discrete system. The remaining step | the discretization of the space variable | will be covered next.
3.2.2 Discretization of the Space Variable
A rst attempt to a fully discrete system is the evaluation of (41) at discrete points in space x = nh. This poses no problem for the second and third term on the right hand side. However, the convolution term cannot be solved exactly on a discrete grid for general c (x) and yh(x; k ? 1). We will therefore focus our attention on a discrete approximation of the term 1
f h (x; k) x
= c (x) x yh(x; k) nh h
1
=
i
n
= FA? c ( )yh(; k) nh 1
x=
1
o
x=nh
:
(43)
An exact solution would always be possible, if yh(x; k) was a bandlimited function in space. This prerequisite is not met, since bb (x) is not bandlimited (see its n. FT in (34)). We may, however, apply the rectangular rule of integration to intervals of length (44) h = NA ; N 2 IN for the convolution (43). This yields the approximate solution 1
1
f d (n; k)
=
f h(x = nh; k) =
h
N X =?N +1
A
Z
?A
c ()yh (nh ? ; k)d 1
c (h)yh((n ? )h; k) = 1
N
X
=?N +1
cd()yd(n ? ; k) ; 1
(45)
where cd(n) and yd(n; k) denote the samples of h c (x) and yh(x; k) at locations x = nh. Choosing the xed relationship h=v T (46) 1
1
0
In case of a transmission line without distortions (see (6)) equation (45) even provides the exact solution. 1
8
between the spatial step size h and the temporal sampling interval T it is shown in [19] that cd(n) and yd(n; k) in (45) may be calculated by means of an inverse Discrete Fourier Transform (DFT) of length 2N . To this aim the Discrete Fourier Transform is applied in the form 1
n
N
o
y^d(; k) = DFT N yd(n; k) = 2
X
n=?N +1 N 1 X
yd(n; k) exp(?j N n ) ;
y^d(; k) exp(j N n ) : yd(n; k) = DFT?N y^d(; k) = 2N ?N Then one gets yd(n; k) = DFT?N y^d(; k) with the DFT coecients y^d(; k) = h1 yh(; k) ; ?N < N : On the other hand cd (n) may be obtained via 2
1
n
o
=
2
+1
n
1
o
(47) (48) (49) (50)
1
cd(n) = DFT?N 1
2
1
8 < X
1
:
i=?1
c ( ? i 2N ) 1
9 =
(51)
;
up to any desired accuracy. Concerning the convergence and the practical evaluation of the in nite series in (51) we refer to [19]. Thus we can write (45) in the form
f d(n; k) =
cd(n)
2N yd(n; k)
1
(52)
with 2N denoting cyclic convolution of length 2N . In summary, the continuous convolution term f h(x; k) from (43) can be evaluated approximately at a spatial grid point (x = nh = nv T ) via a discrete cyclic convolution between the voltage distribution yd(n; k) and a coecient sequence cd(n). The latter can eciently be computed with fast algorithms for the DFT using the spectral components according to (35,50,51). Note, that also the third term on the right hand side of (41) is obtained with the DFT up to any desired accuracy as 0
1
bdb (n) = DFT?N 1
by
2
1
8 < :
1 hi
1
X
?1
=
9 =
bb ( ? i 2N ) : 1
(53)
;
All together, the fully discrete system with the grid points (x; t) = (nh; kT ) is given
yd(n; k) = ? cd(n) 1
2N yd(n; k ? 1) ? c yd(n; k ? 2) + bdb (n) ybd (k ? 1) : 0
1
1
(54)
A realizing structure is shown in g. 1. It is a second order recursive system with respect to the discrete time variable k. In each time step a vector of 2N spatial samples is 9
Figure 1: Structure of the discrete system (see (54)) processed in parallel. The boundary value ybd (k) acts as the input of the system. Its impact on the spatial voltage distribution yd(n; k) is modelled by the sequence bdb (n). Finally we remark that the structure according to g. 1 is not based on a replacement of di erential operators in (2) by di erence operators, as is done in many other numerical methods. Consequently the CFL-condition (cf. section 4.2) governing the stability does notpapply here. This means that the structure is also stable for step sizes T > h=v ; v = 1= lc. A stability criterion for recursive systems of this type can be found in [15]. 1
1
0
10
0
4 Multidimensional Wave Digital Filter (MWDF) In this section a di erent approach based on multidimensional wave digital lter (MWDF) theory will be outlined in brevity. This method has been presented in various publications by Fettweis and Nitsche, e.g. applied to the problem of two conducting plates in [4] and [5]. For a full account of the method we therefore refer to [4, 5, 6], for a survey of MWDF applications to [7]. The key steps of this approach are the derivation of a multidimensional Kirchho circuit from the underlying system of partial di erential equations and its simulation by a multidimensional wave digital lter. Thus a passive continuous system will be simulated by a passive discrete system.
4.1 Derivation of the Kirchho Circuit
In contrast to section 3 now the two partial di erential equations in (1) for voltage u(x; t) and current i(x; t) across a transmission line, having constant parameters r; l; g and c, serve as the basis for the investigations. In order to deal with same quantities a normalization according to i (x; t) := i(x; t) and i (x; t) := u(rx; t) (55) is used where r is an arbitrary positive resistance. One then gets the two equations @ i (x; t) + l @ i (x; t) + r i (x; t) = 0 r @x (56) @t @ i (x; t) + r c @ i (x; t) + r g i (x; t) = 0 (57) @x @t 1
2
2
2
2
2
1
1
2
1
2
2
2
We now go on describing the frequency domain version of the method which is completely equivalent to a direct space-time approach in case of constant parameters r; l; g and c as was shown in [4, 5]. Assuming steady state behaviour one can set i (x; t) = I esT t ; = 1; 2 (58) with complex amplitudes I and vectors sT = [sx; st] ; t = [x; t]T : (59) Inserting (58) and (59) into (56,57) gives the algebraic equations (stl + r)I + sxr I = 0 (60) sxr I + (stcr + gr )I = 0 (61) where the last equation has been multiplied by the factor r . One can easily verify that equations (60) and (61) may be interpreted as mesh equations of a twodimensional Kirchho circuit as depicted in g. 2a. Furthermore g. 2b shows a modi cation of the circuit where an arbitrary positive constant has been introduced. Transformation of the T-con guration of g. 2b into a symmetric lattice structure gives the circuit of g. 2c where the impedances Z ; Z ; Z und Z are given by Z = st(l ? ) ; Z = st(cr ? ) (62) Z = st ? s x r ; Z = st + s x r : (63) 1
2 2
2 2
2 1
2 2
2
2
1
1
3
2
3
2
2
4
11
4
2 2
2
Figure 2: a) Equivalent Kirchho circuit; b) modi ed Kirchho circuit; c) equivalent symmetric lattice structure In order to guarantee passivity of the circuit the following inequalities l ; cr are supposed to hold. With the choice 2 2
(64)
s
= l ; r = cl the impedances Z and Z degenerate to Z = Z = 0: 2
1
(65)
2
1
2
4.2 Derivation of the MWDF
(66)
Similar to the 1D-case the transition from voltage and current quantities u and i to incident and re ected wave quantities a and b, respectively, plays an essential role in the derivation of the 2D WDF. Furthermore the bilinear transformation is applied, this time, however, in a generalized way as the frequency domain equivalent to a generalized trapezoidal rule of integration [4, 5]. 12
In accordance with [4, 5] incident and re ected waves in the frequency domain are introduced as voltage waves
A = U + R I ; B = U ? R I ;
(67)
where the proper orientation of voltage/current and wave quantities can be seen in g. 3. The factor R , called the port resistance, has to be chosen in a suitable way.
Figure 3: Voltage/current and wave quantity orientations of an electrical 1-port The full derivation of the MWDF for the simulation of the Kirchho circuit of g. 2c is beyond the scope of this paper. Here we just give an outline of major steps of this derivation and refer the reader to more profound WDF literature, especially [3] and [4, 5]. As pointed out in [3] the lattice twoport of g. 4a possesses the wave equivalent shown in g. 4b. Note that the twoport of g. 4a constitutes a part of the Kirchho circuit of g. 2c. In order to get the structure of g. 4b, the wave description (67) with orientation
Figure 4: a) Symmetric lattice twoport; b) equivalent wave representation as in g. 3 must be applied to each of the two ports of the lattice. By choosing the two port resistances to be equal to the same positive value R , which is not necessary but convenient for symmetry reasons, the scattering quantities S and S in g. 4b take the 0
1
13
2
form
R ; (68) S = ZZ ? +R R : (69) S = ZZ ? +R Application of the generalized bilinear transformation [4, 5] in this context leads to the approximate substitutions (70) Z = st ? sxr ?! 2Tr tanh 21 (stT ? sxh) Z = st + sxr ?! 2Tr tanh 21 (stT + sxh) (71) 1
2
3
3
0
3
0
4
0
4
0
2
2
4
2
2
where the relationship
2 = 2r =: r (72) T h between spatial and temporal step sizes h and T , respectively, holds. By insertion of the right hand side terms of (70,71) into (68,69) and the appropriate choice R = r the scattering quantities S and S are transformed in the following way 2
0
0
0
1
2
S ?! ?e?stT esx h =: ?zt? zx ; (73) S ?! ?e?stT e?sx h =: ?zt? zx? : (74) As a consequence the lattice twoport of g. 4a possesses the WDF equivalent of g. 4b with S and S being replaced by the right hand side terms of (73,74). Here zt and zx denote temporal delay and spatial shift operators of step sizes T and h, respectively. Thus after discretization via the generalized bilinear transformation, the following relationships in the space-time domain between the corresponding discrete-time/discrete-space wave quantities b0(n; k) = a0(n + 1; k ? 1) (75) 00 00 b (n; k) = a (n ? 1; k ? 1) (76) in g. 4b result. Wave quantities a0; a00; b0, and b00 here are the space-time equivalents to A0; A00; B 0, and B 0 of g. 4b. Applying WDF-concepts, namely system description by wave quantities and bilinear transformation, in a similar way to the rest of the circuit of g. 2c leads to the twodimensional WDF depicted in g. 5. 3-port series adaptors have been used for the representation of port interconnections. The port resistances were chosen as follows: R = r ; R = gr (77) R = 2(l T? ) ; R = 2(crT? ) (78) (79) R = R = r = 2hr = 2T : 1
1
1
2
1
2
1
3
2
5
4
6
2
0
14
2 2
2 2
1
Figure 5: Multidimensional WDF corresponding to the Kirchho circuit of g. 2c The minus sign of the factors ? in the gure in contrast to g. 4b is due to equations (73, 74). Choosing = l; r = l=c (cf. (65)) we get port resistances 1 2
q
2
R =R =0 2
(80)
4
which means that the two 3-port adaptors in g. 5 degenerate to 2-port adaptors. Thus the mere time delay elements in the MWDF structure have no e ect anymore. This special case corresponds to an integration along the characteristics of (1) for a transmission line without distortions by using the rectangular rule of integration. Note that the generalized bilinear transformation corresponds to a generalized rectangular rule of integration in the space-time domain [4, 5]. Of essential importance is how the desired voltage and current at every grid point (x; t) = (nh; kT ) may be obtained from the wave quantities of g. 5. Recalling (67) this can be accomplished in various ways by calculating e.g. b = a ?b = ? b i = i = a 2? (81) R 2R 2R ? b = r a ? b = ?r b u = r i = r a 2R (82) 2R 2R at every space location of interest. Concerning the operation of the structure in g. 5 let us give some additional remarks. In order to calculate voltage and current along the transmission line the MWDF-structure must be applied for every interesting spatial sampling point x = nh; h = NA separately. So for N + 1 locations x 2 [0; A] the structure must be used N + 1 times. A closer look at the structure furthermore shows that in order to calculate wave quantities at (x; t) = (nh; kT ) the appropriate quantities at grid points ((n ? 1)h; (k ? 1)T ), (nh; (k ? 1)T ) 1
2 2
5
5
2
2
0
2
1
2
6
6
0
15
2
1
4
4
4
2
3
3
and ((n + 1)h; (k ? 1)T ) must be known. This means that wave quantities of the former timestep located at the same place as well as at adjacent points contribute to grid point (x; t) = (nh; kT ). This dependency is depicted in g. 6.
Figure 6: Dependence of wave quantities between subsequent time steps Furthermore g. 6 reveals that care has to be taken in the calculation of wave quantities at boundary points, i.e. at x = 0 and x = A, because the structure requires knowledge of quantities outside the actual transmission line. How this problem may be circumvented in an appropriate manner will be shown in section 4.3. A crucial point in simulating partial di erential equations is the ratio of step sizes such that stability is guaranteed. Taking into account inequalities (64) and the coupling between spatial and temporal step sizes by (72) one easily veri es that the upper bound for the temporal step size p T h lc = vh (83) 0
must not be exceeded. This is the well-known Courant-Friedrichs-Levy (CFL) bound which is the stability criterion for many approaches to the numerical solution of the wave equation, i.e. for the situation of a wave travelling freely across a lossless transmission p line (r = g = 0) with phase velocity v = 1= lc. As has been pointed out already the same MWDF structure could have been obtained by a direct space-time approach based on the two partial di erential equations (1). For more information towards this approach we refer to [4, 5]. 0
4.3 Boundary conditions
In the context of the presented WDF treatment the boundary conditions (4) have to be formulated as
u(0; t) = yb (t) u(A; t) = yb (t) 1
2
x = 0; x = A:
(84) (85)
Let us start with the realization of the left boundary condition (x = 0). The crucial point in combining voltage and current calculation at x = 0 with the realization of (84) is the fact that the structure in g. 5 requires knowledge of the state x (?h; kT ) in order to calculate wave quantities a and a for (x = 0; t = (k + 1)T ). Our aim is to nd an expression for x (0; (k +1)T ) incorporating boundary condition (84) without being obliged 2
5
6
2
16
to know x (?h; kT ). To this end we make use of equations (55, 82) and the 3-port series adaptor equation b = a ? (a + a ) (86) with the adaptor coecient of the right series adaptor. This coecient has to be calculated from the appropriate port resistances as it is usual in WDF theory. Equation (86) may also be written as a ? b = (a + a ) : (87) Thus with identity (55) the boundary condition (84) turns into r (a + a ) = y : (88) b 2r Solving for wave quantity a gives a = r2r yb ? a : (89) Moreover a closer look at the MWDF structure in g. 5 reveals the coupling 1 ?1 x a 2 2 (90) = 1 1 a ?2 ?2 x and in inverse form a 1 ?1 x : (91) = ?1 ?1 a x From summation of the two equations (91) 2
6
6
6
4
6
6
6
6
2 6
6
4
4
6
6
1
0
6
0
6
3
2
5
6 4
7 5
6
3
2
1
6 4
1
2 6
2
3
6 6 6 6 4
2 7 76 74 7 5
3
1
76 54
6 4
7 5
2
3
32
2
7 5
4
5
7 5
6
2
x = ?(x + 2a ) 2
1
6
results. Substituting a here by (89) then gives x = ?x ? r4r yb + 2a (92) = ? x + r4r yb + 2x ; where we used the identity a = ?x . Equation (92) constitutes the result we were looking for. It enables us to perform the calculation of voltage and current at x = 0 plus satisfying boundary condition (84). So in every timestep we calculate x ; x and x for x = 0 with the required operations of g. 5 (remember, this poses no problem at all) and then determine x for x = 0 with (92). It is quite obvious that we arrive at the similar expression (93) x = ? x + r4r yb + 2x 6
2
0
1
1
2 6
4
!
0
1
4
2 6
1
4
4
1
3
4
2
!
1
2
0
2 6
17
2
4
in order to satisfy computability and validity of boundary condition (85) at x = A. The outlined boundary treatment can easily be modi ed to accomplish a given current at the boundary. To this end on simply has to use (81) instead of (82) and the 3-port series adaptor relation of the left series adaptor relating b and a and subsequently perform similar steps as described right before. Because of the strict locality of the MWDF algorithm mixed boundary conditions with voltage given at the one side and and current at the other side can be realized. Other types of boundary conditions, such as voltage sources with internal resistances or load resistances can be incorporated into the MWDF algorithm [11]. Terminations with more general impedances are under investigation. Finally we remark that it is equally possible to solve initial-value problems with the MWDF structure of g. 5. Prescribing initial values for voltage and current across the line, i.e. u(x; t = 0) and i(x; t = 0) ; 0 x A ; (94) these values may be transformed into initial values of the individual state variables x to x. 5
5
1
4
18
5 Simulation results As a rst example we consider a transmission line without distortions. The parameters of the transmission line are: A = 2 m, r = 10 /m, l = 10 mH/m, g = 4 mS/m and c = 4 F/m. We prescribe the boundary value 0 t 16 T otherwise
(
V ; yb (t) = 01 V ; 1
(95)
at the beginning of the line and yb (t) = 0; 8 t. The simulation methods described in sections 3 and 4 have been used for the simulation of the voltage distribution across the line. With N = 32 we get (96) h = 0:0625 m and T = vh = 12:5 s as the spatial and temporal step sizes. The result is depicted in g. 7 for times t 2 [ 0; 1:5 ms ]. 2
0
y(x,t) in V
1
0.5
0
−0.5 2
0 1.5 0.5
1 1 t in ms
0.5 1.5
0
Figure 7: Simulation result according to section 3 for the rst example 19
x in m
Figure 8: Boundary value yb (t) and simulation result at x = A=2 according to section 3 for the rst example 1
Furthermore g. 8 shows the the excitation at the line beginning and the simulation result at location x = A=2. In this case where the transmission line causes no distortions, also an analytical solution is possible. It has been used as a reference, against which the simulation results have been compared. For the method of section 3 the simulation error turned out to be less than 10? V, since, as has already been mentioned, equ. (45) provides the exact solution to (43). Fig. 8 clearly demonstrates the absence of numerical dispersion in the calculation using the FTA approach for the chosen spatial and temporal step sizes. On the other hand the maximumdeviation for the MWDF approach is depicted in g. 9 for step sizes T 0 = 2Tm ; h0 = 2hm ; m = 0(1)6 (97) versus the required computation time. All computations were run under MATLAB on a 25 MHz PC-80386. It can be seen that the simultaneous reduction of h and T by a factor of 2 leads to a decrease in the simulation error by a factor of 4. Also shown in g. 9 are the results of a computation method originating from a discretization of di erential operators in (2) in such a way that an explicit di erence equation results. This discretization of (2) was performed using central di erences with respect to grid point (x; t) = (nh; kT ). As a second example we consider the coaxial cable presented in [13] at constant temperature. Its characteristic parameters are: A = 10 km, r = 46:80 m /m, l = 5:60 H/m, g = 0 and c = 120 pF/m. This time we take the raised cosine impulse 12
yb (t) = 1
8 < :
sin
2
t 24 T
V ; 0V ;
0 t 23 T otherwise
(98)
as the excitation of the line. Again the simulation methods of sections 3 and 4 were applied with N = 32 and (99) h = 312:5 m and T = vh 8:10 s : 0
20
Figure 9: Maximum simulation error versus computation time for the rst example. EXPL: explicit di erence equation, MWDF: multidimensional WDF approach The voltage distribution obtained by simulation is shown in g. 10. Additionally g. 11 shows the boundary value yb (t) as well as the simulation result at the selected point x = A=2. Also for this case of a transmission line with distortions a reference solution has been calculated using a technique for numerical inverse Laplace{transformation. It was found to be accurate up to 10? V. Again the maximum simulation error versus required computation time is to be seen in g. 12 for the methods of sections 3 and 4 and the aforementioned explicit numerical method. Starting with step sizes h and T as in (99), the subsequent grid re nement according to (97) for both methods results in a reduction of the simulation error by the factor 4. In terms of the required computation time the MWDF method appears to be superior whereas the method of section 3 produces results that are four times more accurate for the same spatial and temporal sampling grid than the other methods. Moreover, restricting the in nite series in (51) and (53) to the baseband values of the nite Fourier transforms from ?N < N produces still more accurate results for the special boundary value yb (t) of (98). 1
11
1
21
y(x,t) in V
1 0.8 0.6 0.4 0.2 0 -0.2 0
10 8
0.2 6
0.4 4
0.6 2
0.8 1
t in ms
0
x in km
Figure 10: Simulation result according to section 3 for the second example
Figure 11: Boundary value yb (t) and simulation result at x = A=2 according to section 3 for the second example 1
22
Figure 12: Maximum simulation error versus computation time for the second example. FTA: functional transformation approach, EXPL: explicit di erence equation, MWDF: multidimensional WDF approach
23
6 Summary We have presented two methods for the discrete simulation of continuous systems applied to the problem of wave propagation across a uniform transmission line. The method of section 3 applies functional transformations to a single higher order partial di erential equation and subsequently develops a discrete system to match the continuous system's behaviour by discretizing the input/output characteristic which consists of a 2D transfer function. Thus, in the presented form, its use is restricted to partial di erential equations with constant coecients, but as the simulation results show, it produces better results than the MWDF approach of section 4 for a given spatial and temporal sampling grid. Applying functional transformations in a more general way this approach also allows to treat nonlinear partial di erential equations. On the other hand the MWDF method starts from the original set of partial di erential equations and makes use of generalizations to 1D WDF theory in order to obtain a simulating discrete system. As has been shown in [5] the case of non{constant coecients as well as more enhanced problems such as the 3D Maxwell equations [6] may be treated with this approach, too. In case of the problem under investigation the method produces acceptable results at relatively low computational expense. Summarizing, it was shown that both methods are well{suited for the numerical solution of the given problem, each of them exhibiting speci c advantages.
Acknowledgement: The authors wish to thank P. Ste en for many useful discussions and T. Eckart and M. Schetelig for programming work.
24
References [1] R. Bracewell, The Fourier{Transformation and its Applications (McGraw{Hill, New York, 1978). [2] G. Doetsch, Handbuch der Laplace{Transformation, Band III (Birkhauser Verlag, Basel, 1956). [3] A. Fettweis, Wave digital lters: theory and practice, Proc. IEEE 74 (2) (1986) 270-327. [4] A. Fettweis and G. Nitsche, Transformation approach to numerically integrating PDEs by means of WDF principles, Multidimensional Systems and Signal Processing (2) (1991) 127-159. [5] A. Fettweis and G. Nitsche, Numerical integration of partial di erential equations using principles of multidimensional wave digital lters, Journal of VLSI Signal Processing (3) (1991) 7-24. [6] A. Fettweis, Multidimensional wave digital lters for discrete-time modelling of Maxwell's equations, Int. Journal for Numerical Modelling: Electronic Networks, Devices and Fields 5 (1992) 183-201. [7] A. Fettweis, Multidimensional wave-digital principles: from ltering to numerical integration, in: Proceedings Int. Conf. Acoustics, Speech, and Signal Processing , Adelaide, Australia (1994) VI-173{VI-181. [8] E.I. Jury, Theory and Application of the z{Transform{Method (John Wiley, New York, 1964). [9] Z. Kowalczuk, Discrete approximation of continuous-time systems: a survey, IEE Proceedings-G , 140 (4) (1993) 264-278. [10] H. Krau and R. Rabenstein, Numerical solution of the telegraph equation by multidimensional digital lters, in: Proc. Int. Conf. DSP & CAES , Nicosia, Cyprus (1993) 170-175. [11] H. Krau and R. Rabenstein: Application of multidimensional wave digital lters to boundary value problems, to be published in IEEE Signal Processing Letters 2 (7) (1995). [12] C.-C.J Kuo and B.C. Levy, Discretisation and solution of elliptic PDEs { A digital signal processing approach, Proc. IEEE 78 (12) (1990) 1808-1842. [13] H. Mokhtari, A. Nyeck, C. Tosser{Roussey, A. Tosser{Roussey: Finite di erence method and Pspice simulation applied to the coaxial cable in a linear temperature gradient, IEE Proceedings{A 139 (1) (1992) 39-41. [14] A.V. Oppenheim and A.S. Willsky: Signals and Systems (Prentice{Hall Inc., Englewood Cli s, USA, 1983). 25
[15] R. Rabenstein and P. Ste en, Stability of image sequence processing systems, Archiv fur Elektronik und U bertragungstechnik (AEU ) 37 (7/8) (1983) 261-266. [16] R. Rabenstein, Simulation of distributed linear systems, in: Proc. of the EUROSIM'92 Simulation Congress, Capri, Italy (1992) 41-46. [17] R. Rabenstein, Simulation of linear continuous systems with distributed parameters, Simulation Practice and Theory 1 (1993) 93-107. [18] R. Rabenstein, Discrete simulation of linear multidimensional continuous systems, Multidimensional Systems and Signal Processing 5 (1994) 7-40. [19] R. Rabenstein and H. Krau, Discrete simulation of uniform transmission lines by multidimensional digital lters, to be published in: Int. Journal for Numerical Modelling: Electronic Networks, Devices and Fields . [20] R. Rabenstein, Discrete simulation of dynamical boundary value problems, to be published in: Proc. of the EUROSIM'95 Simulation Congress, Vienna, Austria (1995) [21] R. Sauer and I. Szabo, Mathematische Hilfsmittel des Ingenieurs, Teil I (Springer{ Verlag, Berlin, 1967). [22] H.W. Schuler, A signalprocessing approach to simulation, Frequenz 35 (1981) 174184.
26
