oriented spatial pattern formation in a four layer ... - Semantic Scholar

Report 2 Downloads 202 Views
Papers International Journal of Bifurcation and Chaos, Vol. 14, No. 4 (2004) 1209–1221 c World Scientific Publishing Company

ORIENTED SPATIAL PATTERN FORMATION IN A FOUR LAYER CMOS CELLULAR NEURAL NETWORK BERTRAM E. SHI Department of Electrical and Electronic Engineering, Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong [email protected] Received April 4, 2003; Revised May 1, 2003 Cellular neural networks have found applications both as a platform for high speed image processing, as well as a model for complex spatiotemporal dynamical phenomena. This paper describes spatial pattern formation in a CMOS VLSI implementation of a CNN. With four layers of 32 × 64 cells, this CNN contains the most recurrently interconnected 2D layers of any VLSI CNN implementation to date. With over 8000 cells, it is also one of the largest in terms of the total number of cells. The spatial patterns generated by this network are bands of activity with a preferred scale and orientation, both of which can be adjusted via external bias voltages. The VLSI implementation provides a physical substrate for experimenting with complex spatiotemporal phenomena in real time, while also being amenable to theoretical analysis. Our experimental characterization of the observed patterns are in excellent agreement with our theoretical predictions. Keywords: Spatial pattern formation; oriented patterns; Turing patterns; VLSI; cellular neural network; CMOS; nonlinear circuits.

1. Introduction A cellular neural network or cellular nonlinear network (CNN) is a spatial arrangement of locally coupled cells, each of which is a dynamical system. Although originally developed for applications in image processing, recent work has begun to explore its potential as a platform to investigate complex spatiotemporal behaviors [Chua, 1995, 1998]. Numerous theoretical and numerical simulation based studies have established that suitably connected CNNs are capable of exhibiting spatial and spatiotemporal patterns, as well as spatiotemporal chaos. While the CNN is not unique in this respect, as many systems of equations have been shown to exhibit such phenomena, one of the oftcited motivations for using the CNN to model these

phenomena is its compatibility with VLSI technology. VLSI implementations of CNNs provide a well-controlled physical substrate to study complex nonlinear spatiotemporal dynamics in real time. Computer simulation not only is more time consuming, but also introduces the possibility for numerical errors due to discretization in time and amplitude. Despite the compelling advantages, comparatively few papers have reported measurement and experimental characterization of complex spatiotemporal phenomena in VLSI implementations of CNNs. To exhibit bandpass spatial pattern formation, CNNs require either multiple interconnected layers or recurrent connections to second nearest neighbors. However, most VLSI implementations of CNNs have included only a single layer of cells

1209

1210

B. E. Shi

interconnected with their nearest neighbors, since these are less complex but have a wide range of applications in image processing [Espejo et al., 1996; Paasio et al., 1999; Linan et al., 2002]. Work in implementations of two layer networks has been primarily motivated by models of visual processing in the retina or visual cortex [Boahen, 1997; Carmona et al., 2002; Shi, 1999]. Wu and Wen [2001] described the implementation of CNNs with larger neighborhood templates. This paper describes oriented spatial pattern formation observed in a CMOS implementation of a four layer CNN. To our knowledge, this chip contains the most recurrently interconnected twodimensional layers of any VLSI CNN implementation to date. We developed this chip to implement silicon neurons with orientation selective receptive fields similar to those observed in the mammalian visual cortex [Choi et al., 2004]. When the spatially uniform equilibrium point in this network is stable, the network does exhibit the expected orientation tuned receptive fields in response to an externally supplied input. The preferred orientation and the spatial scale of the receptive fields are adjustable via external bias voltages. Although many biological models assume that orientation selectivity arises from the pattern of feedforward connections from the lateral geniculate nucleus to the cortex [Hubel & Weisel, 1962; Ferster & Miller, 2000], the chip achieves the selectivity via recurrent interconnections between neighboring neurons. From an engineering point of view, the purely local interconnections required by the feedback network are easier to implement on chip than the much larger radius of connections required by a feedforward model. Because the connections are recurrent, the receptive fields of the neurons can still extend over a large spatial radius. These recurrent interconnections also raise the possibility of network instability, which gives rise to the strongly spatially oriented patterns of activity even in the absence of any external input that we describe here. In this context, it is interesting to note that patterns of activity in simple neural network models of the visual cortex have been hypothesized to underlie visual hallucinations [Oster, 1970; Markus, 1992]. Our mathematical analysis of this network demonstrates that the pattern formation dynamics of this network can be modeled by a network containing two diffusively coupled layers which

form patterns using a similar mechanism as described in [Turing, 1952] for reaction–diffusion equations. However, whereas spatial discretizations of reaction–diffusion equations include only interconnections between corresponding cells in different layers, this network includes interconnections between cells in one layer and nearest neighbors of corresponding cells in other layers. These additional connections give rise to a convection component in the differential equations governing the dynamics, which is responsible for the strongly oriented nature of the observed patterns. We demonstrate that the experimentally observed patterns from the network are in excellent agreement with our mathematical predictions. By systematically varying the bias voltages of the network that determine the intra-layer and inter-layer couplings, we demonstrate that it is possible to vary both the spatial orientation and spatial scale of the patterns. In addition, for horizontally and vertically oriented patterns, it is possible to vary the spatial bandwidth of the patterns in the direction perpendicular to the dominant orientation, which controls the “coherence” of the oriented bands. The remainder of this paper is organized as follows, Sec. 2 introduces the mathematical model of the array and its CMOS implementation. Section 3 presents a mathematical analysis of the pattern formation dynamics of the CNN, as well as measurement results from the array that are consistent with it. Section 4 concludes with a short summary.

2. Network Architecture and Implementation This section introduces the mathematical description of the network and briefly describes the transistor level circuit design and implementation. Choi et al. [2004] provided a more detailed description of the design and implementation, as well as experimental measurements of the receptive field profiles when the network is tuned to model orientation tuned cortical neurons.

2.1. Network architecture The chip implements a four-layer Cellular Neural Network whose layers are indexed by k ∈ {e+, e−, o+, o−}. Each layer consists of twodimensional array of cells, each with a non-negative state xk (m, n) and output yk (m, n). We group the layers into two pairs that we refer to as EVEN (e+

1211

Oriented Spatial Pattern Formation in a CMOS CNN

and e−) and ODD (o+ and o−). We choose these names because when the chip is tuned to model orientation selective visual cortical neurons, the cells in these layers have even and odd symmetric spatial receptive fields. The two layers within each pair are referred to as ON (+) and OFF (−). Although each cell in the chip also has an associated input uk (m, n), we consider only an autonomous CNN with zero input in the following. The output of each cell is a nonlinear memoryless function of the state of corresponding cells in each pair of layers. ye+ (m, n) = f+ (xe+ (m, n), xe− (m, n)) yo+ (m, n) = f+ (xo+ (m, n), xo− (m, n)) ye− (m, n) = f− (xe+ (m, n), xe− (m, n)) yo− (m, n) = f− (xo+ (m, n), xo− (m, n)) where f+ (x+ , x− ) and f− (x+ , x− ) are positive and negative half-wave rectified versions of x + −x− that saturate at Imax : f+ (x+ , x− ) =

|x+ − x− | − |x+ − x− − Imax | + Imax z

f− (x+ , x− ) = f+ (x− , x+ )

Ase+e+ = Ase−e− = Aso+o+ = Aso−o− = A∇2 Ae+o− = Ae−o+ = Ao+e+ = Ao−e− = A−1 Ae+o+ = Ae−o− = Ao+e− = Ao−e+ = A+1 where

(1)

In analogy with ON and OFF neurons in biology, the ON and OFF outputs encode positive and negative components of a signal: here the difference between the ON and OFF states in each pair of layers. Although this representation doubles hardware requirements, it has compelling advantages for both biological and engineering systems, such as efficient use of energy and communication bandwidth. The state of each layer evolves according to ! X X dxk s Akl ? yl + I = xk • −xk + Akl ? xl + dt l

in the standard autonomous CNN the time derivative of each cell’s state is a linear combination of the states and outputs of the cells in a small neighborhood of the cell. In this network this linear combination is also scaled by the state. Second, it adds a state feedback template Askl , which enables each cell’s state to influence its neighboring cells’ states directly, rather than through the output nonlinearity. Third, the output of each cell is a nonlinear function of the states in cells of two layers, rather than just one. Thus, the output nonlinearity introduces an additional coupling between layers beyond that dictated by the cloning templates. We have introduced these changes because they result in a network which is easily implemented using transistors operating in weak inversion. The nonzero cloning templates implemented by this chip are:

l

(2)

where the summations are over all layers (indexed by l) and I is a small constant bias current applied to all cells. The ? notation denotes correlation over a small neighborhood of each pixel (m, n), e.g. X (Askl ? xl )(m, n) = Askl (o, p)xl (m + o, n + p) . o,p

The coefficient matrices Askl , Akl and Bkl are the state feedback, output feedback and feedforward cloning templates. This network differs from the classical multilayer CNN [Chua & Yang, 1988] in three ways. First,



α1n −(2α1m + 2α1n )

 0 α1m  

0

0

α2n

0

 0

0

 A∇2 =  α1m 0 

0

 A−1 =  α2m 0

α2n

α1n

0





0

 A+1 =  0

0

0

0 0

0



 α2m  0

In each template matrix, the central element indicates the (0, 0) term. The α-parameters, α 1m , α1n , α2m , and α2n are non-negative reals and determine the strength of the interconnections between cells. Making these substitutions, we expand Eq. (2) to dxe+ = xe+ • (−xe+ + A∇2 ? xe+ dt + A+1 ? yo+ + A−1 ? yo− ) dxe− = xe− • (−xe− + A∇2 ? xe− dt + A−1 ? yo+ + A+1 ? yo− ) dxo+ = xo+ • (−xo+ + A∇2 ? xo+ dt + A−1 ? ye+ + A+1 ? ye− ) dxo− = xo− • (−xo− + A∇2 ? xo− dt + A+1 ? ye+ + A−1 ? ye− )

(3)

1212

B. E. Shi

u o-

u oo-

uo+

o-

o-

uo+

uo+

o+

o+

o+

e-

e-

e-

ue-

u ee+

u e+

is diffusively coupled to its nearest neighbor cells in the same layer, where the diffusion constants in the horizontal and vertical directions are controlled by α1m and α1n . The A+1 template determines the coupling to a target cell from cells above and to the right of the corresponding cell in another layer. Similarly, the A−1 template determines the coupling from cells below and to the left. The parameters α2m and α2n control the strength of the inter-layer coupling in the horizontal and vertical directions.

u o-

u ee+

e+

u e+

2.2. Circuit design and implementation

u e+

The circuit uses transistors operating in subthreshold or weak inversion. We represent each state variable array as the drain currents in an array of NMOS transistors, which we refer to as state transistors. The state transistors share the same gate fixed voltage V1 and operate in saturation. We assume that the continuous time dynamics arise due to parasitic capacitances associated with the source node. The multiplication by the state in Eq. (3) arises when we differentiate the exponential relationship between the source voltage and the drain current of a transistor operating in weak inversion.

Fig. 1. Cell interconnections in a 1D version of the four-layer network described by Eq. (2). Interconnections terminating with filled circles are inhibitory. Interconnections terminating with open triangles are excitatory. Interconnections by resistors are diffusive. The mutually inhibitory interconnections between positive and negative layers of the same type are introduced by the output nonlinearity.

Figure 1 shows the interconnections between cells 2 The A∇2 template is a discrete apforFigure a 1D array. proximation to the Laplacian operator. Each cell

x k(m, n + 1) V1

V1

V 1m

x k(m, n) V1

V1 V 1m

V 1n

V 1n

Fig. 2. Two-dimensional diffuser network. The state of each layer xk is represented by the drain currents through the vertical transistors. Figure 2The transistors connecting pairs of nodes implement diffusive intra-layer connections.

21

Oriented Spatial Pattern Formation in a CMOS CNN

V bq y+

x-

y-

x+

y

out

Fig. 3. The circuit mapping the state variables, represented by currents, to outputs, represented by currents. All transistors operate in saturation.

α 2m y

V2 V 2m

V 2n

Fig. 4. The tilted current mirror circuit, which implements the inter-layer connections and mirrors the output of each cell.

As shown in Fig. 2, we implement diffusive interconnections with single NMOS transistors, called diffusers, that connect the source terminals of neighboring state transistors [Andreou & Boahen, 1994; Vittoz & Arreguit, 1993]. By varying the gate voltages of the diffusers, we can change the intra-layer coupling: α1m = eκ(V1m −V1 )/UT

α1n = eκ(V1n −V1 )/UT

where UT is the thermal voltage and 0 < κ ≤ 1 is a process dependent constant. At the edge of the array, diffusers which would connect to cells outside the array are simply deleted, implementing reflective or zero-flux boundary conditions. We compute the output nonlinearity by connecting the drains of corresponding state transistors in each pair of ON/OFF layers to the x + and x− inputs of the circuit shown in Fig. 3 [Zaghloul, 2001]. Raising the control voltage V bq closer to VDD decreases both the current at which the output saturates as well as a small quiescent output current that flows in the absence of any input current. For the results reported here VDD − Vbq = 0.38 V.

23

We implement inter-layer coupling and sense the output of each cell using the tilted current mirror circuit shown in Fig. 4. We connect each output of the ON/OFF nonlinearity circuit in Fig. 3 to the input y of one of these circuits. The source voltage of the diode connected transistor is held constant at V2 . The four outputs with gain α2m and α2n are connected to the source terminals of the state transistors in target cells in other layers. The source terminal of each state transistor is connected to one output from four different tilted current mirror circuits. The bias current models the current through these outputs due to the small quiescent output of the ON/OFF nonlinearity circuit as well as leakage. We change the strength of the connection by adjusting the source voltages V2m and V2n : α2m = e(V2 −V2m )/UT

α 2n y

1213

α2n = e(V2 −V2n )/UT .

We connect the output of the current mirror circuit with unity gain (labeled “out”) to the input of a spiking neuron circuit that outputs a spike train whose rate is proportional to the input current [Culurciello et al., 2003]. Spikes from the neurons in the array are transmitted off chip using the Address Event Representation (AER) communication protocol [Boahen, 2000]. The chip contains four layers with 32 × 64 cells each, over 8000 cells in total. It was fabricated in the TSMC0.25um mixed signal/RF process available through MOSIS. This process contains 5 metal layers and 1 poly layer and uses nonepitaxial wafers. We construct the layout by tiling meta-cells, whose layout is shown in Fig. 5(a). One meta-cell contains the circuits for one cell from each of the four layers as well as four spiking neurons for output and four current mode integrators for input. It contains 172 transistors (48 to implement the analog network and 124 to implement the input/output circuits) and occupies 52 µm by 49 µm. The entire chip layout measures 3.84 mm by 2.54 mm. It contains the CNN array, digital communication circuits, and input/output pads at the periphery for external connections. Figure 5(b) shows a photograph of the packaged and unpackaged chips. Figure 5(c) shows a photomicrograph of the fabricated chip.

3. Analysis and Experimental Results In this section we analyze the four-layer network described above, predicting the dependence of the

1214 Figure B. E. Shi6

Figure 6 Figure 6

(a) (a)

(b) (b)

(a)

(b) (a)

(b)

(c)

Figure 6

(c) (c)

Fig.Figure 5. (a)6The layout of one meta-cell, which is tiled to construct (c)the entire CNN array. The layout is flipped vertically in alternating rows. (b) A photograph of the packaged chip and the unpackaged silicon chip. The minimum division on the ruler 6 is 1 mm.Figure (c) A photomicrograph of the fabricated chip. 26

26

Oriented Spatial Pattern Formation in a CMOS CNN

spatial patterns on the α parameters. Afterwards, we validate these predictions with measurement results from the chip.

3.1. Pattern formation dynamics To analyze the pattern formation dynamics of this network, we examine its evolution around the spatially uniform equilibrium point where the states of all cells are equal to the bias current I. By making a coordinate transformation where the activities in ON and OFF layers are expressed as their difference and sum, we show that the resulting patterns are primarily determined by the dynamics of the interaction between the difference components of the EVEN and ODD layers. These patterns show a strong spatial orientation and scale preference, both dyo+ = dxo+



dyo+ = dxo−



1 0 −1 0

for x ˜o+ > x ˜o− otherwise for x ˜o+ > x ˜o− otherwise

1215

of which are can be controlled by adjusting the external bias voltages that set the α parameters. Because the derivatives of f+ and f− are not continuous at x+ = x− = I, we make a piecewise linear approximation. Define x ˜ k to be a small perturbation about the operating point. For k = e + , we have  d˜ xe+ ≈I −x ˜e+ + A∇2 ? x ˜e+ dt   dyo+ dyo+ x ˜o+ + x ˜o− + A+1 ? dxo+ dxo−   dyo− dyo− x ˜o+ + x ˜o− + A−1 ? dxo+ dxo− where all derivatives are evaluated at the operating point xe+ = xe− = xo+ = xo− = I. Thus, dyo− = dxo+ dyo− dxo−

−1 0  1 = 0 

for x ˜o+ < x ˜o− otherwise for x ˜o+ < x ˜o− otherwise

Defining h+ and h− to be positive and negative half-wave rectifiers, h+ (x+ , x− ) = 0.5(x+ − x− ) + 0.5|x+ − x− |

h− (x+ , x− ) = 0.5(x− − x+ ) + 0.5|x− − x+ |

we can write d˜ xe+ ≈ I(−˜ xe+ + A∇2 ? x ˜e+ + A+1 ? h+ (˜ xo+ , x ˜o− ) + A−1 ? h− (˜ xo+ , x˜o− )) . dt The saturation at Imax does not affect the dynamics in the region around the uniform equilibrium point. Applying a similar procedure to the other layers, we obtain 1 I 1 I 1 I 1 I

dxe+ dt dxe− dt dxo+ dt dxo− dt

= −xe+ + A∇2 ? xe+ + A+1 ? h+ (˜ xo+ , x ˜o− ) + A−1 ? h− (˜ xo+ , x ˜o− ) = −xe− + A∇2 ? xe− + A−1 ? h+ (˜ xo+ , x ˜o− ) + A+1 ? h− (˜ xo+ , x ˜o− )

= −xo− + A∇2 ? xo− + A+1 ? h+ (˜ xe+ , x ˜e− ) + A−1 ? h− (˜ xe+ , x ˜e− )

To facilitate analysis of the network dynamics, we make a change of variables. Define the difference and sum components of the state variables to be x ˜ed = x ˜e+ − x ˜e− and

(4)

= −xo+ + A∇2 ? xo+ + A−1 ? h+ (˜ xe+ , x ˜e− ) + A+1 ? h− (˜ xe+ , x ˜e− )

x ˜es = x ˜e+ + x ˜e−

(5)

and similarly for the odd state variables. Applying this change of variables to Eq. (4), we

obtain: 1 d˜ xed I dt xod 1 d˜ I dt xes 1 d˜ I dt xos 1 d˜ I dt

= −˜ xe+ + (A∇2 ? x˜ed ) + (A∇ ? x ˜od ) = −˜ xod + (A∇2 ? x ˜od ) − (A∇ ? x ˜ed ) = −˜ xes + (A∇2 ? x ˜es ) + (AΠ ? x ˜od ) = −˜ xos + (A∇2 ? x˜os ) − (AΠ ? x˜ed )

(6)

1216

B. E. Shi

where 



0 α2n 0   0 α2m  A∇ =  −α2m 0 −α2n 0   0 α2n 0   0 α2m  AΠ =  α2m 0 α2n 0

and |˜ x| denotes the element-wise absolute value of x ˜. The difference components evolve independently of the sum components. The sum components are driven by the difference components. The dynamics of the sum and difference components are both linear. The nonlinearity is confined to the coupling from the difference components to the sum components. Correlation by A∇ is a discrete approximation to a directional derivative in the direction (2α2m , 2α2n ). Thus, the EVEN and ODD difference components are cross-coupled by a convection term. We see below that the strength and direction of the directional derivative control the existence of patterns, as well as their orientation. The A Π template is a combination of even impulse pairs in the horizontal and vertical directions. The pattern formation properties of this array can be understood by examining the evolution of the spatial Fourier components, which evolve independently because the templates are space invariant. Assume a doubly infinite array and define ˜ k (ωm , ωn ) to be the 2D discrete Fourier transforms X of the state perturbation, e.g. X ˜ k (ωm , ωn ) = X x ˜k (m, n)e−jωm m−jωn n . (7) m,n

The Fourier components of the sum components are stable for any choice of α parameters, since they evolve according to a discrete approximation to a lossy diffusion equation. If we take the discrete Fourier Transform of the last two equations in Eq. (6), the correlations by A∇2 and AΠ correspond to multiplication by A˜∇2 (ωm , ωn ) = −2α1m (1 − cos ωm ) − 2α1n (1 − cos ωn ) and A˜Π (ωm , ωn ) = 2(α2m cos ωm + α2n cos ωn ) .

To avoid clutter in the following, we simply assume the dependence on ωm and ωn . Thus, " # " #" # ˜ es ˜ es −1 + A˜∇2 0 X 1 d X = ˜ os ˜ os I dt X 0 −1 + A˜∇2 X # #" " ˜ |ed| X 0 A˜Π + ˜ |od| X A˜Π 0 ˜ |ed| (X ˜ |od| ) denotes the discrete Fourier where X transform of the absolute value of the even (odd) difference component, which evolves independently of the sum components. Stability is ensured since A˜∇2 < 0 for all (ωm , ωn ) and all α1m , α1n ≥ 0. Thus, the final patterns must be dictated by the evolution of the difference component. Taking the discrete Fourier transform of the first two equations in Eq. (6), we obtain " # " #" # ˜ ed ˜ ed −1 + A˜∇2 A˜∇ X 1 d X = ˜ od ˜ od I dt X −A˜∇ −1 + A˜∇2 X

(8)

where A˜∇ = j2(α2m sin ωm + α2n sin ωn ) and j = √ −1. By defining x ˜d = x ˜ed + j x ˜od , we can express Eq. (8) compactly as the complex valued equation: 1 d ˜ ˜d . Xd = (−1 + A˜∇2 − j A˜∇ )X (9) I dt The final pattern is dominated by the unstable Fourier components, i.e. those for which −1+ A˜∇2 − j A˜∇ > 0. If we define K = −1 + 2rm (1 − cos Ωm ) + 2rn (1 − cos Ωn ) p p rm = α21m + α22m rn = α21n + α22n (10)     α2n α2m Ω = atan Ωm = atan α1m α1n then, −1 + A˜∇2 − j A˜∇ = K − 2rm (1 − cos(ωm − Ωm )) − 2rn (1 − cos(ωn − Ωn )) , which reaches its maximum value of K for (ωm , ωn ) = (Ωm , Ωn ). If K is positive, then Fourier modes corresponding to spatial frequencies near (Ωm , Ωn ) will be unstable. For fixed α1m and α1n , K increases with α2m and α2n since rm , rn , Ωm and Ωn all increase. The magnitude of the unstable Fourier components in the initial state grows exponentially over time. Ultimately, their growth is limited by the saturation

Figure 5 1217

Oriented Spatial Pattern Formation in a CMOS CNN

of y+ and y− , which decreases the effective gain of the cross-coupling, α2m and α2n . In summary, our analysis predicts that the resulting patterns will be dominated by the complex exponential corresponding to the most unstable Fourier component:

10 rn

8 6

x ˜ed + j x ˜od ≈ ej(Ωm m+Ωn n+φ)

4

= cos(Ωm m + Ωn n + φ)

2

+ j sin(Ωm m + Ωn n + φ) where Ωm (Ωn ) is determined by Eq. (10) and φ is a phase parameter that depends upon initial conditions and the boundary conditions of the array. Because the ON and OFF outputs of the array are Figure 5 positive and negative half-wave rectified versions of the state, we should observe alternating regions of ON and OFF activity as we move in the direction θ = atan(Ωn /Ωm ). In addition, we expect activity in the EVEN layers to be 90◦ out of phase with activity in the ODD layers. The constraint that the α parameters be non-negative restricts the range of orientations to between 0 and π/2. However, excluded orientations can be obtained simply by reflecting one of the allowed orientations around either one or both of the coordinate axes. The remainder of this section illustrates these predictions with measured patterns of activity from the array.

3.2. Horizontal and vertical pattern formations We obtain horizontally oriented patterns if we let α2m ≈ 0 by tuning V2m − V2  UT . This simplifies Eq. (10) to K = −1 + 2rn (1 − cos Ωn ) p rm = α1m rn = α21n + α22n

Ωm = 0

Ωn = atan(α2n /α1n )

Patterns exist if K > 0, i.e. rn > (2(1 − cos Ωn ))−1 . Figure 6 shows the parameter region that exhibits pattern formation. Since both α1n and α2n are nonzero, Ωn lies between 0 and π/2. Figure 7 shows the measured output of the array tuned to exhibit horizontally oriented spatial patterns. The bias voltages used to obtain these patterns are given in the first line of Table 1. We represent the activity in the array using two pseudocolor images: one for the EVEN layers and one for the ODD layers. We represent the spike rates

0 0

0.1

0.2

0.3

0.4

0.5 Ωn ⁄ π

Fig. 6. The shaded areas show the regions exhibiting horizontal pattern formation in the Ωn − rn parameter space. Table 1. Pattern (a) (b) (c) (d) (e)

Bias voltages (V ) for patterns in Fig. 8.

V1

V1m

V1n

V2

V2m

V2n

1.500 1.500 1.500 1.500 1.500

1.520 1.404 1.477 1.572 1.450

1.488 1.488 1.439 1.508 1.450

0.195 0.195 0.195 0.195 0.195

0.443 0.443 0.192 0.102 0.215

0.193 0.193 0.427 0.427 0.215

on both the ON and OFF layers by a color assigned in the HSV (hue-saturation-value) color space. The value of all pixels is one, the maximum value. We encode changes in the spike rate as changes in the saturation, where lower spike rates correspond to less saturated (whiter) colors. We encode the relative amount of ON versus OFF activity as changes in the hue. For the EVEN images, the hue changes continuously from blue to red, with blue representing pure OFF activity, red representing pure ON activity and magenta representing equal amounts of ON and OFF activity. For ODD images, the hue changes continuously from magenta to yellow, with magenta representing pure OFF activity, yellow representing pure ON activity and red representing equal amounts of ON and OFF activities. Mathematically, if the ON and OFF spike rates are scaled and clipped so that a spike rate of 0 Hz maps to 0 and spike rates of 1 kHz and above map to 1, then   25 ON − OFF 1 c+ H= 6 max{ON, OFF} S = 0.2 + 0.8 max{ON, OFF} V=1

1218 B. E. Shi Figure 7

Figure Figure7 7

(a)

(b) (b)

(a)

(a) (a)

2

spike(kHz) rate (kHz) spike spike rate rate (kHz)

rate (kHz) spike spike rate (kHz) spike rate (kHz)

22 1.5 1.5 1.5 1

11 0.5 0.5 0.5 0

00

Figure 7

(b) (b)

2

22 1 11 0

00 -1

5

10

15

55

1010 (c) 1515

-1-1 -2

20

25 30 pixel 25 30 30 2020 25 pixel pixel

-2-2

(c)(c) (c)

5

10

15

55

10 (d) 15 10 15

20

20 20

25 30 pixel 25 30 30 25 pixel pixel

(d) (d) (d)

Fig. 7. A horizontally oriented spatial pattern measured from the array. (a) The output activity in the EVEN layers. Blue Figure Figure 7 7 activity. Red represents ON activity. The saturation represents the spike rate, with paler (whiter) colors represents OFF corresponding to lower spike rates. (b) The output activity in the ODD layers. Magenta represents OFF activity. Yellow represents ON activity. The saturation represents the spike rate. (c) A cross-section of the activity in the EVEN array along column 48. The blue trace represents OFF activity. The red trace represents ON activity. (d) A cross-section of the difference in ON and OFF activity in the EVEN (blue) and ODD (magenta) arrays.

where c = 5 for EVEN activity and c = 0 for ODD activity. This color map assigns every combination of scaled and clipped activities to a unique color. However, because the ON/OFF nonlinearity excludes simultaneous ON and OFF activity, the images are dominated by the blue/red or magenta/yellow hues at varying saturations, with a few pale magenta or pale red pixels representing locations where both the ON and OFF activities are close to zero. The vertical cross-section of activity through the EVEN pair of layers plotted in Fig. 7(c) shows clearly alternating regions of ON and OFF activities. Figure 7(d) shows the difference between ON and OFF activities through the same vertical

cross-section in the EVEN and ODD layer pairs. This plot clearly shows the expected 90 ◦ spatial phase offset between the activity in the two pairs of layers. The coherence of the horizontal bands in the pattern is determined by V1m . The unstable modes are contained within a region in the (ω m , ωn ) plane with boundary defined by 2rm (1 − cos ωm ) + 2rn (1 − cos(ωn − Ωn )) = K

For small ω, 2(1 − cos ω) ≈ ω 2 and the boundary can be approximated by the oval defined by 2 rm · ω m + rn · (ωn − Ωn )2 = K

which is centered at (0, Ωn ) and has axes parallel to the ωm and ωn coordinate axes. The ωm axis has 27

2727

Figure 8

Oriented Spatial Pattern Formation in a CMOS CNN

1219

Figure 8

Figure 8 Figure 8

(a)

(b)

(a)

(b)

(a)

(b)

(a)

(b)

(a)

(b)

(c)

(d)

2

(c)

(c)

(d)

ωn (radians/pixel) ωn (radians/pixel) ωn (radians/pixel) ωn (radians/pixel)

(c) (c)

Figure 8

2 1

2 1

1 0 2

0 -1 1 0 -1 -2 0-2 -1 -2 -2 -1 -2 -2

(e)

Figure 8

(d)

(e)

-2 -2

(e) (e)

(d) (d)

-1 0 1 ωm (radians/pixel) -1 1 (f) 0 ωm (radians/pixel) -1 1 (f)0 ωm (radians/pixel) -1 (f) (f)0 1 ωm (radians/pixel)

2 2

2 2

Fig. 8. Measured spatial patterns with varying orientation and spatial frequency. We show only activity in the EVEN array, as Figure 8 (f) OFF activity. Red represents the activity in the OFF array is (e) qualitatively similar except for a 90◦ phase shift. Blue represents ON activity. The saturation represents the spike rate, with paler (whiter) colors corresponding to lower spike rates. (a) A Figure 8oriented pattern with strong coherence of the bands of activity in the horizontal direction. This figure is the same horizontally as that in Fig. 7(a), but is repeated here to facilitate comparison. (b) A horizontally oriented pattern with weak coherence in the horizontal direction. (c) A vertically oriented pattern with high spatial frequency. (d) A vertically oriented pattern with low spatial frequency. (e) A diagonally oriented spatial pattern. (f) A graphical representation of the mean and variance of the spectral energies in the patterns shown here. Each oval is centered at the mean of the spectral energy distribution, denoted by an “×”. Its size represents the variance. Each oval corresponds to one of the patterns in this figure: (a) green, (b) blue, (c) red, (d) cyan and (e) magenta.

28

28

1220

B. E. Shi

length Lm = 2

s

Table 2. Spectral energy mean and variances for patterns in Fig. 8.

2rn (1 − cos Ωn ) − 1 rm

which decreases with increasing rm = α1m (increasing V1m ). A small axis length leads to patterns that exhibit strong horizontal banding because the small range of unstable frequencies leads to patterns which change slowly in the horizontal direction. As Lm increases (V1m decreases), the coherence of the horizontal bands decreases. We can observe this phenomenon comparing the two patterns of EVEN activity in Figs. 8(a) and 8(b), which were measured with the array tuned to the same bias voltages except for a change in V 1m . Table 1 lists the bias voltages used. We can measure the coherence quantitatively by computing the mean and variance of the energy spectra of the two patterns. Let Y (k, l) for k = 0, . . . , M and l = 0, . . . , N be the discrete Fourier transform of the complex valued output difference: y = (ye+ − ye− ) + j(yo+ − yo− ). Define X X  2πk  |Y (k, l)|2 ωm = M |Y |2Total k

ωn =

X X  2πl  |Y (k, l)|2 k

Kmm

l

N

|Y |2Total

X X  2πk 2 |Y (k, l)|2 − (ω m )2 = M |Y |2Total k

Knn

l

l

X X  2πk  2πl  |Y (k, l)|2

− (ω m )(ω n ) M N |Y |2Total P P where |Y |Total = k l |Y (k, l)|2 . Table 2 lists the numerical values of the mean and variance. As predicted, the values computed for the two patterns are nearly identical, except that K mm is larger for Fig. 8(b) than for Fig. 8(a). Figure 8(e) represents the same information in Table 2 graphically. Each set of five numbers is represented by the oval determined by the equation    Kmm Kmn ωm − ω m =1 [ωm − ω m ωn − ω n ] Kmn Knn ωn − ω n Kmn =

k

(a) (b) (c) (d) (e)

ωm

ωn

Kmm

Knn

Kmn

0.011 0.018 1.073 0.439 1.142

1.022 1.036 −0.014 −0.023 1.0326

0.377 1.493 0.238 0.625 0.394

0.333 0.319 0.714 0.840 0.419

0.001 0.005 0.014 0.007 0.045

parameters K. For example if Kmn = 0, the oval has axes which are aligned with and ωn co√ the ω m √ ordinate axes with length Kmm and Knn . We see that the green and blue ovals corresponding to the patterns in Figs. 8(a) and 8(b) are centered at the same point with approximately the same vertical height, but that the blue oval is considerably wider horizontally.

3.3. Vertical pattern formation In exact analogy to the previous discussion, we can obtain vertically oriented patterns by fixing α 2n ≈ 0 (V2n −V2  UT ). Figures 8(c) and 8(d) show two examples of vertical pattern formation using the bias voltages listed in Table 1. As we can see in Table 2 and Fig. 8(f), the mean of the energy spectrum for the pattern in Fig. 8(d) is at a much lower spatial frequency than the mean for the pattern in Fig. 8(c) which is consistent with the broader bands observed in that pattern.

l

X X  2πl 2 |Y (k, l)|2 = − (ω n )2 N |Y |2Total k

Pattern

l

which is centered at (ω m , ω n ) and whose size and orientation represent the variance and covariance

3.4. Diagonal pattern formation By setting both α2m > 0 and α2n > 0, we obtain diagonally oriented patterns, as shown in Fig. 8(e). To obtain these patterns, we set V1m = V1n and V2m = V2n , so that α1m ≈ α1n and α2m ≈ α2n . This implies that Ωm ≈ Ωn . Consistent with this prediction, Table 2 and Fig. 8(f) show that ω m ≈ ω n . Although we do not show data here, we can obtain patterns at orientations besides 0 ◦ , 45◦ and 90◦ by appropriate choice of bias voltages.

4. Conclusions We have described and characterized spatial pattern formation observed in a CMOS implementation of a four-layer CNN. This chip is one of the most complex CNNs reported to date, both in terms of the number of layers and the total number of cells. One of the unique features of the patterns

Oriented Spatial Pattern Formation in a CMOS CNN

we observe from this network is a strong preference for a given orientation and spatial scale. In contrast, patterns observed in more commonly studied reaction–diffusion equations often demonstrate a strong preference for a given spatial scale, but little orientation preference. VLSI implementations of CNNs offer an easily controlled and easily characterized physical substrate to study complex nonlinear spatiotemporal dynamics in real time. Indeed, our experimental results are in excellent agreement with the predictions from our mathematical analysis. Our analysis indicates that the inter-layer coupling in the network gives rise to a convection term in the equations governing the pattern formation dynamics. The balance between the strengths of the convective inter-layer coupling and the diffusive intra-layer coupling determines the spatial scale and orientation of the observed patterns. By controlling external bias voltages that determine these coupling strengths in the chip, we have been able to measure patterns at the expected orientations and spatial scales.

Acknowledgments The author would like to acknowledge Thomas Choi, Kwabena Boahen, Paul Merolla, John Arthur, Kai Hynna and Brian Taba for their assistance in designing and fabricating the chip described here, and S. F. Luk for his assistance in preparing the figures. This work was supported by the Hong Kong Research Grants Council under Grant HKUST6218/01E.

References Andreou, A. G. & Boahen, K. A. [1994] “Neural information processing II,” in Analog VLSI: Signal and Information Processing, eds. Ismail, M. & Fiez, E. (McGraw-Hill, NY). Boahen, K. A. [1997] “The retinomorphic approach: Pixel-parallel adaptive amplification, filtering and quantization,” Anal. Integr. Circuits Sign. Process. 13, 53–68. Boahen, K. A. [2000] “Point-to-point connectivity between neuromorphic chips using address events,” IEEE Trans. Circuits Syst.-II: Anal. Dig. Sign. Process. 47, 416–434. Carmona, R., Jimenez-Garrido, F., Dominguez-Castro, R., Espejo, S. & Rodriguez-Vazquez, A. [2002] “CMOS realization of a 2-layer CNN universal machine chip,” in 7th IEEE Int. Workshop on Cellular Neural Networks and their Applications, pp. 444–451.

1221

Choi, T. Y. W., Shi, B. E. & Boahen, K. A. [2004] “An ON–OFF orientation selective address event representation image transciever chip,” IEEE Trans. Circuits Syst.-I 51, 342–353. Chua, L. O. (ed.) [1995] “Special issue: On nonlinear waves, patterns and spatio-temporal chaos in dynamic arrays,” IEEE Trans. Circuit Syst.-I: Fund. Th. Appl. 42. Chua, L. O. [1998] CNN: A Paradigm for Complexity (World Scientific, Singapore). Chua, L. O. & Yang, L. [1988] “Cellular neural networks: Theory,” IEEE Trans. Circuits Syst. 35, 1257–1272. Culurciello, E., Etienne-Cummings, R. & Boahen, K. A. [2003] “A biomorphic digital image sensor,” IEEE J. Solid-State Circuits 38, 281–294. Hubel, D. H. & Wiesel, T. N. [1962] “Receptive fields, binocular interaction and functional architecture in the cat’s visual cortex,” J. Physiol. 160, 106–154. Espejo, S., Carmona, R., Dominguez-Castro, R. & Rodriguez-Vazquez, A. [1996] “A CNN universal chip in CMOS technology,” Int. J. Circuit Th. Appl. 24, 93–109. Ferster, D. & Miller, K. D. [2000] “Neural mechanisms of orientation selectivity in the visual cortex,” Ann. Rev. Neurosci. 23, 441–471. Linan, G., Espejo, S., Dominguez-Castro, R. & Rodriguez-Vazquez, A. [2002] “ACE4k: An analog I/O 64 ∗ 64 visual microprocessor chip with 7-bit analog accuracy,” Int. J. Circuit Th. Appl. 30, 89–116. Markus, M. [1992] “Neuronal modelling of hallucinations,” in World of Consciousness, 1st Int. Congress of the European College for the Study of Consciousness (ESCS), eds. Schlickting, M. & Laerner, H., pp. 131–140. Oster, G. [1970] “Phosphenes,” Sci. Amer. 22, 82–87. Paasio, A., Kananen, A., Halonen, K. & Porra, V. [1999] “A QCIF resolution binary I/O CNN-UM chip,” J. VLSI Sign. Process. 23, 281–290. Shi, B. E. [1999] “Focal plane implementation of 2D steerable and scalable cortical filters,” J. VLSI Sign. Process. 23, 319–334. Turing, A. M. [1952] “The chemical basis of morphogenesis,” Philos. Trans. Roy. Soc. B237, 37–72. Vittoz, E. & Arreguit, X. [1993] “Linear networks based on transistors,” Electron. Lett. 29, 297–299. Wu, C. Y. & Yen, W. C. [2001] “A new compact neuronbipolar junction transistor (ν BJT) cellular neural network (CNN) structure with programmable large neighborhood symmetric templates for image processing,” IEEE Trans. Circuits Syst.-I: Fund. Th. Appl. 48, 12–27. Zaghloul, K. A. [2001] A Silicon Implementation of a Novel Model for Retinal Processing, Ph.D. thesis, University of Pennsylvania.