In this paper we describe in some detail the architectural features of the CEDAR. multiprocessor. We also discuss strategies for implementation of dense matrix computations, and present performance results on one cluster for a variety of linear system solvers, eigenvalue problem solvers, as well as algorithms for solving linear least squares problems.

We propose a novel algorithm and architecture for minimum variance distortionless response (MVDR) beam-forming. Our approach extracts the least squares residual element directly, and requires only one systolic array for the many look directions.

Two different techniques are proposed for implementing an adaptive beamformer using a systolic QR decomposition array based on real rather than complex Givens rotations. These are compared in terms of their computational requirements, sensitivity to phase and amplitude imbalance, and effect on the reception of a desired signal.

In this paper, by using the matrix decomposition method, the Kalman filter can be formulated as a modified SRIF data processing problem followed by a QR operation. Compared with the conventional SRIF method, this approach simplifies the computational structure, and is more reliable when the system has a singular(or near singular) coefficient matrix. By skewing the order of input matrices, fully pipelined systolic2Kalman filtering operation can be achieved. With the number of processing units of the 0(n ), the system throughput rate is of the 0(n). The numerical properties of the systolic Kalman filtering algorithm under finite word length effect are studied via analysis and computer simulations, and are compared with those of conventional approaches.

Extended Kalman filtering provides for non-linearities in the equations of motion by propagating the state per a differential equation. For the predictor part of the Kalman filter, the Runge-Kutta differential equation solver can be used to extrapolate each new state numerically. This paper explores the systolic implementation of the Runge-Kutta algorithm. For the corrector part of the Kalman filter, this paper outlines an ESL realization of a systolic array backsolver.

The Luk QR decomposition (QRD) and singular value decom-position (SVD) systolic architectures are synthesized into one Gentleman-Kung triangular architecture. Two of these arrays are connected point-to-point, forming a three dimensional architecture suitable for matrix transposition. Overlapping the diagonal processors results in a square array for matrix multiplication using the engagement algorithm. An augmented architecture is described that implements all of the above algorithms with increased throughput for the QRD. Fault tolerance methods for the Luk QRD and SVD algorithms implemented on these new architectures will be presented. The fault tolerance methods to be examined will either detect a transient error and recover the correct solution, or locate the processor where the error occurred, allowing for reconfiguration.

A survey is presented of robust methods for linear deconvolution. Iteratively Reweighted Least Squares (IRLS), Residual Steepest Descent (RSD), linear programming, and other techniques are discussed. Applications of robust regression are discussed and numerical examples are presented.

A high resolution algorithm is presented for resolving multiple incoherent and coherent plane waves incident on an array of sensors. The incident sources may be a mixture of narrowband and broadband sources, and, the array's geometry is unrestricted. This algorithm makes use of a fundamental property associated with the signal eigenvectors of the array's spectral density matrix. Specifically, it is shown that these signal eigenvectors may be represented as linear combinations of steering vectors which identify the directions of the incident plane waves. As such, linear algebraic methods may be efficiently used to estimate the plane wave directions-of-arrivals. Simulation results are presented to illustrate the high resolution performance achieved with this new approach relative to that obtained with smoothed MUSIC and the Wang-Kaveh algorithms.

An advanced signal processing laboratory has been assembled at Martin Marietta in Baltimore, Maryland, based on the SAXPY MATRIX 1 supercomputer. This parallel processor has a 1 Gflop peak computational capacity. It will be used to test sonar algorithms based on vector and matrix computations.

An interference cancellation technique that combines a signal subspace approach with adaptive adjustment of the weight vector is derived. The weight vector is constrained to lie in a signal subspace computed by eigendecomposition of the covariance matrix of the array outputs. The weight vector is rotated to maximize the output signal-to-interference ratio. The technique does not require array calibration as in other subspace methods. The proposed technique applies to receiver arrays for digital communication signals. The case of a BPSK signal modulating an AM transmitter is analyzed in some detail. Numerical examples illustrating the behavior of the algorithm are provided.

An adaptive implementation of the MUSIC (Multiple Signal Classification) algorithm for estimating time-varying wavenumber power spectra is described. The algorithm is based on the iterative computation of the SVD (singular value decomposition) of a data matrix which is changing as new data becomes available.

In adaptive least-squares estimation problems, a desired signal d(n) is estimated using a linear combination of L observation values samples xi (n), x2(n), . . . , xL-1(n) and denoted by the vector X(n). The estimate is formed as the inner product of this vector with a corresponding L-dimensional weight vector W. One particular weight vector of interest is Wopt which minimizes the mean-square between d(n) and the estimate. In this context, the term `mean-square difference' is a quadratic measure such as statistical expectation or time average. The specific value of W which achieves the minimum is given by the prod-uct of the inverse data covariance matrix and the cross-correlation between the data vector and the desired signal. The latter is often referred to as the P-vector. For those cases in which time samples of both the desired and data vector signals are available, a variety of adaptive methods have been proposed which will guarantee that an iterative weight vector Wa(n) converges (in some sense) to the op-timal solution. Two which have been extensively studied are the recursive least-squares (RLS) method and the LMS gradient approximation approach. There are several problems of interest in the communication and radar environment in which the optimal least-squares weight set is of interest and in which time samples of the desired signal are not available. Examples can be found in array processing in which only the direction of arrival of the desired signal is known and in single channel filtering where the spectrum of the desired response is known a priori. One approach to these problems which has been suggested is the P-vector algorithm which is an LMS-like approximate gradient method. Although it is easy to derive the mean and variance of the weights which result with this algorithm, there has never been an identification of the corresponding underlying error surface which the procedure searches. The purpose of this paper is to suggest an alternative approach to providing adaptive solutions to problems in which samples of d(n) are unavailable.

We derive a generalized EM algorithm for the maximum-likelihood estimation of the directions-of-arrival of multiple narrowband signals in noise, and prove that all of the limit points of the algorithm are stable and satisfy the neces-sary maximizer conditions. The deterministic signal model is considered here in which estimates of the unknown signal amplitudes are generated.

Time-Frequency analysis based on the Wigner-Ville Distribution (WVD) is shown to be optimal for a class of signals where the variation of instantaneous frequency is the dominant characteristic. Spectral resolution and instantaneous frequency tracking is substantially improved by using a Modified WVD (MWVD) based on an Autoregressive spectral estimator. Enhanced signal-to-noise ratio may be achieved by using 2D windowing in the Time-Frequency domain. The WVD provides a tool for deriving descriptors of signals which highlight their FM characteristics. These descriptors may be used for pattern recognition and data clustering using the methods presented in this paper.

The Wigner Distribution Function (WDF) is a time-frequency descriptor capable of tracking the time-varying second order statistics in a signal. In this paper, we characterize linear systems in terms of the WDFs of the inputs and outputs. These input/output relations are provided for both continuous-time and discrete-time systems. An application of these results for the identification of a random, time-varying system is suggested.

Parallel algorithms for high level signal processing functions common to tracking and data fusion are presented and stbdied. The algorithms are general techniques for clustering, partitioning and ennumeration and are applicable in a variety of situations. Parallel algorithms were implemented on a hypercube parallel computer and performance data is presented.

A key task in Inverse Synthetic Aperture Radar (ISAR) imaging and many other applications is estimating the power spectrum of a two-dimensional random process from data measurements. Often the data sampling points do not correspond to a uniformly-spaced rectangular lattice. A particular method is reported herein for performing spectrum analysis from data measured on an irregular lattice. The method employs certain optimal weights, termed Generalized Prolate Spheroidal Sequences, that are determined from a generalized matrix eigenvec luor problem. Because the computational burden of the eigenvector solution can be impractical for large sampling lattices, computationally efficient sub-optimal approximations to the optimal eigenvector weights are proposed. These approximate weights result from careful modification of both the optimization criterion and the subspace over which the criterion is optimized. Near-optimal results can be obtained with a significant reduction in computation. A numerical example is presented for a particular ISAR application to verify the utility of the approximations.

A new formulation of the problem of high-resolution radar-imaging is given. We consider a diffuse radar-target for which the reflectivity is modeled as a complex-valued Gaussian ran-dom-process with a mean of zero and an unknown covariance function. Delay-doppler measurements containing additive noise are made by illuminating the target with a sequence of pulses of arbitrary shape. An algorithm for estimating the target's scattering function from these measurements is derived using the method of maximum-likelihood estimation and the concept of an incomplete/complete data-model. The reflectivity process is also estimated.

Pisarenko proposed a method for decomposing a random stationary process into a sum of harmonics in white noise. The numerical determination of the frequencies consists of several parts, one of which is the computation of the zeros of a polynomial which is known to vanish on the unit circle only. We describe how this part of the computations can be formulated as an eigenvalue problem for an orthogonal matrix. Several algorithms for such eigenproblems are reviewed, one of which enables highly parallel computations.

We describe here 'one-sided' Jacobi methods for solving eigenproblems on parallel architectures. These methods are modifications of the procedure introduced by Hestenesl, which uses only columns of the matrix. Thus, only 'local' information is needed in each processor. The Jacobi angle is determined from columns of the factored matrix and columns of the matrix which ultimately becomes the eigenvector matrix. The singular value decomposition may be considered as a special case. The eigenvalue decompositions for symmetric, real normal, and, orthogonal matrices may be regarded as generalizations. Extensions to the complex field are easily derived. We also describe a new formation of Jacobi rotation sets. These require only one send and one receive per set, and a complete set is found in the minimal n-1 steps. We use the ordering of the diagonal elements of the matrix so that, in general, one may obtain only a few of the eigenvalues (and vectors) as required. Implementations of these algorithms may be made on various kinds of parallel architectures, including the hypercube.

The recently exploding interest in parallel computing has given rise to many new ways to implement the Jacobi SVD algorithm. We study five new orderings for the method, and show that they are all equivalent to one another. In addition, we establish that these orderings are also equivalent to the classical cyclic-by-rows ordering and thus share the same convergence properties. To summarize, we have shown that the different researchers came up with essentially the same idea independently and concurrently, and we have proved convergence of their new methods, a hitherto open research problem.

Two floating-point radix-2 schemes using on-line arithmetic for implementing the direct two-angle method for SVDs are presented. The first scheme is an on-line variant of the cosine/sine approach and is the fastest of the schemes considered: it performs the 2x2 SVD step in about 2n clock cycles. However, it requires a relatively large number of modules; this number is reduced when some modules are reused, resulting in a time of 3n clock cycles. The number of modules of this on-line version is still larger than that of the conventional one, but this is compensated by the smaller number of bit-slices per module and by the digit-serial communication among modules. The corresponding speed-up ratios are of 5 and 3 with respect to a conventional arithmetic implementation. The second scheme uses an on-line CORDIC approach and performs the 2x2 SVD in about 7n clock cycles and is advantageous because it is more time-area efficient. It results in a speed-up of about 2.5 with respect to the conventional CORDIC implementation.

ESPRIT is an interesting new method for solving the Direction-of-Arrival estimation problem. It involves some rather tricky matrix manipulations. We show how these calculations can be carried out using only unitary transformations of the data. No inverses or cross-products are required making the new method extremely robust.

This paper discusses the ESPRIT direction finding method of Paulraj, Roy, and Kailath [1-5]. It provides a comparison of ESPRIT with the earlier MUSIC algorithm of Schmidt [6-10], and discusses the computational requirements for ESPRIT. The discussion includes approaches to the parallel computation of nonsymmetric eigensystems, and potential application of ESPRIT to spectrum analysis and to improving the resolution of the Wigner-Ville Distribution.

The paper discusses mapping of a Fast Fourier Transform (FFT), Haar Transform and Hadamard Transform algorithms onto a small, two-dimensional, mesh-connected array of processors. The FFT algorithm is an in-place, decimation in frequency, Cooley-Tuckey algorithm in radix 2 and radix 4 versions applied to multidimensional, complex inputs. The data flow of the algorithms has been implemented on the array using an efficient, regular data transfer pattern, uniform for all the algorithms. The inputs and constants used in the algorithms are prestored in the local memories of the processors. The mapping makes it possible to reduce significantly the number of memory locations needed for the constants. A partitioning scheme has been developed for the algorithms which allows us to execute them with inputs of arbitrary size on a small processor array. Also an algorithm has been proposed for the processor array, which efficiently unscrambles the bit reversed output of the FFT algorithm. The processors of the array have East, West, North, South interconnections with their nearest neighbors. The local memory of the processors is small, on the order of hundreds of locations. The processors are controlled in Single Instruction Multiple Data Stream (SIMD) mode and can be selectively disabled using simple masks, consisting of combinations of rows or columns.

Closed form expressions and recursive equations relating the parameters of an ARMA model (which may be non-minimum phase, non-causal or may even contain all-pass factors) with the cumulants of its output, in response to excitation by a non-Gaussian i.i.d. process are derived. Based on these relationships, system identification and order determination algorithms are developed. The output noise may be colored Gaussian or i.i.d. non-Gaussian. When a state-space representation is adopted, the stochastic realization problem reduces to the balanced realization of an appropriate Hankel matrix formed by cumulant statistics. Using a Kronecker product formulation, an exact expression is presented for identifying state-space quantities when output cumulants are provided, or for computing output cumulants when the state-space triple is known. If a transfer function approach is employed, cumulant based recursions are proposed to reduce the AR parameter estimation problem to the solution of a system of linear equations. Closed form expressions and alternative formulations are given to cover the case of non-causal processes.

This paper describes a novel technique using residue codes to detect errors (caused by either permanent or transient faults) in numerical systolic arrays concurrently with the normal operation of the system. A careful analysis of errors is used to drastically reduce the number of residue generators and checkers necessary. Undetectable errors are avoided by suitably choosing the modulo size of the residue code and by slightly modifying the implementation of the multipliers in the truncating circuits or applying few residue code checkers to the array. Error propagation in the array is analyzed in detail to ensure that an erroneous result gen-erated by any adder or multiplier will always be detected at the outputs of the arrays. VLSI implementations of dif-ferent kinds of adders and multipliers are analyzed to show that errors due to faults inside a single bit slice will always produce a detectable error at the output of the arrays. The procedure can be applied to all the processor arrays which can be derived from signal flow graphs.

Given a map in which each position is associated with a traversability cost, the path planning problem is to find a minimum-cost path from a source position to every other position in the map. The paper proposes a dynamic programming algorithm to solve the problem, and analyzes the exact number of operations that the algorithm takes. The algorithm accesses the map in a highly regular way, so it is suitable for parallel implementation. The paper describes two general methods of mapping the dynamic programming algorithm onto the linear systolic array in the Warp machine developed by Carnegie Mellon. Both methods have led to efficient implementations on Warp. It is concluded that a linear systolic array of powerful cells like the one in Warp is effective in implementing the dynamic programming algorithm for solving the path planning problem.

High performance architectures are using an ever increasing number of processors. The Boolean cube network has many independent paths between any pair of processors. It provides both a high communications bandwidth as well as the ability to emulate many other networks without contention for communication channels. Of particular interest for the Fast Fourier Transform (FFT) is the ability to emulate butterfly networks, which defines the communication pattern of the FFT. Each node of a Boolean cube network of N nodes has a degree of log2N . For a large number of nodes the number of channels required at the chip boundary may be unfeasibly large with several nodes to a chip, and a network with slightly lower connectivity, such as Cube Connected Cycles networks, may be preferable. The communication system is the most critical resource in many high performance architectures, and its effective use imperative. We describe FFT algorithms that use both the storage bandwidth and the communication sys-tem optimally for an architecture such as the Connection Machine that has 65536 processors interconnected in a Boolean cube related network. We also describe the necessary data allocation, and the allocation and use of the twiddle factors.

The key to meeting the extremely high throughput requirements of future military signal and image processing systems is parallelism in algorithms and hardware. This paper will describe the implementation of a core set of algorithms on one possible hardware implementation, designed to achieve high speed and efficient parallelism. This approach and design procedure, while using currently available integrated circuit building blocks, is similar to how this type of processor will be developed in the future using VLSI.

MOSAIC is a heterogeneous multiprocessor architecture designed for a class of real-time signal processing applications. The architecture attempts to reuse both hardware modules and software routines such that fast system building can be achieved. Special purpose signal processing engines and general purpose computers are integrated with a high performance MOSAIC crossbar switch to execute multi-gigaflop applications. The integrated system is programmed in a graphical parallel language. The language supports many signal processing primitives that are essential in real-time applications. Most of these primitives are supported directly by dedicated hardware to enhance system throughput.