Chrestenson transform FPGA embedded factorizations

Chrestenson generalized Walsh transform factorizations for parallel processing imbedded implementations on field programmable gate arrays are presented. This general base transform, sometimes referred to as the Discrete Chrestenson transform, has received special attention in recent years. In fact, the Discrete Fourier transform and Walsh–Hadamard transform are but special cases of the Chrestenson generalized Walsh transform. Rotations of a base-p hypercube, where p is an arbitrary integer, are shown to produce dynamic contention-free memory allocation, in processor architecture. The approach is illustrated by factorizations involving the processing of matrices of the transform which are function of four variables. Parallel operations are implemented matrix multiplications. Each matrix, of dimension N × N, where N = pn, n integer, has a structure that depends on a variable parameter k that denotes the iteration number in the factorization process. The level of parallelism, in the form of M = pm processors can be chosen arbitrarily by varying m between zero to its maximum value of n − 1. The result is an equation describing the generalised parallelism factorization as a function of the four variables n, p, k and m. Applications of the approach are shown in relation to configuring field programmable gate arrays for digital signal processing applications.

requirements, of shuffle operations and of the number of required memory partitions (Corinthios 1994). The factorizations are developed with a view to implementation as embedded architectures of presently available FPGAs (Harmut et al. 2010;Huda et al. 2014).
The algorithms and corresponding architectures relate to general base matrix factorizations (Corinthios 2009). Rotations of a base-p hypercube, where p is an arbitrary integer, produce dynamic memory allocation, in processor architecture. The approach produces factorizations involving the processing of matrices of the transform which are function of four variables. Parallel operations are implemented matrix multiplications. Each matrix, of dimension N × N, where N = p n , n integer, has a structure that depends on a variable parameter k. The level of parallelism, in the form of M = p m processors can be chosen arbitrarily by varying m between zero to its maximum value of n − 1. The result is an equation describing the generalised parallelism factorization as a function of the four variables n, p, k and m. Applications of the approach are shown in relation to configuring field programmable gate arrays for digital signal processing applications.
Hypercube transformations have been applied to diversified problems of information processing. The present paper describes an approach for FPGA parallel processor configuration using an arbitrary number M of general-base processing elements, where M = p m , p being the general radix (base) of factorization. The input data vector dimension N, or input data matrix dimension N × N, where N = p n , the radix, or base, p of factorization of the transformation matrix, the number of processors M, and the span of the matrix, that is, the spacing between data simultaneously accessed are all variable. A unique optimal solution yielding a progressive degree parallel to massively parallel optimal architectures is presented.

Matrix structures
In what follows some definitions relating to the special structure of sparse, permutation and transformation matrices (Corinthios 1994) are employed. In particular matrix span is taken to mean the distance between two successive nonzero elements along a row or a column. A fixed topology processor is one that accesses data in a fixed geometry pattern where data points are equidistant throughout the different iterations, thus requiring no addressing. A shuffle-free algorithm is one that necessitates no data shuffling between iterations. A p k -optimal algorithm is one that requires access of matrix elements which are spaced by a minimum distance of N/p k elements. In addition we adopt the following definitions.

General base processing element
In what follows a general-base processing element PE with a base, or radix, p is a processor that receives simultaneously p input operands and produces simultaneously p output operands. The PE in general applies arithmetic or weighting operations on the input vector to produce the output vector. In matrix multiplication operations for example the PE applies a p × p matrix to the p-element input vector to produce the p-element output vector. The matrix elements may be real or complex.
Due to the diversified general applicability of such a processing element a universal processing element (UPE), which can be constructed in a 3D-type architecture has been proposed (Corinthios 1985). In the present context a UPE may be seen simply as a general base-p processing element PE as defined above, accepting p inputs, weighting them by the appropriate p × p matrix and producing p output operands.

Pilot elements, pilots matrix
Similarly to signals and images an N × N matrix may be sampled and the result is "impulses", that is, isolated elements in the resulting N × N samples matrix. We shall assume uniform sampling of rows and columns yielding p uniformly spaced samples from each of p rows and element alignment along columns, that is, p uniformly spaced samples along columns as well as rows. The samples matrix which we may refer to as a "frame" thus contains p rows of p equally spaced elements each, a rectangular grid of p 2 impulses, which we may refer to as "poles", which we shall call a "dispatch". With N = p n the N 2 elements of the "main" (or "parent") matrix, that is, the original matrix before sampling, may be thus decomposed into N 2 /p 2 = p n−2 such dispatches.
By fixing the row sampling period as well as the column sampling period, the row and column spans of the resulting matrix are known. It therefore suffices to know the coordinates (indices) of the top left element, that is, the element with the smallest of indices, of a dispatch to directly deduce the positions of all its other poles. The top left element acts thus as a reference point, and we shall call it the "pilot element". The other p 2 − 1 elements associated with it may be called its "satellites".
In other words if the element a ij of A is a pilot element, the dispatch consists of the elements c and r being the column and row element spacing (spans), respectively.
A processing element assigned to a pilot element can thus access all p 2 operands of the dispatch, having deduced their positions knowing the given row and column spans.
Since each pilot element of a frame originated from the same position in the parent matrix we can construct a "pilots matrix" by keeping only the pilot elements and forcing to zero all other elements of the parent matrix. The problem then is one of assignment, simultaneous and/or sequential, of the M = p m processors to the different elements of the pilots matrix.

Hypercube dimension reduction
The extraction of a pilots matrix from its parent matrix leads to a dimension reduction of the hypercube representing its elements. The dimension reduction is in the form of a suppression, that is, a forcing to zero, of one of the hypercube digits. Let C = (j n−1 , …, j 1 j 0 ), j k ∊ {0, 1, 2, …, p − 1} be an n-digit base-p hypercube. We will write Ck to designate the hypercube C with the digit k suppressed, that is, forced to zero. Several digits can be similarly suppressed. For example, C 2,4 = j n−1 . . . j 5 0j 3 0j 1 j 0 , and C n−1 = 0j n−2 . . . j 1 j 0 . a i+kc,j+lr ; k = 0, 1, . . . , p − 1, l = 0, 1, . . . , p − 1

Parallel configuration algorithm
A sequence of perfect shuffle operations effected through simple hypercube transformations can be made to broadcast the parallel configuration and access assignments to the different processors. The overall approach is described by the following algorithm.
The parallel dispatch, state assignment and sequencing Algorithm 1 dispatches the M = p m processors for each stage of the matrix factorization. The base-p m-tuple (i m−1 i m−2 …i 1 i 0 ) is assigned to the parallel processors. The (n − m) tuple (j n−1 j n−2 …j m ) is assigned to the sequencing cycles of each processor. The algorithm subsequently applies hypercube transformations as dictated by the type of matrix, the stage of matrix factorization and the number of dispatched processors. It tests optimality to determine the type of scan of matrix elements to be applied and evaluates parameters such as pitch and memory optimal queue length, to be defined subsequently, it accesses the pilot elements and their satellites, proceeding to the parallel dispatch and sequencing of the processing elements.
Each processing element at each step of the algorithm thus accesses from memory its p input operands and writes into memory those of its output operands. The algorithm, while providing an arbitrary, generalised, level of parallelism up to the ultimate massive parallelism, produces optimal multiprocessing machine architecture minimizing addressing, the number of memory partitions as well as the number of required shuffles. Meanwhile it produces virtually wired-in pipelined architecture and properly ordered output.

General matrix decomposition
In developing techniques for the general-base factorization of transformation matrix multiplications it is convenient to effect a decomposition of a matrix into the sum of matrices. To this end let us define an "impulse matrix" as the matrix δ(i, j) of which all the elements are zero except for the element at position (i, j), that is, can be written as the sum where the δ(i, j) matrices are of dimension N × N each. The matrix A can thus be written in the form Furthermore, in the parallel processing of matrix multiplication to a general base p it is convenient to decompose an N × N matrix with N = p n as the sum of dispatches, a dispatch being, as mentioned earlier, a matrix of p 2 elements arranged in a generally rectangular p × p pattern of p columns and p rows. Denoting by σ R and σ C the row and columns spans of a dispatch we can decompose a matrix A into the form More generally we may wish to decompose A in an order different from the uniform row and column scanning as in this last equation. In other words we may wish to pick the dispatches at an arbitrary order rather than in sequence. As mentioned above, we shall call the top left element the pilot element and its p 2 − 1 companions its satellites. In this last equation the pilot elements are those where k = 1 = 0.
To effect a parallel matrix decomposition to a general base-p we use hypercubes described by base-p digits. The order of accessing the different dispatches is made in relation to a main clock. The clock K is represented by the hypercube to base p as Its value at any time is given by At each clock value K a set of M UPE's (PE's) is assigned a set of M dispatches simultaneously. We will reserve the symbols w and z to designate the row and column indices of a pilot element at clock K. In other words, at clock K each selected pilot element shall be designated a w,z , that is, [A] w,z where w and z are functions of K to be defined. They will be determined in a way that optimizes the parallel and sequential operations for the given matrix structure and the number M = p m of available UPE's.
With M = p m base-p processing elements the hypercube representing K shall be rewritten in the form where we have written where the "parentheses" 〈 and 〉 enclose the elements accessed in parallel. In what follows we write P ν,μ to designate the pilot element of processor no. ν at real time clock μ.

Application to the CGW transforms
The lowest order base-p Chrestenson generalised Walsh CGW "core matrix" is the p-point the Discrete Fourier matrix where In the following, for simplicity, the scaling factor 1/ √ p will be dropped. We start by deriving three basic forms of the Chrestenson (generalised Walsh GW) transform in its three different orderings: in natural order CGWN, in Walsh-Paley order CGWP and in Walsh-Kaczmarz order CGWK.

The CGWN transformation matrix
The CGWN transformation matrix W N for N = p n data points is obtained from the generalised-Walsh core matrix W p by the Kroneker multiplication of W p by itself n times.

CGWP transformation matrix
The generalised Walsh transform in the CGWP order is related to the transform in natural order by a digit-reverse ordering. The general-base digit reverse ordering matrix K N (p) can be factored using the general-base perfect shuffle permutation matrix P (p) , also denoted simply P, and Kroneker products where I K is the identity matrix of dimension K.
Operating on a column vector x of dimension K the base-p perfect shuffle permutation matrix of dimension K × K produces the vector The CGWP matrix W N,WP can thus be written in the form

CGWK transformation matrix
The CGWK transformation matrix is related to the CGWP matrix through a p-ary to Gray transformation matrix G N (p) .
The following factorizations lead to shuffle-free optimal parallel-pipelined processors.

CGWN optimal factorization
A fixed topology factorization of the CGWN transformation matrix has the form which can be re-written in the form  and F = CP, noting that the matrix F is p 2 -optimal.

CGWP optimal factorization
We fixed topology factorization of the CGWP matrix has the form Letting we obtain where each matrix Q i , i = 0, 1, …, n − 2, is p 2 -optimal, while Q n−1 is p-optimal.

CGWK optimal factorization
The fixed topology CGWK factorization has the form Letting A quasidiagonal matrix is a matrix containing matrices along its diagonal and null matrices elsewhere.
where Letting we have with The factorization can also be re-written in the form where The matrices Γ i are p 2 -optimal, except for Γ 0 which is maximal span. These are therefore optimal algorithms which can be implemented by an optimal parallel processor, recirculant or pipelined, with no shuffling cycle called for during any of the n iterations.

Application to image processing
The potential in enhanced speed through high-level parallelism of the optimal algorithms is all the more evident within the context of real-time image processing applications. For 2D signals, algorithms of generalised spectral analysis can be applied on sub-images or on successive column-row vectors of the input image. Factorizations of the algorithms of the Chrestenson transform applied on an N × N points matrix X representing an image, with N = p n can be written for the different transform matrices. The CGWN 2D transformation for optimal pipelined architecture can be written in the form where T stands for transpose. The CGWP factorization can be written in the form The CGWK factorization for optimal pipelined architecture can be written in the form These fast algorithms are all p 2 -optimal requiring no shuffling between iterations of a pipelined processor. In applying these factorizations the successive iterations are effected on successive sub-images such that after log p N stages the transform image Y is pipelined at the processor output. Applications include real-time processing of video signals.
The Discrete Fourier transform matrix for N points is the matrix F N defined above in (10) with p replaced by N and the factor 1/ √ p optional: For images the factorization leads to the optimal form and for unidimensional signals the corresponding form for the Discrete Fourier matrix is (36)

Perfect shuffle hypercube transformations
The hypercube transformations approach is illustrated using the important matrices of the Chrestenson generalised Walsh-Paley (CGWP), generalised Walsh-Kaczmarz (CGWK) and the Discrete Fourier transforms.
We note that the matrices C k in the Discrete Fourier transform expansion are closely related to the matrices J i and H i in the Chrestenson generalised Walsh Paley factorization. In fact the following relations are readily established: where the equality ≜ sign means equal by definition.
Therefore, the CGWP matrices Q i are the same as the C i matrices defined above and have the same structure as the F i matrices in the Fourier matrix factorization. Writing the post-multiplication by H k has the effect of permuting the columns of C so that at row w, the pilot element is at column z as determined by the permutation H k , that is, with the special case k = n − 2 producing and that of k = n − 1 yielding Alternatively, we can write z directly as a function of w by using previously developed expressions of permutation matrices. For example, and using the expression defining P, namely, (44) (53) B 0 = CH 0 = CP with k = 1, we can write a relation that defines the pilot elements matrix.
Similarly, and from the definition given in Corinthios (1994): with i = 1 and t = 1 we have

Consider the permutation matrix
Let the base-p hypercube describing the order in a vector x of N = p n elements be represented as the n-tuple.
The application of the matrix R p n on the n-tuple vector x, results in the n-tuple: We note that with respect to x the left k digits and the right m digits are left unchanged while the remaining digits are rotated using a circular shift of one digit to the right.
The pilot-elements matrix β k corresponding to the matrix B k is obtained by restricting the values of w (and hence the corresponding z values) to w = 0, 1, …, p n−1 − 1.
Moreover, we note that if we write and note that G i is similar in structure to C N , we have (59) R N = R p n = I p m × P p j × I p k .
To obtain the pilot elements matrix λ i corresponding to L i we write in order to reveal all satellite elements accompanying each pilot element. We then eliminate all the repeated entries in z′ and the corresponding w values, retaining only pilot elements positions. Alternatively we simply force to zero the digit of weight n − 2 in w and that of weight n − 1 in z.

The CGWP factorization
We presently focus our attention on the matrices In evaluating the pilot elements coordinates we begin by setting the number of processors M = 1. The corresponding w − z relation of the pilot elements are thus evaluated with m = 0. Once this relation has been established it is subsequently used as the reference "w − z conversion template" to produce the pilot element positions for a general number of M = p m processors. A "right" scan is applied to the matrix in order to produce the w − z template with an ascending order of w. In this scanning type the algorithm advances the first index w from zero selecting pilot elements by evaluating their displacement to the right as the second index z. Once the template has been evaluated the value m corresponding to the number of processors to be dispatched is used to perform successive p-ary divisions in proportion to m to assign the M processors with maximum spacing, leading to maximum possible lengths of memory queues. A "down" scan is subsequently applied, where p-ary divisions are applied successively while proceeding downward along the matrix columns, followed by a selection of the desired optimal scan.
The template evaluation and subsequent p-ary divisions for the assignment of the M processors through a right type scan produce the following hypercube assignments. The assignments are as expected functions of the four variables n, p, k and m. The conditions of validity of the different assignments are denoted by numbers and letters for subsequent referencing. With K denoting the main clock, the following hypercube transformations are obtained

Row and column scans for optimal assignment
A processor is considered optimal if it requires a minimum of memory partitions, is shuffle free, meaning the absence of clock times used uniquely for shuffling, and produces an ordered output given an ordered input. It is shown in Corinthios (1994) that p 2 -optimal algorithms and processors lead to a minimum number of p 2 partitions of N/p 2 queue length each. With M = p m base-p processors operating in parallel the number of partitions increases to p m+2 and the queue length of each partition reduces to N/p m+2 .
An optimal multiprocessing algorithm should satisfy such optimality constraints. The horizontal spacing between simultaneously accessed pilot elements defines the input memory queue length. The vertical spacing defines the output memory queue length. With M processors applied in parallel the horizontal spacing between the accessed elements will be referred to as the "input pitch", while the vertical spacing as the "output pitch".
By choosing the pilot elements leading to the maximum possible pitch, which is the highest of the two values: the minimum input pitch and minimum output pitch, optimality in the form of N/p m+2 queue length is achieved.
We note that optimal minimum memory queue length MMQL satisfies The following algorithm, Algorithm 2, describes this approach to state assignment optimality.
In following the algorithm we note that in the validity condition y of the B k matrix y : 1 ≤ m ≤ n − k − 2 the results obtained are such that the digit i 0 of w is of a weight p k . Hence the input pitch is p k while the output pitch which can be deduced from the position of i 0 in z is p n−1 , that is, maximal possible. The input pitch is thus function of k and can be low if k is small. By performing a down scan of B k we obtain the following solution: where now it is i m−1 that leads to a minimum pitch and it has a weight of p n−m−1 in w and p n−m−2 in z. We deduce that the minimum pitch in this solution is p n−m−2 , which is the optimal sought. The same reasoning leads to the optimal assignment for the case These are the only two cases of the matrix that need be thus modified for optimality. All results obtained above for the other validity conditions can be verified to be optimal.

Matrix span
In the above from one iteration to the next the value of k is incremented. In each iteration once the pilot element matrix coordinates (w, z) are determined as shown above each processor accesses p elements spaced by the row span starting with the pilot element and writes its p outputs at addresses spaced by the column span. The row and column spans of a matrix are evaluated as is shown in Corinthios (1994). In particular we note that the matrix has the same column span as that of C, namely σ c (B k ) = σ c (C) = p n−1 . The row span of B k is evaluated by noticing that B k has the same structure as C with its columns permuted in accordance with the order implied by The transformation of the hypercube (i n−1 …i 1 i 0 ) corresponding to H k −1 is one leading to a most significant digit equal to i n−2 . Since this digit changes value from 0 to 1 in a cycle length of p n−2 we deduce that the row span of all the B k matrices is simply Each processing element thus accesses p operands spaced p n−2 points apart and writes their p outputs at points which are p n−1 points apart.

The CGWK factorization
The sampling matrices of the CGWK factorization are more complex in structure than the other generalised spectral analysis matrices. They are defined by (91) Ŵ i = P −1 G i S i+1 (92) L i P −1 G i we have We note that the sampling matrix G i has the same structure in poles and zeros, that is, in the positions of non-zero and zero elements respectively, as that of the matrix C N . We can write for the matrix G i as the pilot elements positions.
Given the definition of the matrix L i a hypercube rotation corresponding to the matrix P −1 would yield the w and z values of L i as: Alternatively, a z-ordered counterpart can be written as: Similarly, the matrix Γ 0 = G 0 S 1 which is obtained from G 0 by permuting its columns according to the order dictated by leads to the m = 0 template assignment and a similar z-ordered state assignment counter part. For we have which leads to the state template assignment (93) Ŵ i = L i S i+1 .