Distributed Memory Parallel Architecture Based on Modular Linear Arrays for 2-D Separable Transforms Computation.
A framework for mapping systematically 2-dimensional (2-D) separable transforms into a parallel architecture consisting of fully pipelined linear array stages is presented. The resulting model architecture is characterized by its generality, high degree of modularity, high throughput, and the exclusive use of distributed memory and control. There is no central shared memory block to facilitate the transposition of intermediate results, as it is commonly the case in row-column image processing architectures. Avoiding shared central memory has positive implications for speed, area, power dissipation and scalability of the architecture. The architecture presented here may be used to realize any separable 2-D transform by only changing the coefficients stored in the processing elements. Pipelined linear arrays for computing the 2-D Discrete Fourier Transform and 2-D separable convolution are presented as examples and their performance is evaluated.
Introduction
Separable transforms play a fundamental role in digital signal and image processing. Nearly every
problem in DSP is based on the transformation of a time or space domain signal to alternative
spaces better suited for efficient storage, transmission, interpretation or estimation. The most
commonly employed 2-dimensional (2-D) signal transforms, such as the Discrete Cosine, Fourier,
Sine Transforms (DCT, DFT, and DST) and the Discrete Hough and Radon Transforms (DHT,
DRT), are known as separable due to a special type of symmetry in their kernel [1]. Separable
transforms are computationally less expensive than non-separable ones, having time complexity in
O(M) per output data element, as opposed to O(M 2 ) for non-separable forms, when applied to a
size M \Theta M input image.
The majority of the available VLSI implementations for separable transforms are based on the
popular row-column approach (see for example Rao and Yip [2], Guo et al. [3], Y.-P. Lee et al. [4],
Bhaskaran and Konstantinides [5], Gertner and Shamash [6], and Chakrabarti and J'aj'a [7]), where
the 2-D transform is performed in three steps: (i) 1-D transformation of the input rows, followed
by (ii) intermediate result transposition (usually implemented with a transposition memory), and
(iii) 1-D transformation of the intermediate result columns, as illustrated in Fig. 1 (a). In this
figure, and throughout this paper, we assume that all image inputs become available in raster-scan
order, that is, as a continuous sequence of row vectors. Arrays that can accept and produce data
in raster-scan avoid expensive host interface memory buffering.
Each of the 1-D row and column processing blocks in Fig. 1 (a) may be realized either with highly
modular structures, such as standard pipelined arrays of identical Processing Elements (PEs), or
by some other type of architecture optimized for the targeted transform. However, one of the
main shortcomings of the conventional row-column architecture is that the central memory block
severely limits the modularity of the overall structure, despite the possible modularity of the 1-
D row and column arrays. In addition, central memory requires moderately complex address
generation circuitry, row and column decoders, sense amplifiers depending on its size, and control
for coordinating the concurrent access of the two row and column architectures to the shared
memory. Large central memories also have negative implications for low-power design.
In this paper we present a method for synthesizing fully pipelined modular arrays for 2-D
separable transforms. The derived computational structures use only small-size FIFO memories
local to the PEs which require neither address generation nor central control and memory-to-array
routing. Faster clock speed, smaller area, and low power can be achieved with architectures having
a regular set of simple PEs with distributed memory and control. Reducing the size of memory
structures has been shown to reduce power consumption [8], largely due to reduced effective bus
capacitances of smaller memories, and absence of sense amplifiers. The architectures synthesized
are general in the sense that any separable transform can be realized by programming in the PEs
memory
bank for
transposition
address
row processing
(a)
row column
raster-
scan
input
(b)
logic
column processing

Figure

1: (a) Conventional architecture for raster-scan 2-D separable transforms relying on shared
central memory bank for the transposition of intermediate results, and (b) proposed alternative
architecture based on two linear arrays where the central memory bank is replaced by FIFO queues
distributed to the PEs.
the appropriate set of kernel coefficients.
The benefits of eliminating memory transposition have been recognized by Lim and Swartzlander
[9] and by Wang and Chang [10]. In [9], the authors treat the DFT problem in the context of
multi-dimensional systolic arrays. Our method is distinct from [9] primarily in that it is restricted
to single input (SI) linear arrays using raster-scan I/O. Fisher and Kung [11] have shown that linear
arrays have bounded clock skew regardless of size, whereas higher than 2-D arrays of arbitrary size
may not be technologically feasible. (See Lee and Kedem [12] for a discussion on general aspects of
mapping algorithms to linear arrays.) Although the DCT arrays presented in [10] avoid transposition
by means of local PE storage, the resulting architectures require M 2 multipliers on an M \Theta M
input, which in many cases may be prohibitive. The arrays developed in this paper require only
multipliers.
Another important characteristic of the proposed architecture is that it uses localized communications
and its interconnection complexity is in O(M) (for M \Theta M separable 2-D transforms).
To the best of our knowledge, parallel architectures that may use less than O(M) PEs also require
O(M) input/output ports and may exhibit non-local communications among PEs i.e. their inter-connection
complexity may grow faster than O(M ). Furthermore, the reduction in the number of
PEs becomes possible by exploiting the specific coefficient symmetries, i.e. these architectures may
not compute any desirable 2-D separable transform, something that the architecture developed here
is capable of.
In addition to the general computational structure, we also derive here as examples arrays for
the 2-D DFT and for 2-D separable convolution. Although purely from a computational efficiency
point of view the DFT is more expensive than the Fast Fourier Transform (FFT), from a VLSI
implementation point of view there are significant reasons why the DFT may be preferable to the
FFT (particularly with a small number of coefficients). The FFT's complex data routing limits
overall speed and is expensive in terms of chip area (see Swartzlander [13] and Thompson [14]). In
this paper we use the 1-D DFT array by Beraldin et al. [15] and derive a modular array structure
for computing the 2-D DFT without transposition memory.
The 2-D separable convolution has been studied extensively for many years (see for example
Dudgeon and Mersereau [16]) and has been used as a basis for constructing general transfer
functions. Abramatic et al. [17] has used separable FIR systems for constructing non-separable
IIR filters, and Treitel and Shanks [18] have used parallel separable filters for approximating non-separable
ones. Despite the limited range of frequency responses of separable filters, they are
attractive due to their low computational cost as compared to non-separable systems. In addi-
tion, separable convolution plays a central role in the Discrete Wavelet Transform (DWT) [19], as
currently most of the basis used for the 2-D DWT are separable.
The rest of this paper is organized as follows. In Section 2 we fix notation, define general
separable 2-D transforms, and discuss alternative views of conventional row-column 2-D processing
that will allow us to formulate equivalent algorithms suitable for raster scan I/O. In Section 3 the
derived algorithms are then transformed systematically to fully pipelined linear array structures
with distributed memory and control using appropriate linear space-time mapping operators. In
Sections 4 and 5 we apply the general method developed in Section 3 to derive modular linear
arrays for the 2-D DFT, and for 2-D separable convolution, respectively.
Transforms, Definitions and Notation
A 2-D transform of an M \Theta M signal x(i; is given by
Notice that in a general 2-dimensional (2-D) transform there are M 2 , 2-D distinct orthonormal
basis functions spanning the transform space, namely g
(where g   denotes
the complex conjugate of g) [1, 20]. Separable and symmetric transforms have the property that
l;m (i; Hence one can express y(l; m) as
l (i)
In order to formulate the transforms in matrix form, we adopt the following notation: Let
A be an M i \Theta M k matrix of elements indexed from zero and interchangeably
also denoted as a(i; k) or a i;k . Furthermore, let column vectors M i \Theta 1 denoted by lower-case
letters as
where the superscript t denotes transposition. We adopt Golub
and Van Loan's [21] Matlab-like notation and denote submatrices as A(u; v), where u and v are
integer column vectors that pick out the rows and columns of A defining the submatrix. Also, index
vectors are specified by colon notation as implying
l
. For example,
is a 2 \Theta 5 submatrix containing rows 2 and 3, and columns 5 through 9 of matrix A.
Entire rows or columns can be specified by a single colon as A(:; k), which corresponds to A's kth
column, and A(i; :) denoting A's ith row. Index vectors with non-unit increments are specified as
denoting a count from p to r in length-q increments.
Let X and Y be the transform input and transform coefficient (output) M \Theta M matrices
respectively, and M \Theta M matrix G be the separable transformation kernel whose columns are
the M , 1-D distinct basis vectors such that G(i;
, then matrix
With this notation in place, Eqn. (2) can then be rewritten as
where
m). Consequently, a separable transform can be computed as
two consecutive matrix-matrix products. The first matrix product can be written as
and can be viewed as operating with matrix G on the rows of input X to produce the rows of the
intermediate result Z. The second matrix product,
can be viewed as operating with matrix G t on the columns of Z to produce the columns of Y . This
interpretation forms the basis for the traditional row-column processing method that necessitates
the transposition of Z. In Section 3 we show how Y a be obtained by processing directly the rows
of Z in the order in which they become available, thus eliminating the need for transposition. We
now give the definition for raster scan ordering.
Definition 2.1 Given the set of matrix elements g, we say
that element X(i; precedes in raster scan order element X(i 0
only if either i 0
Raster scan induces a total ordering on the input set, and is equivalent to row-major ordering.
Let us also define the vectors x and z as x j
and z j
which are column vector views of matrices X and
Z respectively, suitable for raster scan processing. means that the p th row of Z can be
derived by multiplying the p th row of X times G, i.e. Z(p;
This set of equations can be expressed compactly using Kronecker products [20, 22], by multiplying
the 1 \Theta M 2 input vector x t times I
M\Omega G to obtain z
M\Omega G), or equivalently
The matrix-vector product Eqn. (6) is an alternative view of Eqn. (4) that allows us to derive
an algorithm of dimension two (i.e. an algorithm having index space of dimension two) for row
processing, rather than of dimension three as implied by the matrix-matrix product form of Eqn.
(4). Reformulating algorithms into lower dimensional index spaces, as we do here, is an effective
technique for simplifying the mapping of algorithms of high dimension into simple and efficient
mono-dimensional (linear) array structures, by avoiding the complexities of multiprojection [23].
The intermediate result column vector z in Eqn. (6) can be viewed either as being composed
of M consecutive vectors z(i corresponding to the
rows of matrix Z, or as M vectors interspersed by a stride factor of M [21], corresponding to the
columns of matrix Z, and given by 1g. By adopting
the latter view, column 1-D transform Eqn. (5) can be expressed as
I M )z; (7)
where vector y is defined similarly as x and z. 1
3 Array Synthesis Framework for Separable Transforms
In this section we elaborate on the space time mapping for the two components of the general 2-D
separable transform array: one for row, and one for column processing. For generality, we treat 1-D
transformation as a matrix vector product and do not make use of any additional computational
savings that might be possible due to kernel symmetries. Rather than write a three dimensional
Single Assignment Form (SAF) [23] 2 algorithm for row processing as suggested by Eqn. (4) (con-
ventional M \Theta M matrix-matrix product), we work with Eqn. (6) (matrix-vector product). The
nested loop algorithm for Eqn. (6) is only of dimension two, and can be trivially mapped into
a linear array for processing inputs in raster-scan order, due to the way we constructed vectors x
and z.
1 As a check, recall that an alternate definition for a separable transform is that the basis matrix T of size M 2 \Theta M 2
be written as the Kronecker product of the basis M \Theta M matrix G times itself (see Akansu [1]). Substituting Eqn. (6)
into Eqn. (7), we get
t\Omega IM
)x. Using the property
(A\Omega B)(C\Omega
(AC\Omega BD) [1, 20], we have
that
An algorithm in SAF is one in which all variables are assigned a single value throughout the entire execution, so
statements of the form x(i) / x(i)   b are not allowed.
If we define the
G as ~
Eqn. (6) is given below
Algorithm 3.1
Inputs x(i) and ~
Initially z
z k (i) / ~
Outputs z
Algorithm 3.1 is defined over 2-D index space I = f(i; g. Using linear
space-time mapping methods, under space transformation this algorithm can be
transformed into an array with M Processing Elements (PEs), but with efficiency of only 50% due
to the block diagonal sparse structure of ~
G. By applying a transformation on index variable i
defined
cM , we obtain an equivalent algorithm which has a rectangular Constant
Bounded Index Space (CBIS) [27] given by I 0
g. This new
algorithm can be mapped under into a linear array with M , fully utilized (100% efficient)
PEs.
Localization of broadcast variables is a complex and rich problem [23, 28, 29]. For maintaining
focus we do not address here the details required for writing Algorithm 3.1 in localized form where
all variables have a common domain of support, but rather give below as Algorithm 3.2 one possible
localized algorithm suitable for raster scan processing. (Note that the variables ~
~
z(i; are the localized equivalents of x(k), G and z(k) respectively).
Algorithm 3.2
Inputs x(k), G(i;
for
for
~
~
~
Outputs z(b k
Pictorially, the localized Algorithm 3.2 can be represented by its Dependence Graph (DG) [23] as
shown in Fig. 2 (a) for 4, where inputs x(k) appear along the horizontal k-axis and are shown
with the corresponding elements in matrix X. The operations modeled by a DG node are shown in
Fig. 2 (b). In this figure, arrows represent data dependencies which are seen to be local among DG
nodes for variables x and z. Not shown in the figure are the dependence vectors associated with
variable G, which are horizontal and of length M ; note that the coefficients of matrix G appear in
this DG replicated M times as
one term for each DG section.
Using the linear space map the DG of Fig. 2 (a) is projected along the i axis giving
rise to the linear array in Fig. 2 (c). The linear time map (schedule) used is which is
consistent with raster scan processing. Input stream in raster scan is fed to PE 0 , one data point
per schedule time period, and the output stream becomes available after an initial delay of three
time units at also in raster scan. So for example, at is available at
and at at PE 0 again. Note that although outputs
do not become available at a single PE, every time instant only one PE produces a z value, and
consequently outputs can be collected by means of a bus connecting all PE outputs. 3
The local PE control needed to realize the output of partial results z in Algorithm 3.2, (i.e.
recognize the inter-row boundaries at k 2 in the DG of Fig. 2 (a)) can be
implemented simply by attaching a bit-tag (called decode bit in Fig. 2 (d)) to each input data token
X. This bit is set only for the last element of every row of matrix X, and the PEs accumulate
locally the result of the multiply-accumulate operation on z only when the bit of the inputs is not
set. When a PE finds that an input in X(:; arrives that has the decode bit set, it places
the result on an output link rather than on the internal z register, and then clears the internal z
register. A column of transform coefficients G are stored in a bank of M registers, and are fed to
the multiplier in sequence with wraparound.
Next we derive a SAF algorithm for the transformation which can be viewed either as
(i) matrix G t operating on column vectors Z(:; k) to produce column vectors Y (:;
or as (ii) matrix Z operating on row vectors G t (i; :) to produce row vectors Y (i; :). Using the
latter interpretation, for which in turn is equal to
?From this expression follows the nested-loop SAF algorithm given below which is very similar to
Algorithm 3.1 except from the fact that the basic multiply-add operation is of the type "scalar
times vector plus vector." 4
Algorithm 3.3
Alternatively, outputs can also be systolically pipelined along with inputs X so that they become available at a
single PEM \Gamma1 . See Beraldin et al. [15] for a discussion on output pipelining and systolization for the DFT.
4 So-called saxpy operation in [21].
z
x
(c) row array
2-bit cntr
z
x
output
(d) internal PE structure
to column array
input
(a) row DG
(b) DG node

Figure

2: (a) Dependence Graph (DG) for transformation on rows of X, with
node operations; (c) the resulting row array after applying the linear space-time mapping operators
along with the I/O processor-time schedule, and (d) the internal structure
of processing element PE 2 .
Inputs: matrices Z and G
Initially Y
Outputs:
scalar multiplications, or one multiplication of scalar G(1; 0) by vector Z(1; :) resulting to vector
which is then added to vector Y 0 (0; :).
The most important aspect of this formulation is that it allows us to operate directly on the rows
of intermediate matrix Z thereby exploiting the raster scan order in which Z is becoming available
from the row array and, as a result, avoiding any intermediate matrix transposition. Furthermore,
although Eqn. (5) is a matrix-matrix product defined over a 3-D index space, we have constructed
Algorithm 3.3 as a matrix-vector product over a 2-D index space on variables (i; k) by using vector
data objects instead of scalars.
Similarly to Algorithm 3.2, Algorithm 3.3 can be localized by propagating variable Z(k; :) along
the i-axis and accumulating output Y (i; :) along the k-axis. The DG for localized Algorithm 3.3
is shown in Fig. 3 (a), where basic DG node operation consist of M scalar multiplications and M
scalar additions and data dependencies among DG nodes consist of M parallel vectors (depicted
as one thick arc). Fig. 3 (b) shows the operation modeled by a DG node (saxpy).
We map DG shown in Fig. 3 (a) using the same space transformation applied to the DG of
Algorithm 3.2, namely resulting in the column linear array shown in Fig. 3 (c).
Scheduling for this DG is accomplished by means of a function mapping a computation within
DG node (i; k) to a schedule time period
is an index that
picks the lth element in row Z(k; :), and is the same schedule vector used for the row
array. So for example, if the DG node depicted in Fig. 3 (b) is (i; computation
scheduled for execution on PE 1 at
at and the first computation of DG node (i; namely
is scheduled on PE 2 at Note that this timing function has two resolution
levels, one that assigns nodes in the DG to macro time steps, and one that assigns individual
multiply-add operations within a DG node to actual time steps (clock cycles).
The regularity of this scheduling allows a simple implementation with no additional local PE
control. And, very importantly, the width-M data dependencies for Z(k; :) pointing upwards in
the column DG of Fig. 3 (a) do not map to M physical links. Rather, only a single physical link
suffices to propagate the inputs z since the sequence of z values is produced at a rate of one scalar
token per schedule time period by the row array. And in addition, z tokens are consumed by the
vector array at the same rate they are produced.
(a) column DG (b) DG node
Y (i,0)+.
Y (i,1)+.
Y (i,2)+.
Y (i,3)+.
2-bit cntr.
clock/4
y
z
output
(d) PE structure
(c) column array
z
from row
array
dependencies
map to single physical link
FIFO
FIFO

Figure

3: (a) Dependence Graph (DG) for transformation DG on columns of Z, with
nodes represent scalar-vector multiply and vector-vector add; (c) the resulting column array
after applying the linear space-time mapping operators along with the I/O
processor-time schedule, and (d) the internal structure of processing element PE 2 .
Under array 5 has identical structure as the row array in Fig. 2 (c),
and is shown in Fig. 3 (c) along with its I/O schedule. The difference is in the amount of memory
required by the column array for holding entire rows of Z; Y rather than single elements. In
Fig. 3 (d) we show the internal PE structure, where we see that I/O registers for x and z in the row
array have been replaced by length-M FIFO queues for Z and Y in the column array, respectively.
3.1 Performance Characterization
We now calculate the latency and pipelining period [23] of the proposed array. Fig. 4 shows the
complete I/O schedule for both arrays (row and column) on a transform of size 4, with two
consecutive input matrices X 0 and X 1 . In this figure we see that both arrays are perfectly space
and time matched: output stream z generated by the row array can be immediately consumed by
the column array at the exact same rate of production. As a consequence no intermediate buffering
is required.
The problem latency (total time) for a single separable transform with input X 0 is calculated
as follows: Assuming that the row array starts computation at we know that under schedule
the time at which the first intermediate result is produced is i.e. at
4. From that point on, one intermediate result becomes available per schedule time
period. Hence, the column array requires 28 time steps for processing the entire stream of z tokens
(time between accepting z(0) and producing Y (3; 3)) , or in general MT4
time units. Since the first input to the column array becomes available at , the latency (i.e.
the total time required for the parallel architecture to compute the 2-D separable transform for a
single image
When considering a continuous stream of problems with input matrices fX
relevant measure of performance is the block pipelining period i.e. the number of time steps between
initiations of two consecutive problems [23]. Referring again to Fig. 4, note first that the last data
token of the first input to enter the row array, X 0 does so at
at in the figure), and thus the time at which the first element of the next input X 1 (0; may
enter the row array is t Next, the time at which the column array may start
accepting inputs from intermediate result Z 1 (of problem instance
Assuming that the row array does accept X 1 (0; 0) at the earliest possible time
then units later at the first data token of the second input stream,
becomes available to the column array, which is exactly t column the time at which
the column array is ready to accept it. Consequently, the resulting block pipelining period is
Loosely speaking we will name this second array column array, or vector array to remind us that it produces
the same result Y as the column array implementing the conventional row-column approach in Figure 1 (a). Note
however that in our array, Z is provided in raster scan order i.e. as a sequence of rows, not columns.
. X(0,3) X(1,0) . X(1,3) X(2,0) . X(2,3) X(3,0) . X(3,3)
. Y(0,3) Y(1,0) . Y(1,3) Y(2,0) . Y(2,3) Y(3,0) . Y(3,3)
X0 inputs to row array
. X(0,3) X(1,0) . X(1,3) X(2,0) . X(2,3) X(3,0) . X(3,3)
X1 inputs to row array
Intermediate results Z0
Intermediate results Z1
Column array outputs Y0
time
last X0 input at 1first intermediate result of
Z1 to column array at
Figure

4: Overall time schedule for both arrays considering two successive inputs X 0 and X 1 , with
as in the example). Since there is a single input line to the row
array and thus only one token can be fed per time step the theoretically minimum fi, or maximum
throughput, is attained.
Efficiency is given by the ratio of the total number of computations to the product of latency
times number of processors. Single input efficiency is 50%, and for a fully pipelined array is exactly
100%, implying that one output Y (i; becomes available every schedule time period.
Since the performance metrics depend on the particular schedule chosen here, there are a number
of available tradeoffs that can be explored. For instance, the non-systolic schedule
used for the column array (with the same timing function) leads to a reduced latency of M 2
while the block pipelining period remains minimal, . However, in this case M outputs (a
column of Y ) are becoming available simultaneously at every time step and should be collected in
parallel (Single Input Parallel Output, or SIPO, architecture).
4 Arrays for the 2-D Discrete Fourier Transform
In this Section we derive a parallel architecture for the 2-D Discrete Fourier Transform (DFT) which
is a special case of the computational structures derived in Section 3. Therefore the arrays presented
here are similar to those shown in Figs. 2 and 3, except that by using a DFT algorithm known as
Goertzel's DFT [15] we may reduce kernel coefficient storage and simplify control complexity.
4.1 Goertzel's DFT
The DFT of a 1-D, M-sample sequence x(k) is given by the set of M coefficients
and the DFT of a 2-D, M \Theta M-sample sequence x(i; k) is given by
Clearly, as kernel g lm (i;
li
M the DFT is separable on 1-D kernel g l li
M .
The 1-D DFT array by Beraldin et al. [15] based on the so-called Goertzel Algorithm is derived
by applying Horner's polynomial multiplication rule. The following recursive equation is equivalent
to the 1-D DFT
where z(m; k) denotes the m th DFT coefficient after the k th recursive step, and 0 - m;
Consider, for example, the computation of DFT coefficient z(2) with which, according to
Eqn. (9) is given by
4 x(3). Now, by unraveling
the recursion in Eqn. (11) we have that z(2;
or just z(2;
4 x(0), which indeed is equal to the DFT
coefficient z(2), since W \Gamma6
4 , and W \Gamma8
1. With Goertzel's algorithm, a DFT coefficient
is computed using a single twiddle factor W \Gammam
M , rather than M twiddle factors as with the direct
application of Eqn. (9).
4.2 Modular DFT Array Derivation
Let x and z be the input and intermediate result M 2 \Theta 1 vectors representing the corresponding
raster scan matrices X and Z, as in Section 2. Examination of Eqn. (11) reveals that, together with
the above definitions for x and z, an algorithm very similar to the localized SAF Algorithm 3.2 (see
Section can be formulated, where W \Gammai
M plays the role of G(i; k), and in the statement for updating
z addition takes place before multiplication. Hence, a localized SAF algorithm for Eqn. (11) is given
by
Algorithm 4.1
Inputs x(k), W \Gammai
for
6 In this equation, m is the "frequency" index, and k the "time" index.
~
~
Outputs z(b k
where ~ x(i; k); ~ g(i;
z(i; are the localized equivalents of variables x(k), W \Gammai
M and z(k) re-
spectively. The DG for this algorithm is similar in structure to the one shown in Fig. 2 (a), except
that instead of M 2 G coefficients propagated by length-M horizontal dependencies along the k-
direction, we have only M complex twiddle factors propagated by length-one dependencies along
the k-direction. Also, the order of internal PE operations is reversed.
Algorithm 4.1 is mapped under space transformation
resulting in the row array shown in Fig. 5 (a). As shown in Fig. 5 (b), the internal PE structure
of the DFT row array is simpler than the array for general transforms in that a single complex
twiddle factor is always used as a multiplicand. This sort of simplification can be further exploited
for designing a specialized multiplier for that particular twiddle factor. For example, depending
on the value of W \Gammai
M , or using approximations to W \Gammai
M which are factors of powers of two, a fast
shift-add accumulator can be used rather than an array multiplier.
The DG for the DFT column array, which is the vector version of Eqn. (11), is constructed
similarly as the one for the general case shown in Fig. 3 (a), with DG node shown in Fig. 3 (b).
The vector DG for Goertzel's DFT column algorithm is similar to the general one in Fig. 3 (a),
with the exception that there are only M twiddle factors propagating in the k-direction, and the
operations in the PEs are exchanged.
In order to illustrate a different solution having better latency performance than the arrays of
Section 3, we schedule the column DFT DG with function T 0
mapping
computation l of DG node (i; k) at time period
instead of
Mapping the vector DG under (S; T 0
results in the column array shown in Fig. 5 (c),
with internal PE structure shown in Fig. 5 (d).
In this column array, by using schedule T 0
the need for input FIFO queues holding intermediate
results Z is eliminated, resulting in total latency and memory of nearly one half that of the general
version in Section 3. Schedule T 0
is also applicable to the general arrays of Section 3 and is specific
neither to the DFT nor to Goertzel's formulation. However, to realize the performance gains M
parallel output ports are needed in order to collect in parallel one full column of Y results every
clock cycle.
y
z
output
(d)
(c) column array
z
x
output
(b)
(a) row array
FIFO

Figure

5: Parallel architecture for 2-D DFT with 4. (a) The row array and processor-time I/O
schedule, (b)the internal structure of PE 2 . (c) The column array and processor-time I/O schedule,
(d) the internal PE structure.
4.3 Performance Characterization
The latency of the overall architecture under schedules for the row and T 0
the column arrays is M 2 +M , which is faster than the that of the general architecture developed
in Section 3 (latency of 2M 2 ) by almost a factor of two, due to the broadcasting of intermediate
results z(i) to all PEs in the column array and the availability of M parallel output ports. As a
result, entire columns of the final result Y become available simultaneously, as shown in Fig.5 (c)
at times pipelining period is identical to that of the array in Section 3,
i.e. . However, the smaller latency obtained with broadcasting comes at the expense of a
potentially larger duration for the each time step, i.e. a reduced maximum rate at which we can
run the clock, as M increases.
5 Arrays for 2-D separable Convolution
5.1 2-D Convolution
Linear 2-D convolution on M i \Theta M k input x(i; is given by
1. Notice that, in contrast to the general transform
case of Eqn. (1), the summation limits in convolution range over a region that depends on the
support of both functions g(\Delta; \Delta) and x(\Delta; \Delta). Since we deal here with the case where the kernel g
is of smaller support than the input x, the summation limits are defined over the kernel's support
region. 7
Separable filters have the property that g(i
There is a substantial difference in complexity between non-separable and separable filters, the
former requiring L i L k multiplications and L i additions per output sample, while the latter
requiring multiplications and additions per sample.
Eqn. (13) implies that input x(i; can be first processed in raster-scan order by performing
convolutions on its rows, resulting in intermediate result z(i;
The final result is then obtained by performing M k , 1-D column-wise convolutions on z(i; k), as
k). For simplicity of exposition, but without loss of generality, from now
on we assume
7 The method is not restricted to small convolution kernels, but we focus here on VLSI arrays with a relatively
small number of PEs.

Figure

6 (a) shows an example of the DG for row-wise convolution on raster-scan input x(i; k),
3. In this figure we have omitted the data dependence vector arrowheads
indicating the data flow direction for broadcast variables x and g, and have included only those
pointing upwards for the localized accumulation variables z. In contrast to the general transform
case in Section 3 (see Fig. 2), in convolution there is a slight "spill-over" of outputs at row
interface boundaries. In our proposed solution, we let input rows be separated by
slots as they enter the array, and, as a result, there is no need to introduce any control mechanism
for row interface management. Typically, the inefficiency incurred by leaving empty computations
between rows in the schedule is small, in the order of 1%. 8
Let the linear space-time mapping operators
tA , be given by
and indexing processors and t is denoting the schedule time periods.
The resulting array, along with its I/O schedule is shown in Fig. 6 (b), where every
input token is made simultaneously available to all PEs. The arrows in the I/O schedule indicate
the progression of partial results, that lead to the computation of output z(0; 2). No row interface
control is needed as long as are padded at the end of each input row. For example, at
time instants (following input x(0; 5)), two zeros are inserted into the input stream
clearing the filter for the subsequent insertion of row x(1; k),
Using the notation introduced in Section 2, we can write
as
which is convolution based on saxpy row operations. At the vector level, the DG in Fig. 7 (a)
accepts row inputs
produces row outputs Y (i; :) (filter coefficients g(i 0
are still
scalar values). The internal DG node structure shown in Fig. 7 (b) has M k
each one accepting as arguments an element from row
:) and a filter coefficient g(i 0
subsequently adding the product to one y partial result that is propagating from bottom to top,
leading to the computation of the corresponding element of Y (i; :) In Fig. 7 (c) we show the
complete 2-D convolution pipelined architecture consisting of two linear array stages: the first,
shown on the left and called the row array, for processing the DG of Fig. 6 (a), followed by another
linear array stage, called the column array, for processing the vector DG of Fig. 7 (a). Each PE is
8 An alternative to spacing rows by slots is to introduce control tags carried by the last and first elements
of every row, resulting in perfect row pipelining. However, the efficiency improvement of perfect row pipelining over
the solution presented here is negligible.
row 0
row 1 row 2
(a)
(b)

Figure

(a) Row-wise 1-D convolution on x(i; resulting in z(i; k), with M
3. Note that rows are spaced by slots at the input. (b) array for
row-wise convolution and processor-time snapshots of array operation.
y(i,0)+. y(i,1)+. y(i,7)+.
(a) (b)
FIFO
y
(d)
FIFO

Figure

7: (a) Column 1-D convolution DG in vector form. Dark nodes represent
scalar computations and thick edges represent data dependencies among vector variables with
(b) The internal structure of DG node (i; i 0
(c) The
parallel architecture for computing 2-D separable convolutions: the first stage linear array performs
the row and the second the column processing respectively.(d) The internal row and column array
PE structure.
identified by two indices (p; q), where q indicates whether a PE is for row processing or for
column processing
The space transformation for the column array is and the timing function mapping a
computation to a schedule time period is given by
The choice of this timing function implies that PE i 0 ;1 can serialize the execution of the M k
scalar operations modeled by DG node (i; i 0
) in the following
2)g. The two arrays are perfectly space and time matched because
this order corresponds exactly to the order in which outputs
are produced from PE 0;0
of the row array (see Fig. 6). As a result, each intermediate output
is broadcasted to
the column array PEs via a bus as soon as it is produced and there is no need for introducing any
additional control or memory to interface the two pipelined array stages.

Figure

7 (d) shows the internal structure of the PEs of both arrays. During every schedule time
period, PE p;q performs the following three operations: (a) multiplication of the value on the input
bus by the local filter coefficient; (b) addition of the product to the incoming partial result from
(c) propagation of the result to PE p\Gamma1;q . With regards to partial result accumulation
in the column array, since every DG node produces outputs this is also the number of
registers in the FIFO queues required between every pair of PEs (see Fig. 7 (c)). Note, however,
that buses for broadcasting as well as links for moving partial results y(i; into and out of the
FIFOs carry scalar quantities, so their bit-width should be only one word (and does not depend on
M k or L).
5.2 Performance Characterization: Latency and Memory
Total computation time for a single input problem instance x(i; k) is as follows. The row array
computation time is T 9 The column array computation time is
Since the communication between both structures is perfectly pipelined, there is only a delay of one
time unit in producing the first output y(0; 0). As a result, total computation time is T vector + 1.
As an example, consider typical numbers for image processing, where M
which result in 96.6% efficiency (loss due to imperfect row pipelining), while with
efficiency is 98.3%. 11 This relatively small efficiency loss
is the result of having arrays with virtually no control circuitry. Alternatively, this efficiency loss
rows, each requiring processing time of Mk +Lk \Gamma 1. The insertion of zeros between rows is the cause of the
steps.
the total number of computations is and since there are 2L processors the
efficiency is given by
could be reduced by introducing control at the PE-level for managing the interfaces between rows
in the row array (as opposed to insertion of zeros).
The block pipelining period on an input sequence fX
implying that a complete output y(i; k) is available every fi time steps, substantially improving
efficiency. Efficiency for the case when M
1024 and
Next we prove a Proposition on the lower bound on memory of any architecture that processes
inputs under raster scan order, and show that the 2-D convolution array proposed here nearly
attains this lower bound. This Proposition shows that this array does not require larger memory
storage than the conventional transposition memory solutions, a result holding also for the
structures presented in Sections 3 and 4.
In the following I/O model, we draw a boundary around the circuit under consideration, and
assume that the input stream x(i; k) crosses this boundary in raster scan, or row-major order. 12
Once an input has entered the circuit, it is held in local storage until used in some computation for
the last time, and then it is discarded. It is further assumed that once input x(i; k) has appeared
in the input stream, the circuit cannot request x(i; k) from the host again, and has to manage it
throughout its lifetime. (This analysis is similar to that in [30].)
Proposition 5.1 The lower bound on memory for 2-D convolution by separable L-tap and by non-separable
filters with row-scan input is
Proof: Consider output y(i; k). By Eqn. (12) (non-separable) or Eqn. (13) (separable) we know
that computation of y(i; requires the set of L 2 input elements L+1)g.
Under the row-scan order imposed by Def. 2.1, the last element in this set to enter the circuit
boundary is x(i; k). The minimum time x(i; spends inside the boundary can be calculated by
looking for the last output element that requires x(i; its computation, which can be shown
to be y(i 1). Making the assumption that the circuit is operating under perfect
pipelining, and hence that one output element is available every single time period, the minimum
time needs to spend inside the boundary is
are rows, each with Finally,
during time t, at least (M k new input elements entered the boundary, and therefore
minimum memory is lower bounded by t storage elements. 2
Total memory requirements for the proposed convolution array of Fig. 7 (c) are dominated by
the column FIFOs, accounting for a total of (M k memory registers, and meeting
exactly the lower memory bound in Proposition 5.1. In addition, there is a total of 2L registers
for locally storing filter coefficients (which could be reduced to L by having corresponding PEs in
12 The derivation also holds if we substitute column-major for row-major order, as long as we are consistent.
each array share their coefficient), and L registers for the accumulation of z(i; results in
the row array. The array proposed in Section 4 also meets this lower memory bound. In the case
of a global transform (i.e. DFT) with equal support regions for kernel and input, Proposition 5.1
reduces to the trivial case where storage for an entire M \Theta M frame is required.
6 Conclusions
In this paper we have shown that the basis for constructing modular architectures for 2-D separable
transforms is an array for the 1-D transform, and that the resulting overall architecture can reflect
the modularity of that array. The only restriction required by the 1-D transform array is that it
accepts input and generates output in raster scan order. Then a regular architecture consisting of
two fully pipelined linear array stages with no shared memory can be derived, based on the key idea
that the DG for the first linear array (performing the "row" processing) is structurally equivalent
to the DG of the second linear array (for the "column" processing), but with the difference that the
latter DG operates at a coarser level of granularity (vector instead of scalar computations). This
idea can be exploited even for separable transforms having a complex computation structure, since,
in principle, avoiding transposition and shared memory does not depend on the 1-D transform
structural details. Rather, it depends only on (i) there being a raster scan order for the 1-D
transform and (ii) a change in abstraction level from scalar to vector operations. The 2-D DFT and
the 2-D convolution have been used as case studies to demonstrate the generality and effectiveness
of the proposed architecture. We further conjecture that the same method can be used to derive
pipelined array structures for multi-dimensional separable transforms of any dimension.
We have also applied this methodology to the Discrete Wavelet Transform, a time-scale transform
that represents finite energy signals with a set of basis functions that are shifts and translates
of a single prototype mother wavelet. Modular arrays for computing the 1-D DWT have been proposed
in the literature [24, 25, 26], but 2-D solutions usually rely on transposition, or some form of
complex data routing mechanism that destroys modularity. Since most of the current applications
and research on wavelets is based on separable FIR wavelet filters, we have been able to apply the
results from Section 5 to the 1-D arrays presented in [24, 25].
All Dependence Graphs and space-time transformation pairs derived here, can be supplied as
input to the rapid prototyping tool DG2VHDL developed in our group [31, 32] that can generate
automatically high quality Hardware Description Language models for the individual PEs and the
overall processor array architectures. The generated VHDL models can be used to perform behavioral
simulation of the array and then synthesized (using industrial strength behavioral compilers,
such as Synopsys [35]) to rapidly generate hardware implementations of the algorithm targeting a
variety of devices, including the increasingly popular Field Programmable Gate Arrays
[33, 34], or the Application Specific Integrated Circuits (ASICS). Using the ideas presented here,
modular linear array cores for the 2-D Discrete Cosine Transform (DCT) which supporting a variety
of different I/O patterns have been synthesized and compared in [36].



--R

"Multiresolution Signal Decomposition"
"Discrete Cosine Transform: Algorithms, Advantages, and Applica- tions"
"An Novel CORDIC-based Array Architecture for the Multidimensional Discrete Hartley Transform"
"A Cost-Effective Architecture for 8x8 Two-Dimensional DCT/IDCT Using Direct Method"
"Image and Video Compression Standars"
"VLSI Architectures for Multidimensional Fourier Transform Processing"
"VLSI Architectures for Multidimensional Transforms"
"Low-Power Video Encoder/Decoder Using Wavelet/TSVQ With Conditional Replenishment"
"Multidimensional Systolic Arrays for Multidimensional DFTs"
"Highly Parallel VLSI Architectures for the 2-D DCT and IDCT Computations"
"Synchronizing Large Systolic Arrays"
"Synthesizing Linear Array Algorithms from Nested For Loop Algo- rithms"
"VLSI Signal Processing Systems"
"Fourier Transforms in VLSI"
"Efficient One-Dimensional Systolic Array Realization of the Discrete Fourier Transform"
"Multidimensional Digital Signal Processing"
"Design of 2-D Recursive Filters with Separable Denominator Transfer Functions"
"The Design of Multistage Separable Planar Filters"
"Parallel Algorithms and Architectures for Discrete Wavelet Transforms"
"Fundamentals of Digital Image Processing"
"Matrix Computations"
"Algorithms for Discrete Fourier Transform and Convolution"
VLSI Array Processors.
"On the Synthesis of Regular VLSI Architectures for the 1-D Discrete Wavelet Transform"
"Distributed Memory and Control VLSI Architectures for the 1-D Discrete Wavelet Transform"
"1-D Discrete Wavelet Transform: Data Dependence Analysis and Synthesis of Distributed Memory and Control Array Architectures"
"On Time Mapping of Uniform Dependence Algorithms into Lower Dimensional Processor Arrays"
"Systolic Array Implemenation of Nested Loop Pro- grams"
"Synthesizing Systolic Arrays with Control Signals from Recurrence Equa- tions"
"Computational Aspects of VLSI"
"DG2VHDL: A tool to facilitate the high level synthesis of parallel processing array architectures"
"DG2VHDL: a tool to facilitate the synthesis of parallel VLSI architectures"
"Using DG2VHDL to Synthesize an FPGA Implementation of the 1-D Discrete Wavelet Transform"
Synthesis of Array Architectures of Block Matching Motion Estimation: Design Exploration using the tool DG2VHDL.
Behavioral Synthesis.
Design and Synthesis of maximum throughput parallel array architectures for real-time image transforms
--TR
VLSI array processors
VLSI Architectures for multidimensional fourier transform processing
Synthesizing Linear Array Algorithms from Nested FOR Loop Algorithms
Fundamentals of digital image processing
Discrete cosine transform: algorithms, advantages, applications
VLSI Architectures for Multidimensional Transforms
Behavioral synthesis
VLSI Signal Processing Systems
Image and Video Compression Standards
Multiresolution Signal Decomposition
Multidimensional Digital Signal Processing
On Time Mapping of Uniform Dependence Algorithms into Lower Dimensional Processor Arrays
Synchronizing large VLSI processor arrays
Parallel architectures and algorithms for discrete wavelet transforms
VLSI architectures;parallel processing;2-D separable transforms