Probabilistic Complexity Analysis for a Class of ... - CiteSeerX
WINOGRAD AND NAWAB
1
Probabilistic Complexity Analysis for a Class of Approximate DFT Algorithms
Joseph M. Winograd and S. Hamid Nawab | To appear in J. VLSI Signal Processing, 1996. |
J. M. Winograd and S. H. Nawab are with the ECS Department, Boston University, Boston, MA, 02215. November 19, 1996
DRAFT
WINOGRAD AND NAWAB
2
Abstract We present a probabilistic complexity analysis of a class of multi-stage algorithms which incrementally re ne DFT approximations. Each stage of any algorithm in this class improves the results of the previous stage by a xed increment in one of three dimensions: SNR, frequency resolution, or frequency coverage. However, the complexity of each stage is probabilistically dependent upon certain characteristics of the input signal. Assuming that an algorithm has to be terminated before its arithmetic cost exceeds a given limit, we have formulated a method for predicting the probability of completion of each of the algorithm's stages. This analysis is useful for low-power and real-time applications where FFT algorithms cannot meet the speci ed limits on arithmetic cost.
I. Introduction
While the palette of transforms available to the DSP system designer continues to broaden, the utility of the DFT across a broad range of applications remains unparalleled. This fact can be attributed in part to the existence of a large class of FFT algorithms [1] which are able to evaluate an N -point DFT using O(N log N ) operations. In some systems, most notably those which operate in real-time and lowpower contexts, the use of the FFT can be precluded by strict design constraints regarding the number of arithmetic operations to be employed for DFT evaluation. For these applications, it is necessary to consider algorithms which compute approximations to the DFT in order to reduce computational cost. Many di erent approximate DFT algorithms have been proposed. The most well-known are the \pruning"-type algorithms which obtain computational eciency by excluding some subset of input and/or output points. Algorithms of this type include the FFT pruning algorithms [2]-[4], Goertzel's algorithm [5], and others [6] [7]. Advantages of pruning algorithms include the possibility of using the ecient FFT structure and the ease with which the error introduced through the approximation may be quanti ed. In contrast to pruning approaches, one may consider sacri cing the precision with which the DFT is computed. For example, such DFT approximations have been obtained using the summation by parts approach [8], the Poorman's approach [9], and the quantization and backward di erencing (QBD) approach [10]. We have recently introduced [11] a class of approximate DFT algorithms which use both QBD approximation and pruning. We refer to these algorithms as DFT incremental-re nement (DFT-IR) algorithms. November 19, 1996
DRAFT
WINOGRAD AND NAWAB
3
The term incremental re nement denotes the fact that each of these algorithms consists of multiple stages, where each stage improves upon the DFT approximation produced by the previous stage. The spectral degradation in the DFT approximation after each stage may be characterized in terms of commonly used input-independent metrics for spectral quality: SNR, frequency resolution, and frequency coverage. The arithmetic complexity of each stage, however, depends upon the nature of the input signal. This has two important implications: If an algorithm is to be terminated when the arithmetic cost reaches a speci ed level, the resulting
spectral degradation in the algorithm's output depends on the nature of the input signal. If an algorithm is to be terminated when the spectral degradation has been reduced to a speci ed
level, the arithmetic cost of the algorithm depends on the nature of the input signal. It follows that in order to select an algorithm for a particular application, we must characterize the dependence on the input signal in a more quantitative fashion. In this paper, we focus upon the fundamental issues underlying algorithm selection from the class of DFT-IR algorithms. In particular, we consider applications in which the algorithm is to be terminated after a xed arithmetic cost. Our starting point is to assume that the input signals in the application may be characterized by a Gaussian-distributed stationary process with a known autocorrelation. For such cases, we derive the probability of completing any particular algorithm stage, and thus producing a corresponding level of spectral degradation. This, for example, enables a system designer to identify what levels of spectral degradation the di erent algorithms may produce within the arithmetic bound with a speci ed probability P . The designer may then select an algorithm for which the associated spectral degradation is deemed acceptable. If none of the candidate spectral degradations is considered acceptable for the application, the designer may choose to lower the value of the probability P or, alternatively, to increase the limit on the arithmetic cost. The emphasis in this paper is on the probabilistic complexity analysis underlying the approach to algorithm selection we have just outlined. In Sec. II, we provide some background on the class of DFTIR algorithms. In particular, we point out that when the input is a stochastic process, the arithmetic complexity associated with the completion of each stage has a probabilistic component. The distribution November 19, 1996
DRAFT
WINOGRAD AND NAWAB
4
of this probabilistic component is presented and examined in Sec. III. Finally, in Sec. IV we present experimental veri cation of our theoretical predictions of the probability of completion of any given stage of DFT-IR algorithms. II. DFT-IR Algorithms
We have recently presented [11] a framework for obtaining algorithms that produce sequences of successively improving approximations to the N -point DFT of a real-valued xed-point binary input signal. Each algorithm in this framework is speci ed by the values of a set of control parameters which describe the manner in which the approximations are improved in successive stages of processing. The values of these control parameters also dictate the level of spectral degradation present in the approximations at each stage. The number of arithmetic operations required by each stage, however, is dependent on the input signal. This important characteristic justi es a probabilistic approach to the analysis of the algorithms' arithmetic complexity. In this section, we brie y outline the structure of these algorithms, present characterizations of the spectral degradation in their results, and discuss their computational requirements. A. Structure
Every DFT-IR algorithm can be viewed as a cascade of stages, each of which takes a DFT approximation
X^i?1 (k) and produces an improved approximation X^ i (k). Examples of this structure are illustrated using a block-diagram format in Fig. 1(a)-(c). The re nement process is \jump-started" with the computation of an initial approximation X^0 (k), de ned later in this section, whose computation is carried out using the block of type J in Fig. 1. Each subsequent stage performs one of three di erent updates, and each type of update improves the previous approximation in a di erent way. The blocks of type S in Fig. 1 perform SNR updates. That is, each improves the SNR of the previous approximation by performing additional
computation. Similarly, the blocks of types R and C represent resolution updates and coverage updates respectively. The speci c arithmetic operations that are performed in any particular update depend upon the sequence of blocks (not their order) that precede it. We indicate this dependence in Fig. 1 by denoting successive instances of a particular update block by S 0 , S 00 , etc. November 19, 1996
DRAFT
WINOGRAD AND NAWAB
5
Every unique sequence of updates corresponds to a di erent DFT-IR algorithm. For the N -point transform of a Q-bit input signal, the total number of di erent DFT-IR algorithms, which we denote by
NA , can be shown to be Q + 3N=2 ? 3)! NA = (Q ? (1)!( N=2 ? 1)!(N ? 1)!
(1)
We represent each of these algorithms using a set of control parameters, si , ri , and ci . Table 1 lists the control parameters associated with the sequence of stages shown in Fig. 1(a). For each i, the control parameter values essentially represent the number of updates of the corresponding type that are present up to and including the ith stage. For example, si ? 1 is equal to the number of SNR updates performed through the ith stage. The o set by one in each of the control parameters accounts for the initial approximation produced by the \jump-start" stage. The operation performed by each stage of a DFT-IR algorithm can be stated in terms of the values of the input data and the control parameters associated with that stage. We begin by de ning the three updates mathematically. Their implementation is discussed in Sec. II-B. The SNR update improves an approximation by the incorporation of an additional bit level of the input signal. It is de ned by 8 > > >
> > :
P ?1 X^i?1 (k) + nri=0 gsi (n)Gsi ;n (k); 1 k ci ;
X^i?1 (k);
otherwise;
(2)
where Gq;n (k) is de ned as 8 > > >
?e?j 2N kn =(1 ? e?j 2N k );
> > :21?q
q=1
e?j N kn =(1 ? e?j N k ); 2 q Q; 2
2
(3)
and gq (n) is the rst circular backward di erence of the bit vector xq (n), or 8 > > >
> > :xq (n) ? xq (n ? 1); 1 n N ? 1;
(4)
where xq (n) denotes the qth bit of the two's complement binary fraction representation of x(n). In this signal representation, the value of each element of the Q-bit signal x(n) is related to the corresponding November 19, 1996
DRAFT
WINOGRAD AND NAWAB
6
values its component bit vectors by:
x(n) = ?x1 (n) +
Q X q=2
21?q xq (n)
(5)
The resolution update improves the frequency resolution of the approximation by including an additional time sample of the input signal. It is de ned as 8 > > >
> > ^i?1 (k); :X
otherwise.
(6)
The coverage update improves the approximation by adding to it an additional frequency sample. It is de ned as
8 > > >
> > :
P P ?1 X^ i?1 (ci ) + sqi=1 nri=0 gq (n)Gq;n (ci ); k = ci ;
X^ i?1 (k);
otherwise.
(7)
The \jump-start" stage computes an initial approximation, X^0 (k). It is de ned as 8 > > >
(8) > > :0; otherwise. The amount of spectral degradation present in the ith successive DFT approximation can be characterized in terms of the values of the control parameters. As an alternative to the recursive update form,
X^i (k) can be expressed as X^ i (k) =
si rX i ?1 X q=1 n=0
gq (n)Gq;n (k); 1 k ci :
(9)
From this equation, it can be seen that ci dictates the spectral bandwidth, 2ci =N radians, over which
X^i (k) is evaluated. The approximate transformation truncates gq (n) to a length of ri samples, e ectively reducing the frequency resolution of the transformation so that not more than ri distinct frequencies can be resolved. The SNR of the approximate transform X^i (k) is reduced by signal quantization to approximately 6si dB. B. Arithmetic Complexity
The DFT-IR algorithms implement the update equations (2), (6), and (7) without multiplications using a technique [12] based on the summation of pre-computed vectors. The use of pre-computed partial results November 19, 1996
DRAFT
WINOGRAD AND NAWAB
7
to perform linear combinations is generally referred to as distributed arithmetic [13]. The application of various distributed arithmetic techniques to DFT processing has been considered by others [14] [15], though not in the context of approximate DFT algorithms. In the vector summation approach [12] used by the DFT-IR algorithms, the complex values of Gq;n (k) are stored in memory and are added or subtracted from X^ i?1 (k) according to the value of gq (n) as dictated by the update equations. All summations corresponding to gq (n) = 0 are skipped, resulting in a signi cant reduction in computation. The total number of real additions, i , required for evaluating all the stages up to and including the ith stage is
i = si ri + 2 (si ; ri )ci ;
(10)
where
(si ; ri ) =
si rX i ?1 X q=1 n=0
jgq (n)j:
(11)
The si ri term in (10) accounts for the backward di erencing operations required to produce gq (n) from
xq (n) over the region included through the ith stage of processing (recall that gq (n) is the backward di erenced vector de ned in Eq. (4)). The second term in (10) re ects the number of additions required to evaluate only those terms of the update equations for which gq (n) 6= 0. The quantity (si ; ri ), de ned in Eq. (11), is the total number of non-zero elements in the portion of the backward di erenced signal vectors gq (n) included through the ith stage. We consequently refer to it as the non-zero count for stage i. For notational convenience, we abbreviate the quantity (si ; ri ) by
i , though it should be understood that there is a dependence on the control parameters associated with the ith stage. The non-zero count is related to the input signal x(n) through Eqs. (4), (5), and (11). It represents the total signal-dependent contribution to the arithmetic cost of completing the ith stage of processing. For any two input signals the non-zero count may take on di erent values, resulting in a di erent arithmetic cost for performing the same sequence of stages on those signals. A complete characterization of the arithmetic complexity of the DFT-IR algorithms requires that the signal-dependence of the non-zero count be determined more precisely. Our approach to this problem is based on a probabilistic analysis. When the input signal x(n) is modeled as a stochastic process, the November 19, 1996
DRAFT
WINOGRAD AND NAWAB
8
non-zero count, i , is a random variable. The total arithmetic cost of completing all the stages up to and including the ith stage, i , is related to i through Eq. (10). Thus, if we wish to determine the probability distribution of the arithmetic cost for a class of signals, it is required that we rst obtain the probability distribution of the non-zero count. C. Spectral Degradation and Arithmetic Bounds
In the context of this paper, we are especially interested in determining the e ect of terminating any particular DFT-IR algorithm when its arithmetic cost reaches a speci ed bound B . As discussed in the introduction, we characterize this e ect in terms of the probabilities of completion associated with each of the individual algorithm stages. We de ne the probability of completion, Pi , of a DFT-IR algorithm to be the probability with which all stages up to and including stage i of that algorithm are completed using not more than B arithmetic operations. That probability can be expressed as
Pi = Prob fi B g
(12)
where i is the arithmetic complexity measure given in Eq. (10). In turn, this leads to the conclusion:
Pi = Prob i B ?2csi ri i
(13)
In order to determine the probability Pi , we once again need to determine the probability distribution of the non-zero count i . That is, we need to characterize the same random variable which arose in the context of the arithmetic complexity measure i . III. Probabilistic Complexity Analysis
The arithmetic complexity associated with the completion of the ith stage of any given DFT-IR algorithm has a probabilistic component. As discussed in the previous section, this component depends upon the random variable i (de ned in Eq. (11)). In this section, we begin by presenting and discussing1 an approximate probability distribution for i . We then illustrate some important properties of the mean 1
Further details of the derivation are supplied in the appendices.
November 19, 1996
DRAFT
WINOGRAD AND NAWAB
9
of this distribution. Finally, we discuss the implications of these results for analyzing the complexity of DFT-IR algorithms. A. Derivation of Distribution
We have derived an approximate expression for the probability distribution of i . While this expression is not in closed form, it can be evaluated at any point in polynomial time. The assumptions which lead to this expression being considered approximate are: (i) It is assumed that the values of gq (0), 1 q Q, make insigni cant contributions to i . (ii) The interactions2 between elements of gq (n) of third and higher order are assumed to be negligible. (iii) Any sample of gj (n) is is assumed to not be signi cantly correlated with any sample of gk (n) when
j 6= k. Each of these assumptions in some way facilitates a tractable analysis of the distribution of the nonzero count while introducing some error into its results. Assumption (i) is used because while gq (n0 ) =
xq (n0 ) ? xq (n0 ? 1) for n0 > 0, the sample gq (0) depends upon the di erence xq (0) ? xq (N ? 1). Our probabilistic analysis is simpli ed by ignoring this exception. This is reasonable because generally gq (0) is only one of many samples in gq (n) which contribute to i . Our second assumption is based on the observation that interactions between elements of gq (n) diminish with increasing order. Thus, restricting our consideration to the rst and second order moments can be expected to capture the most important characteristics of the distribution. The implications of this approximation have been investigated in detail previously. Bounds on its error are known [16] and it has been shown experimentally to be reasonable [17]. While we are not aware of any formal analysis to support assumption (iii), our empirical results on the probability of completion (as reported in Sec. IV) indicate that this assumption did not have a signi cant e ect. Our expression for the approximate probability distribution for i is:
p i (k) = Prob f i = kg = p i1 (k) p i2 (k) p isi (k)
(14)
We say that random variables interact when the expected value of their product is not equal to the product of their expected values. 2
November 19, 1996
DRAFT
WINOGRAD AND NAWAB
10
where "
i ?1 Rzz (q; n) ? 2q ri + rX ( r ? n ) i k q (1 ? q ) n=1 !2?w min(2 X;k) 2 r ? 2 ? i q p q (1 ? q ) w=max(0;k?ri +2) w k ? w
p iq (k) = kq (1 ? q )ri ?k
1 ? q p q (1 ? q )
!w 3 5
(15)
We use iq to denote the non-zero count over 0 n ri ? 1 of the backward di erence vector gq (n) and p iq (k) to denote its distribution. The quantities q and Rzz (q; n) respectively represent the non-zero probability and the autocorrelation of the random variable z (q; n) = jgq (n)j.3 They are related to the distribution of the stochastic process x~(n) from which the input signal x(n) is quantized as
q =
ZZ
(x;y): q (x)6=q (y)
and
8 ZZ > > > > > > > > (x;y): > > > q (x)6=q (y) > > ZZZ > > > >
px~(n?1)~x(n) (x; y) dx dy
px~(n?1)~x(n) (x; y) dx dy px~(n?1)~x(n)~x(n+1) (x; y; z ) dx dy dz
(x;y;z):
> q (x)6=q (y)6=q (z) > > ZZZZ > > > > > x~(?1)~x(0)~x(n?1)~x(n) ( > > > > > (w;x;y;z): > > > :q (w)6=q (x)^
p
(16)
n = 0; n = 1 _ n = N ? 1;
(17)
w; x; y; z ) dw dx dy dz 2 n N ? 2:
q (y)6=q (z)
The functions q (x) represent the input/output relation for the qth bit of a two's complement binary quantizer. They are de ned in Eq. (21) of Appendix A. The derivation of Eqs. (16) and (17) is also given in the appendices. Careful inspection of Eq. (15) reveals that the distribution of iq has a similar form to that of the binomial distribution. In fact, when there are no second order interactions between distinct elements of
gq (n), a binomial distribution with \probability of success" q is obtained. This equation can be viewed as providing the corrections to the binomial distribution that are required when the elements of gq (n) are correlated. 3
We note that (
) is made stationary in time by assumption (i).
z q; n
November 19, 1996
DRAFT
WINOGRAD AND NAWAB
11
We also note that the result given in Eq. (14) holds for any stationary process x~(n). Obtaining the values of q and Rzz (q; n) requires, however, that the distribution of that process be known. We consider only stationary zero-mean Gaussian processes, for which the distribution is completely speci ed by the autocorrelation function. Determination of q and Rzz (q; n) for these processes requires the evaluation of the multivariate normal integral. There are only a few special cases [18] for which closed form solutions exist, so for this task we rely on numerical methods [19] [20]. To illustrate the distribution of the non-zero count which results from our analysis, we consider as the input to a DFT-IR algorithm a zero-mean stationary Gaussian-distributed process with the autocorrelation shown in Fig. 2. This is the autocorrelation function associated with a reported long-term average spectrum of male speech [21]. The distribution obtained for i when only the most signi cant bit of the input signal is used is shown in Fig. 3(a). The distributions obtained when using the two, three, and four most signi cant bits are shown in Fig. 3(b), (c), and (d) respectively.4 We observe from the gure that as we include additional bits, the mean of the distribution tends to increase in proportion to its range. This is a general property of the distribution of i and it is explored next. B. Analysis of the Mean
When the ith stage of a DFT-IR algorithm consists of either a resolution update or an SNR update (i.e. ri > ri?1 or si > si?1 ), the number of elements of the backward di erence vectors gq (n) that are included in the approximation is increased. This results in the probability distribution of i di ering from the distribution of i?1 and, in particular, the mean value of the distribution must increase. In order to study the rate at which the mean value increases across stages, we consider the non-zero density ^i , which we de ne as
^i = i =si ri
(18)
The non-zero density represents the proportion of non-zero values over the portion of the backward differenced vectors under consideration, and takes on values in the range [0; 1]. In all cases, we illustrate the distribution obtained for a stage in which only the rst 32 samples of q ( ) are used (i.e. i = 32). 4
g
n
r
November 19, 1996
DRAFT
WINOGRAD AND NAWAB
12
Using the linearity property of expectation, the mean of the non-zero density can be shown to be: si X E [^ i ] = E s ri = s1 q : (19) i i i q=1 where E [] denotes the expectation operator. We see that E [^ i ] is independent of frequency resolution,
ri , and frequency coverage, ci , and that it only changes at stages for which additional bits of the input signal are added through an SNR update (i.e. si is increased). The nature of this change depends on the non-zero probability q , de ned in Eq. (16), for the backward di erence vectors that are added in that SNR update. We consider, then, the value of q . For the most signi cant backward di erenced vector, the integration in Eq. (16) can be evaluated [18] to yield a simple expression for 1 :
Rx~x~ (1) 1 = 1 arccos R x~x~ (0)
(20)
We use Rx~x~ (n) to denote the autocorrelation of the Gaussian input process x~(n). It should be noted that the non-zero probability depends only on the correlation between adjacent samples of the input process. This is also the case for the non-zero probability of the backward di erenced vectors associated with bits of the input signal of lesser signi cance. There is no closed form expression available for their evaluation, but high-precision numerical evaluation [19] is possible. An important property of the non-zero probability for the backward di erence vectors is illustrated by Fig. 4. This gure shows the non-zero probability associated with the four most signi cant vectors gq (n) as a function of the correlation between adjacent samples of x(n).5 It can be seen from this gure that for any positive autocorrelation, the non-zero probability in the backward di erence vectors gq (n) increases monotonically as one moves from the most signi cant towards the least signi cant bits. This implies that the mean of the non-zero density, given in Eq. (19), must also increase as additional bits are used through application of the SNR update. C. Implications for Complexity
One useful way of characterizing the arithmetic complexity associated with completing a given algorithm stage is by its expected arithmetic cost. This cost can be obtained by using the mean value of the non-zero Only positive correlations are shown. In cases where the rst autocorrelation value is negative, a frequency reversal technique [12] can be applied to make it positive. 5
November 19, 1996
DRAFT
WINOGRAD AND NAWAB
13
count in the cost analysis of Eq. (10). The results of the mean value analysis presented above suggests that the sequence in which the various updates are performed has an in uence on the rate at which the expected arithmetic cost increases. In particular, they imply that when an SNR update is performed, the rate of growth in the expected arithmetic cost across the stages that follow it will increase. To illustrate this property, we compare the expected arithmetic complexity of the stages of two di erent 256-point DFT-IR algorithms (Algorithm 1 and Algorithm 2). The control parameters, associated spectral degradation, and expected cost of four stages taken from each of these algorithms are given in Table 2. For this example, we assume that the signal under analysis has the same probability distribution as the one described in Sec. III-A. The correlation between adjacent samples of this process is 0:815. The upper half of Table 2 contains data corresponding to selected stages of Algorithm 1 for which frequency coverage and SNR are held constant, and the frequency resolution is doubled between each of the stages shown. The expected cost of completing each stage, which is determined by using the expected value of the non-zero count in Eq. (10), can be seen to increase with each stage. The non-zero density, however, remains constant due to the fact that no SNR updates are performed between these stages. The lower half of Table 2 contains similar information for four stages of Algorithm 2. Across the stages, the frequency coverage and frequency resolution are held constant while the SNR is increased in consecutive stage. The expected total cost for completing each of the stages is shown, as is the expected non-zero density. In comparison with Algorithm 1, the expected cost grows at a faster rate even though there are far fewer updates being performed. In fact, in the absence of SNR updates, Algorithm 1 does not reach the expected arithmetic cost of stage 73 of Algorithm 2 until more than 50 additional stages are performed. This illustrates how the presence of SNR updates increases the rate of growth of the expected completion cost for the stages that follow it. IV. Application
As discussed in the introduction of this paper, a system designer would be aided in the process of algorithm selection from the class of DFT-IR algorithms by being provided with the probabilities of completion of the various stages of those algorithms. In this section, we present experimental results to November 19, 1996
DRAFT
WINOGRAD AND NAWAB
14
support our theoretical predictions (based on the results from Sec. III) of the probability of completion of various DFT-IR algorithm stages. The algorithms were applied to Gaussian-distributed stationary data whose autocorrelation corresponds to a reported [21] long-term average for speech signals. For our experiments, we assume that the signal under analysis is contained in a 128-point frame which is to undergo a 256 point DFT. The signals were synthetically generated using the autocorrelation function shown in Fig. 2 sampled at a rate of 32 kHz. The variance of the input process was normalized to 0.25, representing a well-scaled input to a quantizer with amplitude range [?1; 1]. Signal quantization of 16bits was applied to obtain a signal suitable for DFT-IR analysis. The probabilities of completion for the stages of selected DFT-IR algorithms were determined empirically by monitoring the relative frequency with which those stages were completed in 50,000 Monte Carlo trials. The DFT-IR algorithms used in the analysis were the Algorithms 1 and 2 introduced in Sec. III-C. The stages of those algorithms for which completion statistics were tabulated are those given in Table 2. The probabilities of completion as determined by this experiment are listed in the fourth column of Table 3. The probabilities of completion for the same algorithm stages were also evaluated using the theoretical results of Secs. II-C and III-A. Numerical methods [19] [20] for the evaluation of the multivariate normal integral were employed to obtain the non-zero probability q and autocorrelations Rzz (q; n) as required by Eqs. (16)-(17). The probabilities of completion were computed by using the resulting moments in Eqs. (13), (14), and (15). The probabilities of completion obtained through the application of our theoretical results are listed in the third column of Table 3. Our experimental results indicate that reasonable estimates of the probabilities of completion are obtained through application of the theoretical analysis which we have presented. This tends to validate the approximations that were made in our derivations. The agreement between the theoretically determined probabilities and those that were obtained experimentally can be seen to diminish as the SNR of the approximation is increased. This e ect can be attributed to our assumption of independence between bit levels. Relaxation of this assumption is possible. However, the inclusion of all of the probabilistic interactions between di erent bit levels leads to a combinatorial growth in the complexity of obtaining the distribution of the non-zero count. This is described further in Appendix A. November 19, 1996
DRAFT
WINOGRAD AND NAWAB
15
V. Concluding Remarks
We have presented a probabilistic analysis of the arithmetic complexity of the DFT-IR class of approximate DFT algorithms. Through judicious use of suitable assumptions we have derived expressions for obtaining the probabilities of completion for the stages of any DFT-IR algorithm in the presence of a xed upper bound on arithmetic complexity. These results may be used by system designers to select appropriate DFT-IR algorithms for meeting the arithmetic complexity constraints of any given application. These results also lead to further interesting questions regarding the design process. For example, one may consider the development of procedures for eciently obtaining the DFT-IR algorithm whose stage sequence may be considered as optimal with respect to application speci c requirements. VI. Acknowledgments
This work was sponsored in part by the Department of the Navy, Oce of the Chief of Naval Research, contract number N00014-93-1-0686 as part of the Advanced Research Projects Agency's RASSP program. Appendices
A. Exact Distribution of the Non-Zero Count
In this appendix we derive the exact distribution of the non-zero count, i , when the signal under analysis is derived from a stochastic process. Our analysis is applicable to processes with general distributions. Assume the signal under analysis, x(n), to be the result of amplitude quantization of the discrete-time stochastic process x~(n). The values of the bit vectors xq (n) of which x(n) is comprised are related to this stochastic process by
xq (n) = q (~x(n))
(21)
where q (x) is the relation between the input value to the quantizer and the qth bit of its output. For two's complement binary PCM encoding, with the quantizer input range normalized to [?1; 1], that relation can be expressed as
8 > >
> :
November 19, 1996
u(?x);
q = 1; P q?1
u(x ? 1 + 21?q ) + 2i=0 ?2 u21?q (x + 1 ? 21?q (2i + 1)); q 2;
(22) DRAFT
WINOGRAD AND NAWAB
16
where u(x) is the unit step function and ul (x) is the unit pulse function which is equal to 1 for 0 x < l and 0 elsewhere. The quantizer relation q (x) is shown in Fig. 5 for several di erent values of q. Let px~(n)~x(m) () represent the joint probability of x~(n) and x~(m) and <x>y represent x modulo y. Considering Eqs. (4) and (21), the non-zero probability of z (q; n) is seen to be
q (n) = p(z (q; n) = 1) = pxq (N )xq (n) (xq (N ) 6= xq (n)) = px~(N )~x(n) (q (~x(N )) 6= q (~x(n))) =
ZZ
(x;y): q (x)6=q (y)
px~(N )~x(n) (x; y) dx dy:
(23)
The correlation between elements of z (q; n) can also be expressed in terms of the distribution of x~(n). That is,
Rzz (q; n; q0 ; n0 ) = E [z (q; n)z (q0; n0 )] = pz(q;n)z(q0 ;n0 ) (z (q; n) = 1 ^ z (q0 ; n0 ) = 1) 8 ZZ > > > > > > > > (x;y): > > >q (x)6=q (y)^ > > > q0 (x)6=q0 (y) > > ZZZ > > > >
(x;y;z):
px~(N )~x(n) (x; y) dx dy
n0 = n;
px~(N )~x(n)~x(n+1) (x; y; z ) dx dy dz
n0 ? n = 1 _ n0 ? n = N ? 1;
> q (x)6=q (y)^ > > > q0 (y)6=q0 (z) > > ZZZZ > > > > > > > > > > (w;x;y;z): > > > :q (w)6=q (x)^
px~(N )~x(n)~x(n0 ?1)~x(n0 ) (w; x; y; z ) dw dx dy dz 2 n0 ? n N ? 2;
q0 (y)6=q0 (z)
(24) with n0 n. Higher order moments of z (q; n) are similarly de ned, with the j th order moments of z (q; n) relying on 2j th order joint distributions of x~(n). Introducing y(q; n), which we de ne as
y(q; n) = pz (q; n) ? q (n) ; q (n)(1 ? q (n)) November 19, 1996
(25) DRAFT
WINOGRAD AND NAWAB
17
we normalize z (q; n) to zero mean and unit variance. The correlations between elements of y(q; n) are related to those of z (q; n) by 0 0 0 0 Ryy (q; n; q0 ; n0 ) = p Rzz (q; n; q ; n ) ? q (0 n) q (n ) 0 : q (n)(1 ? q (n)) q0 (n )(1 ? q0 (n ))
(26)
We denote by z the whole of z (q; n), de ned over the range 0 n ri ? 1 and 1 q si . Recall that each of the elements of z (q; n) are either 0 or 1 in value and i is the number of ones present in z. The joint distribution of the elements of z can be elegantly expressed in terms of q (n) and the moments of
y(q; n) using a little known result on correlated binomial random variables [16]: pz (z) = f (z)
si rY i ?1 Y q=1 n=0
q (n)z(q;n) (1 ? q (n))1?z(q;n)
(27)
where
f (z) = 1 +
X
qri +n
1
Probabilistic Complexity Analysis for a Class of Approximate DFT Algorithms
Joseph M. Winograd and S. Hamid Nawab | To appear in J. VLSI Signal Processing, 1996. |
J. M. Winograd and S. H. Nawab are with the ECS Department, Boston University, Boston, MA, 02215. November 19, 1996
DRAFT
WINOGRAD AND NAWAB
2
Abstract We present a probabilistic complexity analysis of a class of multi-stage algorithms which incrementally re ne DFT approximations. Each stage of any algorithm in this class improves the results of the previous stage by a xed increment in one of three dimensions: SNR, frequency resolution, or frequency coverage. However, the complexity of each stage is probabilistically dependent upon certain characteristics of the input signal. Assuming that an algorithm has to be terminated before its arithmetic cost exceeds a given limit, we have formulated a method for predicting the probability of completion of each of the algorithm's stages. This analysis is useful for low-power and real-time applications where FFT algorithms cannot meet the speci ed limits on arithmetic cost.
I. Introduction
While the palette of transforms available to the DSP system designer continues to broaden, the utility of the DFT across a broad range of applications remains unparalleled. This fact can be attributed in part to the existence of a large class of FFT algorithms [1] which are able to evaluate an N -point DFT using O(N log N ) operations. In some systems, most notably those which operate in real-time and lowpower contexts, the use of the FFT can be precluded by strict design constraints regarding the number of arithmetic operations to be employed for DFT evaluation. For these applications, it is necessary to consider algorithms which compute approximations to the DFT in order to reduce computational cost. Many di erent approximate DFT algorithms have been proposed. The most well-known are the \pruning"-type algorithms which obtain computational eciency by excluding some subset of input and/or output points. Algorithms of this type include the FFT pruning algorithms [2]-[4], Goertzel's algorithm [5], and others [6] [7]. Advantages of pruning algorithms include the possibility of using the ecient FFT structure and the ease with which the error introduced through the approximation may be quanti ed. In contrast to pruning approaches, one may consider sacri cing the precision with which the DFT is computed. For example, such DFT approximations have been obtained using the summation by parts approach [8], the Poorman's approach [9], and the quantization and backward di erencing (QBD) approach [10]. We have recently introduced [11] a class of approximate DFT algorithms which use both QBD approximation and pruning. We refer to these algorithms as DFT incremental-re nement (DFT-IR) algorithms. November 19, 1996
DRAFT
WINOGRAD AND NAWAB
3
The term incremental re nement denotes the fact that each of these algorithms consists of multiple stages, where each stage improves upon the DFT approximation produced by the previous stage. The spectral degradation in the DFT approximation after each stage may be characterized in terms of commonly used input-independent metrics for spectral quality: SNR, frequency resolution, and frequency coverage. The arithmetic complexity of each stage, however, depends upon the nature of the input signal. This has two important implications: If an algorithm is to be terminated when the arithmetic cost reaches a speci ed level, the resulting
spectral degradation in the algorithm's output depends on the nature of the input signal. If an algorithm is to be terminated when the spectral degradation has been reduced to a speci ed
level, the arithmetic cost of the algorithm depends on the nature of the input signal. It follows that in order to select an algorithm for a particular application, we must characterize the dependence on the input signal in a more quantitative fashion. In this paper, we focus upon the fundamental issues underlying algorithm selection from the class of DFT-IR algorithms. In particular, we consider applications in which the algorithm is to be terminated after a xed arithmetic cost. Our starting point is to assume that the input signals in the application may be characterized by a Gaussian-distributed stationary process with a known autocorrelation. For such cases, we derive the probability of completing any particular algorithm stage, and thus producing a corresponding level of spectral degradation. This, for example, enables a system designer to identify what levels of spectral degradation the di erent algorithms may produce within the arithmetic bound with a speci ed probability P . The designer may then select an algorithm for which the associated spectral degradation is deemed acceptable. If none of the candidate spectral degradations is considered acceptable for the application, the designer may choose to lower the value of the probability P or, alternatively, to increase the limit on the arithmetic cost. The emphasis in this paper is on the probabilistic complexity analysis underlying the approach to algorithm selection we have just outlined. In Sec. II, we provide some background on the class of DFTIR algorithms. In particular, we point out that when the input is a stochastic process, the arithmetic complexity associated with the completion of each stage has a probabilistic component. The distribution November 19, 1996
DRAFT
WINOGRAD AND NAWAB
4
of this probabilistic component is presented and examined in Sec. III. Finally, in Sec. IV we present experimental veri cation of our theoretical predictions of the probability of completion of any given stage of DFT-IR algorithms. II. DFT-IR Algorithms
We have recently presented [11] a framework for obtaining algorithms that produce sequences of successively improving approximations to the N -point DFT of a real-valued xed-point binary input signal. Each algorithm in this framework is speci ed by the values of a set of control parameters which describe the manner in which the approximations are improved in successive stages of processing. The values of these control parameters also dictate the level of spectral degradation present in the approximations at each stage. The number of arithmetic operations required by each stage, however, is dependent on the input signal. This important characteristic justi es a probabilistic approach to the analysis of the algorithms' arithmetic complexity. In this section, we brie y outline the structure of these algorithms, present characterizations of the spectral degradation in their results, and discuss their computational requirements. A. Structure
Every DFT-IR algorithm can be viewed as a cascade of stages, each of which takes a DFT approximation
X^i?1 (k) and produces an improved approximation X^ i (k). Examples of this structure are illustrated using a block-diagram format in Fig. 1(a)-(c). The re nement process is \jump-started" with the computation of an initial approximation X^0 (k), de ned later in this section, whose computation is carried out using the block of type J in Fig. 1. Each subsequent stage performs one of three di erent updates, and each type of update improves the previous approximation in a di erent way. The blocks of type S in Fig. 1 perform SNR updates. That is, each improves the SNR of the previous approximation by performing additional
computation. Similarly, the blocks of types R and C represent resolution updates and coverage updates respectively. The speci c arithmetic operations that are performed in any particular update depend upon the sequence of blocks (not their order) that precede it. We indicate this dependence in Fig. 1 by denoting successive instances of a particular update block by S 0 , S 00 , etc. November 19, 1996
DRAFT
WINOGRAD AND NAWAB
5
Every unique sequence of updates corresponds to a di erent DFT-IR algorithm. For the N -point transform of a Q-bit input signal, the total number of di erent DFT-IR algorithms, which we denote by
NA , can be shown to be Q + 3N=2 ? 3)! NA = (Q ? (1)!( N=2 ? 1)!(N ? 1)!
(1)
We represent each of these algorithms using a set of control parameters, si , ri , and ci . Table 1 lists the control parameters associated with the sequence of stages shown in Fig. 1(a). For each i, the control parameter values essentially represent the number of updates of the corresponding type that are present up to and including the ith stage. For example, si ? 1 is equal to the number of SNR updates performed through the ith stage. The o set by one in each of the control parameters accounts for the initial approximation produced by the \jump-start" stage. The operation performed by each stage of a DFT-IR algorithm can be stated in terms of the values of the input data and the control parameters associated with that stage. We begin by de ning the three updates mathematically. Their implementation is discussed in Sec. II-B. The SNR update improves an approximation by the incorporation of an additional bit level of the input signal. It is de ned by 8 > > >
> > :
P ?1 X^i?1 (k) + nri=0 gsi (n)Gsi ;n (k); 1 k ci ;
X^i?1 (k);
otherwise;
(2)
where Gq;n (k) is de ned as 8 > > >
?e?j 2N kn =(1 ? e?j 2N k );
> > :21?q
q=1
e?j N kn =(1 ? e?j N k ); 2 q Q; 2
2
(3)
and gq (n) is the rst circular backward di erence of the bit vector xq (n), or 8 > > >
> > :xq (n) ? xq (n ? 1); 1 n N ? 1;
(4)
where xq (n) denotes the qth bit of the two's complement binary fraction representation of x(n). In this signal representation, the value of each element of the Q-bit signal x(n) is related to the corresponding November 19, 1996
DRAFT
WINOGRAD AND NAWAB
6
values its component bit vectors by:
x(n) = ?x1 (n) +
Q X q=2
21?q xq (n)
(5)
The resolution update improves the frequency resolution of the approximation by including an additional time sample of the input signal. It is de ned as 8 > > >
> > ^i?1 (k); :X
otherwise.
(6)
The coverage update improves the approximation by adding to it an additional frequency sample. It is de ned as
8 > > >
> > :
P P ?1 X^ i?1 (ci ) + sqi=1 nri=0 gq (n)Gq;n (ci ); k = ci ;
X^ i?1 (k);
otherwise.
(7)
The \jump-start" stage computes an initial approximation, X^0 (k). It is de ned as 8 > > >
(8) > > :0; otherwise. The amount of spectral degradation present in the ith successive DFT approximation can be characterized in terms of the values of the control parameters. As an alternative to the recursive update form,
X^i (k) can be expressed as X^ i (k) =
si rX i ?1 X q=1 n=0
gq (n)Gq;n (k); 1 k ci :
(9)
From this equation, it can be seen that ci dictates the spectral bandwidth, 2ci =N radians, over which
X^i (k) is evaluated. The approximate transformation truncates gq (n) to a length of ri samples, e ectively reducing the frequency resolution of the transformation so that not more than ri distinct frequencies can be resolved. The SNR of the approximate transform X^i (k) is reduced by signal quantization to approximately 6si dB. B. Arithmetic Complexity
The DFT-IR algorithms implement the update equations (2), (6), and (7) without multiplications using a technique [12] based on the summation of pre-computed vectors. The use of pre-computed partial results November 19, 1996
DRAFT
WINOGRAD AND NAWAB
7
to perform linear combinations is generally referred to as distributed arithmetic [13]. The application of various distributed arithmetic techniques to DFT processing has been considered by others [14] [15], though not in the context of approximate DFT algorithms. In the vector summation approach [12] used by the DFT-IR algorithms, the complex values of Gq;n (k) are stored in memory and are added or subtracted from X^ i?1 (k) according to the value of gq (n) as dictated by the update equations. All summations corresponding to gq (n) = 0 are skipped, resulting in a signi cant reduction in computation. The total number of real additions, i , required for evaluating all the stages up to and including the ith stage is
i = si ri + 2 (si ; ri )ci ;
(10)
where
(si ; ri ) =
si rX i ?1 X q=1 n=0
jgq (n)j:
(11)
The si ri term in (10) accounts for the backward di erencing operations required to produce gq (n) from
xq (n) over the region included through the ith stage of processing (recall that gq (n) is the backward di erenced vector de ned in Eq. (4)). The second term in (10) re ects the number of additions required to evaluate only those terms of the update equations for which gq (n) 6= 0. The quantity (si ; ri ), de ned in Eq. (11), is the total number of non-zero elements in the portion of the backward di erenced signal vectors gq (n) included through the ith stage. We consequently refer to it as the non-zero count for stage i. For notational convenience, we abbreviate the quantity (si ; ri ) by
i , though it should be understood that there is a dependence on the control parameters associated with the ith stage. The non-zero count is related to the input signal x(n) through Eqs. (4), (5), and (11). It represents the total signal-dependent contribution to the arithmetic cost of completing the ith stage of processing. For any two input signals the non-zero count may take on di erent values, resulting in a di erent arithmetic cost for performing the same sequence of stages on those signals. A complete characterization of the arithmetic complexity of the DFT-IR algorithms requires that the signal-dependence of the non-zero count be determined more precisely. Our approach to this problem is based on a probabilistic analysis. When the input signal x(n) is modeled as a stochastic process, the November 19, 1996
DRAFT
WINOGRAD AND NAWAB
8
non-zero count, i , is a random variable. The total arithmetic cost of completing all the stages up to and including the ith stage, i , is related to i through Eq. (10). Thus, if we wish to determine the probability distribution of the arithmetic cost for a class of signals, it is required that we rst obtain the probability distribution of the non-zero count. C. Spectral Degradation and Arithmetic Bounds
In the context of this paper, we are especially interested in determining the e ect of terminating any particular DFT-IR algorithm when its arithmetic cost reaches a speci ed bound B . As discussed in the introduction, we characterize this e ect in terms of the probabilities of completion associated with each of the individual algorithm stages. We de ne the probability of completion, Pi , of a DFT-IR algorithm to be the probability with which all stages up to and including stage i of that algorithm are completed using not more than B arithmetic operations. That probability can be expressed as
Pi = Prob fi B g
(12)
where i is the arithmetic complexity measure given in Eq. (10). In turn, this leads to the conclusion:
Pi = Prob i B ?2csi ri i
(13)
In order to determine the probability Pi , we once again need to determine the probability distribution of the non-zero count i . That is, we need to characterize the same random variable which arose in the context of the arithmetic complexity measure i . III. Probabilistic Complexity Analysis
The arithmetic complexity associated with the completion of the ith stage of any given DFT-IR algorithm has a probabilistic component. As discussed in the previous section, this component depends upon the random variable i (de ned in Eq. (11)). In this section, we begin by presenting and discussing1 an approximate probability distribution for i . We then illustrate some important properties of the mean 1
Further details of the derivation are supplied in the appendices.
November 19, 1996
DRAFT
WINOGRAD AND NAWAB
9
of this distribution. Finally, we discuss the implications of these results for analyzing the complexity of DFT-IR algorithms. A. Derivation of Distribution
We have derived an approximate expression for the probability distribution of i . While this expression is not in closed form, it can be evaluated at any point in polynomial time. The assumptions which lead to this expression being considered approximate are: (i) It is assumed that the values of gq (0), 1 q Q, make insigni cant contributions to i . (ii) The interactions2 between elements of gq (n) of third and higher order are assumed to be negligible. (iii) Any sample of gj (n) is is assumed to not be signi cantly correlated with any sample of gk (n) when
j 6= k. Each of these assumptions in some way facilitates a tractable analysis of the distribution of the nonzero count while introducing some error into its results. Assumption (i) is used because while gq (n0 ) =
xq (n0 ) ? xq (n0 ? 1) for n0 > 0, the sample gq (0) depends upon the di erence xq (0) ? xq (N ? 1). Our probabilistic analysis is simpli ed by ignoring this exception. This is reasonable because generally gq (0) is only one of many samples in gq (n) which contribute to i . Our second assumption is based on the observation that interactions between elements of gq (n) diminish with increasing order. Thus, restricting our consideration to the rst and second order moments can be expected to capture the most important characteristics of the distribution. The implications of this approximation have been investigated in detail previously. Bounds on its error are known [16] and it has been shown experimentally to be reasonable [17]. While we are not aware of any formal analysis to support assumption (iii), our empirical results on the probability of completion (as reported in Sec. IV) indicate that this assumption did not have a signi cant e ect. Our expression for the approximate probability distribution for i is:
p i (k) = Prob f i = kg = p i1 (k) p i2 (k) p isi (k)
(14)
We say that random variables interact when the expected value of their product is not equal to the product of their expected values. 2
November 19, 1996
DRAFT
WINOGRAD AND NAWAB
10
where "
i ?1 Rzz (q; n) ? 2q ri + rX ( r ? n ) i k q (1 ? q ) n=1 !2?w min(2 X;k) 2 r ? 2 ? i q p q (1 ? q ) w=max(0;k?ri +2) w k ? w
p iq (k) = kq (1 ? q )ri ?k
1 ? q p q (1 ? q )
!w 3 5
(15)
We use iq to denote the non-zero count over 0 n ri ? 1 of the backward di erence vector gq (n) and p iq (k) to denote its distribution. The quantities q and Rzz (q; n) respectively represent the non-zero probability and the autocorrelation of the random variable z (q; n) = jgq (n)j.3 They are related to the distribution of the stochastic process x~(n) from which the input signal x(n) is quantized as
q =
ZZ
(x;y): q (x)6=q (y)
and
8 ZZ > > > > > > > > (x;y): > > > q (x)6=q (y) > > ZZZ > > > >
px~(n?1)~x(n) (x; y) dx dy
px~(n?1)~x(n) (x; y) dx dy px~(n?1)~x(n)~x(n+1) (x; y; z ) dx dy dz
(x;y;z):
> q (x)6=q (y)6=q (z) > > ZZZZ > > > > > x~(?1)~x(0)~x(n?1)~x(n) ( > > > > > (w;x;y;z): > > > :q (w)6=q (x)^
p
(16)
n = 0; n = 1 _ n = N ? 1;
(17)
w; x; y; z ) dw dx dy dz 2 n N ? 2:
q (y)6=q (z)
The functions q (x) represent the input/output relation for the qth bit of a two's complement binary quantizer. They are de ned in Eq. (21) of Appendix A. The derivation of Eqs. (16) and (17) is also given in the appendices. Careful inspection of Eq. (15) reveals that the distribution of iq has a similar form to that of the binomial distribution. In fact, when there are no second order interactions between distinct elements of
gq (n), a binomial distribution with \probability of success" q is obtained. This equation can be viewed as providing the corrections to the binomial distribution that are required when the elements of gq (n) are correlated. 3
We note that (
) is made stationary in time by assumption (i).
z q; n
November 19, 1996
DRAFT
WINOGRAD AND NAWAB
11
We also note that the result given in Eq. (14) holds for any stationary process x~(n). Obtaining the values of q and Rzz (q; n) requires, however, that the distribution of that process be known. We consider only stationary zero-mean Gaussian processes, for which the distribution is completely speci ed by the autocorrelation function. Determination of q and Rzz (q; n) for these processes requires the evaluation of the multivariate normal integral. There are only a few special cases [18] for which closed form solutions exist, so for this task we rely on numerical methods [19] [20]. To illustrate the distribution of the non-zero count which results from our analysis, we consider as the input to a DFT-IR algorithm a zero-mean stationary Gaussian-distributed process with the autocorrelation shown in Fig. 2. This is the autocorrelation function associated with a reported long-term average spectrum of male speech [21]. The distribution obtained for i when only the most signi cant bit of the input signal is used is shown in Fig. 3(a). The distributions obtained when using the two, three, and four most signi cant bits are shown in Fig. 3(b), (c), and (d) respectively.4 We observe from the gure that as we include additional bits, the mean of the distribution tends to increase in proportion to its range. This is a general property of the distribution of i and it is explored next. B. Analysis of the Mean
When the ith stage of a DFT-IR algorithm consists of either a resolution update or an SNR update (i.e. ri > ri?1 or si > si?1 ), the number of elements of the backward di erence vectors gq (n) that are included in the approximation is increased. This results in the probability distribution of i di ering from the distribution of i?1 and, in particular, the mean value of the distribution must increase. In order to study the rate at which the mean value increases across stages, we consider the non-zero density ^i , which we de ne as
^i = i =si ri
(18)
The non-zero density represents the proportion of non-zero values over the portion of the backward differenced vectors under consideration, and takes on values in the range [0; 1]. In all cases, we illustrate the distribution obtained for a stage in which only the rst 32 samples of q ( ) are used (i.e. i = 32). 4
g
n
r
November 19, 1996
DRAFT
WINOGRAD AND NAWAB
12
Using the linearity property of expectation, the mean of the non-zero density can be shown to be: si X E [^ i ] = E s ri = s1 q : (19) i i i q=1 where E [] denotes the expectation operator. We see that E [^ i ] is independent of frequency resolution,
ri , and frequency coverage, ci , and that it only changes at stages for which additional bits of the input signal are added through an SNR update (i.e. si is increased). The nature of this change depends on the non-zero probability q , de ned in Eq. (16), for the backward di erence vectors that are added in that SNR update. We consider, then, the value of q . For the most signi cant backward di erenced vector, the integration in Eq. (16) can be evaluated [18] to yield a simple expression for 1 :
Rx~x~ (1) 1 = 1 arccos R x~x~ (0)
(20)
We use Rx~x~ (n) to denote the autocorrelation of the Gaussian input process x~(n). It should be noted that the non-zero probability depends only on the correlation between adjacent samples of the input process. This is also the case for the non-zero probability of the backward di erenced vectors associated with bits of the input signal of lesser signi cance. There is no closed form expression available for their evaluation, but high-precision numerical evaluation [19] is possible. An important property of the non-zero probability for the backward di erence vectors is illustrated by Fig. 4. This gure shows the non-zero probability associated with the four most signi cant vectors gq (n) as a function of the correlation between adjacent samples of x(n).5 It can be seen from this gure that for any positive autocorrelation, the non-zero probability in the backward di erence vectors gq (n) increases monotonically as one moves from the most signi cant towards the least signi cant bits. This implies that the mean of the non-zero density, given in Eq. (19), must also increase as additional bits are used through application of the SNR update. C. Implications for Complexity
One useful way of characterizing the arithmetic complexity associated with completing a given algorithm stage is by its expected arithmetic cost. This cost can be obtained by using the mean value of the non-zero Only positive correlations are shown. In cases where the rst autocorrelation value is negative, a frequency reversal technique [12] can be applied to make it positive. 5
November 19, 1996
DRAFT
WINOGRAD AND NAWAB
13
count in the cost analysis of Eq. (10). The results of the mean value analysis presented above suggests that the sequence in which the various updates are performed has an in uence on the rate at which the expected arithmetic cost increases. In particular, they imply that when an SNR update is performed, the rate of growth in the expected arithmetic cost across the stages that follow it will increase. To illustrate this property, we compare the expected arithmetic complexity of the stages of two di erent 256-point DFT-IR algorithms (Algorithm 1 and Algorithm 2). The control parameters, associated spectral degradation, and expected cost of four stages taken from each of these algorithms are given in Table 2. For this example, we assume that the signal under analysis has the same probability distribution as the one described in Sec. III-A. The correlation between adjacent samples of this process is 0:815. The upper half of Table 2 contains data corresponding to selected stages of Algorithm 1 for which frequency coverage and SNR are held constant, and the frequency resolution is doubled between each of the stages shown. The expected cost of completing each stage, which is determined by using the expected value of the non-zero count in Eq. (10), can be seen to increase with each stage. The non-zero density, however, remains constant due to the fact that no SNR updates are performed between these stages. The lower half of Table 2 contains similar information for four stages of Algorithm 2. Across the stages, the frequency coverage and frequency resolution are held constant while the SNR is increased in consecutive stage. The expected total cost for completing each of the stages is shown, as is the expected non-zero density. In comparison with Algorithm 1, the expected cost grows at a faster rate even though there are far fewer updates being performed. In fact, in the absence of SNR updates, Algorithm 1 does not reach the expected arithmetic cost of stage 73 of Algorithm 2 until more than 50 additional stages are performed. This illustrates how the presence of SNR updates increases the rate of growth of the expected completion cost for the stages that follow it. IV. Application
As discussed in the introduction of this paper, a system designer would be aided in the process of algorithm selection from the class of DFT-IR algorithms by being provided with the probabilities of completion of the various stages of those algorithms. In this section, we present experimental results to November 19, 1996
DRAFT
WINOGRAD AND NAWAB
14
support our theoretical predictions (based on the results from Sec. III) of the probability of completion of various DFT-IR algorithm stages. The algorithms were applied to Gaussian-distributed stationary data whose autocorrelation corresponds to a reported [21] long-term average for speech signals. For our experiments, we assume that the signal under analysis is contained in a 128-point frame which is to undergo a 256 point DFT. The signals were synthetically generated using the autocorrelation function shown in Fig. 2 sampled at a rate of 32 kHz. The variance of the input process was normalized to 0.25, representing a well-scaled input to a quantizer with amplitude range [?1; 1]. Signal quantization of 16bits was applied to obtain a signal suitable for DFT-IR analysis. The probabilities of completion for the stages of selected DFT-IR algorithms were determined empirically by monitoring the relative frequency with which those stages were completed in 50,000 Monte Carlo trials. The DFT-IR algorithms used in the analysis were the Algorithms 1 and 2 introduced in Sec. III-C. The stages of those algorithms for which completion statistics were tabulated are those given in Table 2. The probabilities of completion as determined by this experiment are listed in the fourth column of Table 3. The probabilities of completion for the same algorithm stages were also evaluated using the theoretical results of Secs. II-C and III-A. Numerical methods [19] [20] for the evaluation of the multivariate normal integral were employed to obtain the non-zero probability q and autocorrelations Rzz (q; n) as required by Eqs. (16)-(17). The probabilities of completion were computed by using the resulting moments in Eqs. (13), (14), and (15). The probabilities of completion obtained through the application of our theoretical results are listed in the third column of Table 3. Our experimental results indicate that reasonable estimates of the probabilities of completion are obtained through application of the theoretical analysis which we have presented. This tends to validate the approximations that were made in our derivations. The agreement between the theoretically determined probabilities and those that were obtained experimentally can be seen to diminish as the SNR of the approximation is increased. This e ect can be attributed to our assumption of independence between bit levels. Relaxation of this assumption is possible. However, the inclusion of all of the probabilistic interactions between di erent bit levels leads to a combinatorial growth in the complexity of obtaining the distribution of the non-zero count. This is described further in Appendix A. November 19, 1996
DRAFT
WINOGRAD AND NAWAB
15
V. Concluding Remarks
We have presented a probabilistic analysis of the arithmetic complexity of the DFT-IR class of approximate DFT algorithms. Through judicious use of suitable assumptions we have derived expressions for obtaining the probabilities of completion for the stages of any DFT-IR algorithm in the presence of a xed upper bound on arithmetic complexity. These results may be used by system designers to select appropriate DFT-IR algorithms for meeting the arithmetic complexity constraints of any given application. These results also lead to further interesting questions regarding the design process. For example, one may consider the development of procedures for eciently obtaining the DFT-IR algorithm whose stage sequence may be considered as optimal with respect to application speci c requirements. VI. Acknowledgments
This work was sponsored in part by the Department of the Navy, Oce of the Chief of Naval Research, contract number N00014-93-1-0686 as part of the Advanced Research Projects Agency's RASSP program. Appendices
A. Exact Distribution of the Non-Zero Count
In this appendix we derive the exact distribution of the non-zero count, i , when the signal under analysis is derived from a stochastic process. Our analysis is applicable to processes with general distributions. Assume the signal under analysis, x(n), to be the result of amplitude quantization of the discrete-time stochastic process x~(n). The values of the bit vectors xq (n) of which x(n) is comprised are related to this stochastic process by
xq (n) = q (~x(n))
(21)
where q (x) is the relation between the input value to the quantizer and the qth bit of its output. For two's complement binary PCM encoding, with the quantizer input range normalized to [?1; 1], that relation can be expressed as
8 > >
> :
November 19, 1996
u(?x);
q = 1; P q?1
u(x ? 1 + 21?q ) + 2i=0 ?2 u21?q (x + 1 ? 21?q (2i + 1)); q 2;
(22) DRAFT
WINOGRAD AND NAWAB
16
where u(x) is the unit step function and ul (x) is the unit pulse function which is equal to 1 for 0 x < l and 0 elsewhere. The quantizer relation q (x) is shown in Fig. 5 for several di erent values of q. Let px~(n)~x(m) () represent the joint probability of x~(n) and x~(m) and <x>y represent x modulo y. Considering Eqs. (4) and (21), the non-zero probability of z (q; n) is seen to be
q (n) = p(z (q; n) = 1) = pxq (N )xq (n) (xq (N ) 6= xq (n)) = px~(N )~x(n) (q (~x(N )) 6= q (~x(n))) =
ZZ
(x;y): q (x)6=q (y)
px~(N )~x(n) (x; y) dx dy:
(23)
The correlation between elements of z (q; n) can also be expressed in terms of the distribution of x~(n). That is,
Rzz (q; n; q0 ; n0 ) = E [z (q; n)z (q0; n0 )] = pz(q;n)z(q0 ;n0 ) (z (q; n) = 1 ^ z (q0 ; n0 ) = 1) 8 ZZ > > > > > > > > (x;y): > > >q (x)6=q (y)^ > > > q0 (x)6=q0 (y) > > ZZZ > > > >
(x;y;z):
px~(N )~x(n) (x; y) dx dy
n0 = n;
px~(N )~x(n)~x(n+1) (x; y; z ) dx dy dz
n0 ? n = 1 _ n0 ? n = N ? 1;
> q (x)6=q (y)^ > > > q0 (y)6=q0 (z) > > ZZZZ > > > > > > > > > > (w;x;y;z): > > > :q (w)6=q (x)^
px~(N )~x(n)~x(n0 ?1)~x(n0 ) (w; x; y; z ) dw dx dy dz 2 n0 ? n N ? 2;
q0 (y)6=q0 (z)
(24) with n0 n. Higher order moments of z (q; n) are similarly de ned, with the j th order moments of z (q; n) relying on 2j th order joint distributions of x~(n). Introducing y(q; n), which we de ne as
y(q; n) = pz (q; n) ? q (n) ; q (n)(1 ? q (n)) November 19, 1996
(25) DRAFT
WINOGRAD AND NAWAB
17
we normalize z (q; n) to zero mean and unit variance. The correlations between elements of y(q; n) are related to those of z (q; n) by 0 0 0 0 Ryy (q; n; q0 ; n0 ) = p Rzz (q; n; q ; n ) ? q (0 n) q (n ) 0 : q (n)(1 ? q (n)) q0 (n )(1 ? q0 (n ))
(26)
We denote by z the whole of z (q; n), de ned over the range 0 n ri ? 1 and 1 q si . Recall that each of the elements of z (q; n) are either 0 or 1 in value and i is the number of ones present in z. The joint distribution of the elements of z can be elegantly expressed in terms of q (n) and the moments of
y(q; n) using a little known result on correlated binomial random variables [16]: pz (z) = f (z)
si rY i ?1 Y q=1 n=0
q (n)z(q;n) (1 ? q (n))1?z(q;n)
(27)
where
f (z) = 1 +
X
qri +n
