Convolution Algorithms

Fast Convolution by the FFT One of the main applications of the FFT is to do convolution more efficiently than the direct calculation from the definition which is: y ( n ) = ∑ m = 0 n h ( m ) x ( n - m ) which can also be written as: y ( n ) = ∑ m = 0 n x ( m ) h ( n - m ) This is often used to filter a signal x(n) with a filter whose impulse response is h(n). Each output value y(n) requires N multiplications and N-1 additions if y(n) has N terms. So, for N output values, on the order of N2 arithmetic operations are required. Because the DFT converts convolution to multiplication: D F T { y ( n ) } = D F T { h ( n ) } D F T { x ( n ) } can be calculated with the FFT and bring the order of arithmetic operations down to Nlog(N) which can be significate with large N. This approach, which is called “fast convolutions", is a form of block processing since a whole block of x(n) must be available to calculate even one output value, y(n). So, a time delay of one block length is always required. Another problem is the filtering use of convolution is usually non-cyclic and the convolution implemented with the DFT is cyclic. This is dealt with by appending zeros to x(n) and h(n) such that the output of the cyclic convolution gives one block of the output of the desired non-cyclic convolution. For filtering and some other applications, one want “on going" convolution where the filter response h(n) may be finite in length or duration, but the input x(n) is of arbitrary length. Two methods have traditionally used to break the input into blocks and use the FFT to convolve the block so that the output that would have been calculated by directly implementing or can be constructed efficiently. These are called “overlap-add" and “over-lap save".

Fast Convolution by Overlap-Add In order to use the FFT to convolve (or filter) a long input sequence x(n) with a finite length-M impulse response, h(n), we partition the input sequence in segments or blocks of length L. Because convolution (or filtering) is linear, the output is a linear sum of the result of convolving the first block with h(n) plus the result of convolving the second block with h(n), plus the rest. Each of these block convolutions can be calculated by using the FFT. The output is the inverse FFT of the product of the FFT of x(n) and the FFT of h(n). Since the number of arithmetic operation to calculate the convolution directly is on the order of M2 and, if done with the FFT, is on the order of Mlog(M), there can be a great savings by using the FFT for large M.The reason this procedure is not totally straightforward, is the length of the output of convolving a length-L block with a length-M filter is of length L+M-1. This means the output blocks cannot simply be concatenated but must be overlapped and added, hence the name for this algorithm is “Overlap-Add".The second issue that must be taken into account is the fact that the overlap-add steps need non-cyclic convolution and convolution by the FFT is cyclic. This is easily handled by appending L-1 zeros to the impulse response and M-1 zeros to each input block so that all FFTs are of length M+L-1. This means there is no aliasing and the implemented cyclic convolution gives the same output as the desired non-cyclic convolution.The savings in arithmetic can be considerable when implementing convolution or performing FIR digital filtering. However, there are two penalties. The use of blocks introduces a delay of one block length. None of the first block of output can be calculated until all of the first block of input is available. This is not a problem for “off line" or “batch" processing but can be serious for real-time processing. The second penalty is the memory required to store and process the blocks. The continuing reduction of memory cost often removes this problem.The efficiency in terms of number of arithmetic operations per output point increases for large blocks because of the Mlog(M) requirements of the FFT. However, the blocks become very large (L>>M), much of the input block will be the appended zeros and efficiency is lost. For any particular application, taking the particular filter and FFT algorithm being used and the particular hardware being used, a plot of efficiency vs. block length, L should be made and L chosen to maximize efficiency given any other constraints that are applicable.Usually, the block convolutions are done by the FFT, but they could be done by any efficient, finite length method. One could use “rectangular transforms" or “number-theoretic transforms". A generalization of this method is presented later in the notes.

Fast Convolution by Overlap-Save An alternative approach to the Overlap-Add can be developed by starting with segmenting the output rather than the input. If one considers the calculation of a block of output, it is seen that not only the corresponding input block is needed, but part of the preceding input block also needed. Indeed, one can show that a length M+L-1 segment of the input is needed for each output block. So, one saves the last part of the preceding block and concatenates it with the current input block, then convolves that with h(n) to calculate the current output

Block Processing, a Generalization of Overlap Methods Convolution is intimately related to the DFT. It was shown in The DFT as Convolution or Filtering that a prime length DFT could be converted to cyclic convolution. It has been long known that convolution can be calculated by multiplying the DFTs of signals. An important question is what is the fastest method for calculating digital convolution. There are several methods that each have some advantage. The earliest method for fast convolution was the use of sectioning with overlap-add or overlap-save and the FFT , , . In most cases the convolution is of real data and, therefore, real-data FFTs should be used. That approach is still probably the fastest method for longer convolution on a general purpose computer or microprocessor. The shorter convolutions should simply be calculated directly.

Introduction The partitioning of long or infinite strings of data into shorter sections or blocks has been used to allow application of the FFT to realize on-going or continuous convolution , . This section develops the idea of block processing and shows that it is a generalization of the overlap-add and overlap-save methods , . They further generalize the idea to a multidimensional formulation of convolution , . Moving in the opposite direction, it is shown that, rather than partitioning a string of scalars into blocks and then into blocks of blocks, one can partition a scalar number into blocks of bits and then include the operation of multiplication in the signal processing formulation. This is called distributed arithmetic and, since it describes operations at the bit level, is completely general. These notes try to present a coherent development of these ideas.

Block Signal Processing In this section the usual convolution and recursion that implements FIR and IIR discrete-time filters are reformulated in terms of vectors and matrices. Because the same data is partitioned and grouped in a variety of ways, it is important to have a consistent notation in order to be clear. The nth element of a data sequence is expressed h(n) or, in some cases to simplify, hn. A block or finite length column vector is denoted h̲n with n indicating the nth block or section of a longer vector. A matrix, square or rectangular, is indicated by an upper case letter such as H with a subscript if appropriate.

Block Convolution The operation of a finite impulse response (FIR) filter is described by a finite convolution as y ( n ) = ∑ k = 0 L - 1 h ( k ) x ( n - k ) where x(n) is causal, h(n) is causal and of length L, and the time index n goes from zero to infinity or some large value. With a change of index variables this becomes y ( n ) = ∑ k = 0 n h ( n - k ) x ( k ) which can be expressed as a matrix operation by y 0 y 1 y 2 ⋮ = h 0 0 0 ⋯ 0 h 1 h 0 0 h 2 h 1 h 0 ⋮ ⋮ x 0 x 1 x 2 ⋮ . The H matrix of impulse response values is partitioned into N by N square sub matrices and the X and Y vectors are partitioned into length-N blocks or sections. This is illustrated for N=3 by H 0 = h 0 0 0 h 1 h 0 0 h 2 h 1 h 0 H 1 = h 3 h 2 h 1 h 4 h 3 h 2 h 5 h 4 h 3 etc. x ̲ 0 = x 0 x 1 x 2 x ̲ 1 = x 3 x 4 x 5 y ̲ 0 = y 0 y 1 y 2 etc. Substituting these definitions into gives y ̲ 0 y ̲ 1 y ̲ 2 ⋮ = H 0 0 0 ⋯ 0 H 1 H 0 0 H 2 H 1 H 0 ⋮ ⋮ x ̲ 0 x ̲ 1 x ̲ 2 ⋮ The general expression for the nth output block is y ̲ n = ∑ k = 0 n H n - k x ̲ k which is a vector or block convolution. Since the matrix-vector multiplication within the block convolution is itself a convolution, is a sort of convolution of convolutions and the finite length matrix-vector multiplication can be carried out using the FFT or other fast convolution methods. The equation for one output block can be written as the product y ̲ 2 = H 2 H 1 H 0 x ̲ 0 x ̲ 1 x ̲ 2 and the effects of one input block can be written H 0 H 1 H 2 x ̲ 1 = y ̲ 0 y ̲ 1 y ̲ 2 . These are generalize statements of overlap save and overlap add , . The block length can be longer, shorter, or equal to the filter length.

Block Recursion Although less well-known, IIR filters can be implemented with block processing , , , , . The block form of an IIR filter is developed in much the same way as for the block convolution implementation of the FIR filter. The general constant coefficient difference equation which describes an IIR filter with recursive coefficients al, convolution coefficients bk, input signal x(n), and output signal y(n) is given by y ( n ) = ∑ l = 1 N - 1 a l y n - l + ∑ k = 0 M - 1 b k x n - k using both functional notation and subscripts, depending on which is easier and clearer. The impulse response h(n) is h ( n ) = ∑ l = 1 N - 1 a l h ( n - l ) + ∑ k = 0 M - 1 b k δ ( n - k ) which can be written in matrix operator form 1 0 0 ⋯ 0 a 1 1 0 a 2 a 1 1 a 3 a 2 a 1 0 a 3 a 2 ⋮ ⋮ h 0 h 1 h 2 h 3 h 4 ⋮ = b 0 b 1 b 2 b 3 0 ⋮ In terms of N by N submatrices and length-N blocks, this becomes A 0 0 0 ⋯ 0 A 1 A 0 0 0 A 1 A 0 ⋮ ⋮ h ̲ 0 h ̲ 1 h ̲ 2 ⋮ = b ̲ 0 b ̲ 1 0 ⋮ From this formulation, a block recursive equation can be written that will generate the impulse response block by block. A 0 h ̲ n + A 1 h ̲ n - 1 = 0 for n ≥ 2 h ̲ n = - A 0 - 1 A 1 h ̲ n - 1 = K h ̲ n - 1 for n ≥ 2 with initial conditions given by h ̲ 1 = - A 0 - 1 A 1 A 0 - 1 b ̲ 0 + A 0 - 1 b ̲ 1 This can also be written to generate the square partitions of the impulse response matrix by H n = K H n - 1 for n ≥ 2 with initial conditions given by H 1 = K A 0 - 1 B 0 + A 0 - 1 B 1 ane K=-A0-1A1. This recursively generates square submatrices of H similar to those defined in and and shows the basic structure of the dynamic system. Next, we develop the recursive formulation for a general input as described by the scalar difference equation and in matrix operator form by 1 0 0 ⋯ 0 a 1 1 0 a 2 a 1 1 a 3 a 2 a 1 0 a 3 a 2 ⋮ ⋮ y 0 y 1 y 2 y 3 y 4 ⋮ = b 0 0 0 ⋯ 0 b 1 b 0 0 b 2 b 1 b 0 0 b 2 b 1 0 0 b 2 ⋮ ⋮ x 0 x 1 x 2 x 3 x 4 ⋮ which, after substituting the definitions of the sub matrices and assuming the block length is larger than the order of the numerator or denominator, becomes A 0 0 0 ⋯ 0 A 1 A 0 0 0 A 1 A 0 ⋮ ⋮ y ̲ 0 y ̲ 1 y ̲ 2 ⋮ = B 0 0 0 ⋯ 0 B 1 B 0 0 0 B 1 B 0 ⋮ ⋮ x ̲ 0 x ̲ 1 x ̲ 2 ⋮ . From the partitioned rows of , one can write the block recursive relation A 0 y ̲ n + 1 + A 1 y ̲ n = B 0 x ̲ n + 1 + B 1 x ̲ n Solving for y̲n+1 gives y ̲ n + 1 = - A 0 - 1 A 1 y ̲ n + A 0 - 1 B 0 x ̲ n + 1 + A 0 - 1 B 1 x ̲ n y ̲ n + 1 = K y ̲ n + H 0 x ̲ n + 1 + H ˜ 1 x ̲ n which is a first order vector difference equation , . This is the fundamental block recursive algorithm that implements the original scalar difference equation in . It has several important characteristics. The block recursive formulation is similar to a state variable equation but the states are blocks or sections of the output , , , . The eigenvalues of K are the poles of the original scalar problem raised to the N power plus others that are zero. The longer the block length, the “more stable" the filter is, i.e. the further the poles are from the unit circle , , , , . If the block length were shorter than the denominator, the vector difference equation would be higher than first order. There would be a non zero A2. If the block length were shorter than the numerator, there would be a non zero B2 and a higher order block convolution operation. If the block length were one, the order of the vector equation would be the same as the scalar equation. They would be the same equation. The actual arithmetic that goes into the calculation of the output is partly recursive and partly convolution. The longer the block, the more the output is calculated by convolution and, the more arithmetic is required. It is possible to remove the zero eigenvalues in K by making K rectangular or square and N by N This results in a form even more similar to a state variable formulation , . This is briefly discussed below in section 2.3. There are several ways of using the FFT in the calculation of the various matrix products in and in and . Each has some arithmetic advantage for various forms and orders of the original equation. It is also possible to implement some of the operations using rectangular transforms, number theoretic transforms, distributed arithmetic, or other efficient convolution algorithms , , , , , . By choosing the block length equal to the period, a periodically time varying filter can be made block time invariant. In other words, all the time varying characteristics are moved to the finite matrix multiplies which leave the time invariant properties at the block level. This allows use of z-transform and other time-invariant methods to be used for stability analysis and frequency response analysis , . It also turns out to be related to filter banks and multi-rate filters , , .

Block State Formulation It is possible to reduce the size of the matrix operators in the block recursive description to give a form even more like a state variable equation , , . If K in has several zero eigenvalues, it should be possible to reduce the size of K until it has full rank. That was done in and the result is z ̲ n = K 1 z ̲ n - 1 + K 2 x ̲ n y ̲ n = H 1 z ̲ n - 1 + H 0 x ̲ n where H0 is the same N by N convolution matrix, N1 is a rectangular L by N partition of the convolution matrix H, K1 is a square N by N matrix of full rank, and K2 is a rectangular N by L matrix. This is now a minimal state equation whose input and output are blocks of the original input and output. Some of the matrix multiplications can be carried out using the FFT or other techniques.

Block Implementations of Digital Filters The advantage of the block convolution and recursion implementations is a possible improvement in arithmetic efficiency by using the FFT or other fast convolution methods for some of the multiplications in or , . There is the reduction of quantization effects due to an effective decrease in the magnitude of the eigenvalues and the possibility of easier parallel implementation for IIR filters. The disadvantages are a delay of at least one block length and an increased memory requirement. These methods could also be used in the various filtering methods for evaluating the DFT. This the chirp z-transform, Rader's method, and Goertzel's algorithm.

Multidimensional Formulation This process of partitioning the data vectors and the operator matrices can be continued by partitioning and and creating blocks of blocks to give a higher dimensional structure. One should use index mapping ideas rather than partitioned matrices for this approach , .

Periodically Time-Varying Discrete-Time Systems Most time-varying systems are periodically time-varying and this allows special results to be obtained. If the block length is set equal to the period of the time variations, the resulting block equations are time invariant and all to the time varying characteristics are contained in the matrix multiplications. This allows some of the tools of time invariant systems to be used on periodically time-varying systems. The PTV system is analyzed in , , , , the filter analysis and design problem, which includes the decimation–interpolation structure, is addressed in , , , and the bandwidth compression problem in . These structures can take the form of filter banks .

Multirate Filters, Filter Banks, and Wavelets Another area that is related to periodically time varying systems and to block processing is filter banks , . Recently the area of perfect reconstruction filter banks has been further developed and shown to be closely related to wavelet based signal analysis , , , . The filter bank structure has several forms with the polyphase and lattice being particularly interesting. An idea that has some elements of multirate filters, perfect reconstruction, and distributed arithmetic is given in , , . Parks has noted that design of multirate filters has some elements in common with complex approximation and of 2-D filter design , and is looking at using Tang's method for these designs.

Distributed Arithmetic Rather than grouping the individual scalar data values in a discrete-time signal into blocks, the scalar values can be partitioned into groups of bits. Because multiplication of integers, multiplication of polynomials, and discrete-time convolution are the same operations, the bit-level description of multiplication can be mixed with the convolution of the signal processing. The resulting structure is called distributed arithmetic , . It can be used to create an efficient table look-up scheme to implement an FIR or IIR filter using no multiplications by fetching previously calculated partial products which are stored in a table. Distributed arithmetic, block processing, and multi-dimensional formulations can be combined into an integrated powerful description to implement digital filters and processors. There may be a new form of distributed arithmetic using the ideas in , .

Direct Fast Convolution and Rectangular Transforms A relatively new approach uses index mapping directly to convert a one dimensional convolution into a multidimensional convolution , . This can be done by either a type-1 or type-2 map. The short convolutions along each dimension are then done by Winograd's optimal algorithms. Unlike for the case of the DFT, there is no savings of arithmetic from the index mapping alone. All the savings comes from efficient short algorithms. In the case of index mapping with convolution, the multiplications must be nested together in the center of the algorithm in the same way as for the WFTA. There is no equivalent to the PFA structure for convolution. The multidimensional convolution can not be calculated by row and column convolutions as the DFT was by row and column DFTs. It would first seem that applying the index mapping and optimal short algorithms directly to convolution would be more efficient than using DFTs and converting them to convolution to be calculated by the same optimal algorithms. In practical algorithms, however, the DFT method seems to be more efficient . A method that is attractive for special purpose hardware uses distributed arithmetic . This approach uses a table look up of precomputed partial products to produce a system that does convolution without requiring multiplications . Another method that requires special hardware uses number theoretic transforms , , to calculate convolution. These transforms are defined over finite fields or rings with arithmetic performed modulo special numbers. These transforms have rather limited flexibility, but when they can be used, they are very efficient.

Number Theoretic Transforms for Convolution

Results from Number Theory A basic review of the number theory useful for signal processing algorithms will be given here with specific emphasis on the congruence theory for number theoretic transforms , , , , .

Number Theoretic Transforms Here we look at the conditions placed on a general linear transform in order for it to support cyclic convolution. The form of a linear transformation of a length-N sequence of number is given by X ( k ) = ∑ n = 0 N - 1 t ( n , k ) x ( n ) for k=0,1,⋯,(N-1). The definition of cyclic convolution of two sequences is given by y ( n ) = ∑ m = 0 N - 1 x ( m ) h ( n - m ) for n=0,1,⋯,(N-1) and all indices evaluated modulo N. We would like to find the properties of the transformation such that it will support the cyclic convolution. This means that if X(k), H(k), and Y(k) are the transforms of x(n), h(n), and y(n) respectively, Y ( k ) = X ( k ) H ( k ) . The conditions are derived by taking the transform defined in of both sides of equation which gives Y ( k ) = ∑ n = 0 N - 1 t ( n , k ) ∑ m = 0 N - 1 x ( m ) h ( n - m ) = ∑ m = 0 N - 1 ∑ n = 0 N - 1 x ( m ) h ( n - m ) t ( n , k ) . Making the change of index variables, l=n-m, gives = ∑ m = 0 N - 1 ∑ l = 0 N - 1 x ( m ) h ( l ) t ( l + m , k ) . But from , this must be Y ( k ) = ∑ n = 0 N - 1 x ( n ) t ( n , k ) ∑ m = 0 N - 1 x ( m ) t ( m , k ) = ∑ m = 0 N - 1 ∑ l = 0 N - 1 x ( m ) h ( l ) t ( n , k ) t ( l , k ) . This must be true for all x(n), h(n), and k, therefore from and we have t ( m + l , k ) = t ( m , k ) t ( l , k ) For l=0 we have t ( m , k ) = t ( m , k ) t ( 0 , k ) and, therefore, t(0,k)=1. For l=m we have t ( 2 m , k ) = t ( m , k ) t ( m , k ) = t 2 ( m , k ) For l=pm we likewise have t ( p m , k ) = t p ( m , k ) and, therefore, t N ( m , k ) = t ( N m , k ) = t ( 0 , k ) = 1 . But t ( m , k ) = t m ( 1 , k ) = t k ( m , 1 ) , therefore, t ( m , k ) = t m k ( 1 , 1 ) . Defining t(1,1)=α gives the form for our general linear transform as X ( k ) = ∑ n = 0 N - 1 α n k x ( n ) where α is a root of order N , which means that N is the smallest integer such that αN=1. Theorem 1 The transform supports cyclic convolution if and only if α is a root of order N and N-1 is defined. This is discussed in , . Theorem 2 The transform supports cyclic convolution if and only if N | O ( M ) where O ( M ) = g c d { p 1 - 1 , p 2 - 1 , ⋯ , p l - 1 } and M = p 1 r 1 p 2 r 2 ⋯ p l r l . This theorem is a more useful form of Theorem 1. Notice that Nmax=O(M). One needs to find appropriate N, M, and α such that N should be appropriate for a fast algorithm and handle the desired sequence lengths. M should allow the desired dynamic range of the signals and should allow simple modular arithmetic. α should allow a simple multiplication for αnkx(n). We see that if M is even, it has a factor of 2 and, therefore, O(M)=Nmax=1 which implies M should be odd. If M is prime the O(M)=M-1 which is as large as could be expected in a field of M integers. For M=2k-1, let k be a composite k=pq where p is prime. Then 2p-1 divides 2pq-1 and the maximum possible length of the transform will be governed by the length possible for 2p-1. Therefore, only the prime k need be considered interesting. Numbers of this form are know as Mersenne numbers and have been used by Rader . For Mersenne number transforms, it can be shown that transforms of length at least 2p exist and the corresponding α=-2. Mersenne number transforms are not of as much interest because 2p is not highly composite and, therefore, we do not have FFT-type algorithms. For M=2k+1 and k odd, 3 divides 2k+1 and the maximum possible transform length is 2. Thus we consider only even k. Let k=s2t, where s is an odd integer. Then 22t divides 2s2t+1 and the length of the possible transform will be governed by the length possible for 22t+1. Therefore, integers of the form M=22t+1 are of interest. These numbers are known as Fermat numbers . Fermat numbers are prime for 0≤t≤4 and are composite for all t≥5. Since Fermat numbers up to F4 are prime, O(Ft)=2b where b=2t and we can have a Fermat number transform for any length N=2m where m≤b. For these Fermat primes the integer α=3 is of order N=2b allowing the largest possible transform length. The integer α=2 is of order N=2b=2t+1. This is particularly attractive since α to a power is multiplied times the data values in . The following table gives possible parameters for various Fermat number moduli. t b M = F t N 2 N 2 N m a x α for Nmax 3 8 2 8 + 1 16 32 256 3 4 16 2 16 + 1 32 64 65536 3 5 32 2 32 + 1 64 128 128 2 6 64 2 64 + 1 128 256 256 2

This table gives values of N for the two most important values of α which are 2 and 2. The second column give the approximate number of bits in the number representation. The third column gives the Fermat number modulus, the fourth is the maximum convolution length for α=2, the fifth is the maximum length for α=2, the sixth is the maximum length for any α, and the seventh is the α for that maximum length. Remember that the first two rows have a Fermat number modulus which is prime and second two rows have a composite Fermat number as modulus. Note the differences. The books, articles, and presentations that discuss NTT and related topics are , , , , , , , , , , , , . A recent book discusses NT in a signal processing context .