Skip to content Skip to navigation Skip to collection information

Connexions

You are here: Home » Content » Automatic Generation of Prime Length FFT Programs » Bilinear Forms for Circular Convolution

Navigation

Lenses

What is a lens?

Definition of a lens

Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

Who can create a lens?

Any individual member, a community, or a respected organization.

What are tags? tag icon

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

This content is ...

Affiliated with (What does "Affiliated with" mean?)

This content is either by members of the organizations listed or about topics related to the organizations listed. Click each link to see a list of all content affiliated with the organization.
  • Rice Digital Scholarship

    This collection is included in aLens by: Digital Scholarship at Rice University

    Click the "Rice Digital Scholarship" link to see all content affiliated with them.

  • Featured Content display tagshide tags

    This collection is included inLens: Connexions Featured Content
    By: Connexions

    Click the "Featured Content" link to see all content affiliated with them.

    Click the tag icon tag icon to display tags associated with this content.

Also in these lenses

  • UniqU content

    This collection is included inLens: UniqU's lens
    By: UniqU, LLC

    Click the "UniqU content" link to see all content selected in this lens.

  • Lens for Engineering

    This module and collection are included inLens: Lens for Engineering
    By: Sidney Burrus

    Click the "Lens for Engineering" link to see all content selected in this lens.

Recently Viewed

This feature requires Javascript to be enabled.

Tags

(What is a tag?)

These tags come from the endorsement, affiliation, and other lenses that include this content.
 

Bilinear Forms for Circular Convolution

Module by: Ivan Selesnick, C. Sidney Burrus. E-mail the authors

Summary: This collection of modules is from a Rice University, ECE Department Technical Report written around September 1994. It grew out of the doctoral and post doctoral research of Ivan Selesnick working with Prof. C. Sidney Burrus at Rice. Earlier reports on this work were published in the ICASSP and ISCAS conference proceedings in 1992-94 and a fairly complete report was published in the IEEE Transaction on Signal Processing in January 1996.

Bilinear Forms for Circular Convolution

A basic technique in fast algorithms for convolution is that of interpolation. That is, two polynomials are evaluated at some common points and these values are multiplied [1], [3], [4]. By interpolating these products, the product of the two original polynomials can be determined. In the Winograd short convolution algorithms, this technique is used and the common points of evaluation are the simple integers, 0, 1, and -1-1. Indeed, the computational savings of the interpolation technique depends on the use of special points at which to interpolate. In the Winograd algorithm the computational savings come from the simplicity of the small integers. (As an algorithm for convolution, the FFT interpolates over the roots of unity.) This interpolation method is often called the Toom-Cook method and it is given by two matrices that describe a bilinear form.

We use bilinear forms to give a matrix formulation of the split nesting algorithm. The split nesting algorithm combines smaller convolution algorithms to obtain algorithms for longer lengths. We use the Kronecker product to explicitly describe the way in which smaller convolution algorithms are appropriately combined.

The Scalar Toom-Cook Method

First we consider the linear convolution of two nn point sequences. Recall that the linear convolution of hh and xx can be represented by a matrix vector product. When n=3n=3:

h 0 h 1 h 0 h 2 h 1 h 0 h 2 h 1 h 2 x 0 x 1 x 2 h 0 h 1 h 0 h 2 h 1 h 0 h 2 h 1 h 2 x 0 x 1 x 2
(1)

This linear convolution matrix can be written as h0H0+h1H1+h2H2h0H0+h1H1+h2H2 where HkHk are clear.

The product ∑k=0n-1hkHkx∑k=0n-1hkHkx can be found using the Toom-Cook algorithm, an interpolation method. Choose 2n-12n-1 interpolation points, i1,⋯,i2n-1i1,⋯,i2n-1, and let AA and CC be matrices given by

A = i 1 0 ⋯ i 1 n - 1 ⋮ i 2 n - 1 0 ⋯ i 2 n - 1 n - 1 and C = i 1 0 ⋯ i 1 2 n - 2 ⋮ i 2 n - 1 0 ⋯ i 2 n - 1 2 n - 2 - 1 . A = i 1 0 ⋯ i 1 n - 1 ⋮ i 2 n - 1 0 ⋯ i 2 n - 1 n - 1 and C = i 1 0 ⋯ i 1 2 n - 2 ⋮ i 2 n - 1 0 ⋯ i 2 n - 1 2 n - 2 - 1 .
(2)

That is, AA is a degree n-1n-1 Vandermonde matrix and CC is the inverse of the degree 2n-22n-2 Vandermonde matrix specified by the same points specifying AA. With these matrices, one has

∑ k = 0 n - 1 h k H k x = C A h * A x ∑ k = 0 n - 1 h k H k x = C A h * A x
(3)

where ** denotes point by point multiplication. The terms AhAh and AxAx are the values of H(s)H(s) and X(s)X(s) at the points i1,⋯i2n-1i1,⋯i2n-1. The point by point multiplication gives the values Y(i1),⋯,Y(i2n-1)Y(i1),⋯,Y(i2n-1). The operation of CC obtains the coefficients of Y(s)Y(s) from its values at these points of evaluation. This is the bilinear form and it implies that

H k = Cdiag ( A e k ) A . H k = Cdiag ( A e k ) A .
(4)

Example 1

If n=2n=2, then equations Equation 3 and Equation 4 give

h 0 0 h 1 h 0 0 h 1 x = C A h * A x h 0 0 h 1 h 0 0 h 1 x = C A h * A x
(5)

When the interpolation points are 0, 1,and -1-1,

A = 1 0 1 1 1 - 1 and C = 1 0 0 0 . 5 - . 5 - 1 . 5 . 5 A = 1 0 1 1 1 - 1 and C = 1 0 0 0 . 5 - . 5 - 1 . 5 . 5
(6)

However, AA and CC do not need to be Vandermonde matrices as in Equation 2. For example, see the two point linear convolution algorithm in the appendix. As long as AA and CC are matrices such that Hk=Cdiag(Aek)AHk=Cdiag(Aek)A, then the linear convolution of xx and hh is given by the bilinear form y=C{Ah*Ax}y=C{Ah*Ax}. More generally, as long as AA, BB and CC are matrices satisfying Hk=Cdiag(Bek)AHk=Cdiag(Bek)A, then y=C{Bh*Ax}y=C{Bh*Ax} computes the linear convolution of hh and xx. For convenience, if C{Ah*Ax}C{Ah*Ax} computes the nn point linear convolution of hh and xx (both hh and xx are nn point sequences), then we say “(A,B,C)(A,B,C) describes a bilinear form for nn point linear convolution."

Similarly, we can write a bilinear form for cyclotomic convolution. Let dd be any positive integer and let X(s)X(s) and H(s)H(s) be polynomials of degree φ(d)-1φ(d)-1 where φ(·)φ(·) is the Euler totient function. If AA, BB and CC are matrices satisfying CΦdk=Cdiag(Bek)ACΦdk=Cdiag(Bek)A for 0≤k≤φ(d)-10≤k≤φ(d)-1, then the coefficients of Y(s)=〈X(s)H(s)〉Φd(s)Y(s)=〈X(s)H(s)〉Φd(s) are given by y=C{Bh*Ax}y=C{Bh*Ax}. As above, if y=C{Bh*Ax}y=C{Bh*Ax} computes the dd-cyclotomic convolution, then we say “(A,B,C)(A,B,C) describes a bilinear form for Φd(s)Φd(s) convolution."

But since 〈X(s)H(s)〉Φd(s)〈X(s)H(s)〉Φd(s) can be found by computing the product of X(s)X(s) and H(s)H(s) and reducing the result, a cyclotomic convolution algorithm can always be derived by following a linear convolution algorithm by the appropriate reduction operation: If GG is the appropriate reduction matrix and if (A,B,F)(A,B,F) describes a bilinear form for a φ(d)φ(d) point linear convolution, then (A,B,GF)(A,B,GF) describes a bilinear form for Φd(s)Φd(s) convolution. That is, y=GF{Bh*Ax}y=GF{Bh*Ax} computes the coefficients of 〈X(s)H(s)〉Φd(s)〈X(s)H(s)〉Φd(s).

Circular Convolution

By using the Chinese Remainder Theorem for polynomials, circular convolution can be decomposed into disjoint cyclotomic convolutions. Let pp be a prime and consider pp point circular convolution. Above we found that

S p = R p - 1 1 C Φ p R p S p = R p - 1 1 C Φ p R p
(7)

and therefore

S p k = R p - 1 1 C Φ p k R p . S p k = R p - 1 1 C Φ p k R p .
(8)

If (Ap,Bp,Cp)(Ap,Bp,Cp) describes a bilinear form for Φp(s)Φp(s) convolution, then

S p k = R p - 1 1 C p diag 1 B p R p e k 1 A p R p S p k = R p - 1 1 C p diag 1 B p R p e k 1 A p R p
(9)

and consequently the circular convolution of hh and xx can be computed by

y = R p - 1 C { B R p h * A R p x } y = R p - 1 C { B R p h * A R p x }
(10)

where A=1⊕ApA=1⊕Ap, B=1⊕BpB=1⊕Bp and C=1⊕CpC=1⊕Cp. We say (A,B,C)(A,B,C) describes a bilinear form for pp point circular convolution. Note that if (D,E,F)(D,E,F) describes a (p-1)(p-1) point linear convolution then ApAp, BpBp and CpCp can be taken to be Ap=DAp=D, Bp=EBp=E and Cp=GpFCp=GpF where GpGp represents the appropriate reduction operations. Specifically, GpGp is given by Equation 42 from Preliminaries.

Next we consider pepe point circular convolution. Recall that Spe=Rpe-1⊕i=0eCΦpiRpeSpe=Rpe-1⊕i=0eCΦpiRpe as in Equation 27 from Preliminaries so that the circular convolution is decomposed into a set of e+1e+1 disjoint Φpi(s)Φpi(s) convolutions. If (Api,Bpi,Cpi)(Api,Bpi,Cpi) describes a bilinear form for Φpi(s)Φpi(s) convolution and if

A = 1 ⊕ A p ⊕ ⋯ ⊕ A p e B = 1 ⊕ B p ⊕ ⋯ ⊕ B p e C = 1 ⊕ C p ⊕ ⋯ ⊕ C p e A = 1 ⊕ A p ⊕ ⋯ ⊕ A p e B = 1 ⊕ B p ⊕ ⋯ ⊕ B p e C = 1 ⊕ C p ⊕ ⋯ ⊕ C p e
(11)

then (ARpe,BRpe,Rpe-1C)(ARpe,BRpe,Rpe-1C) describes a bilinear form for pepe point circular convolution. In particular, if (Dd,Ed,Fd)(Dd,Ed,Fd) describes a bilinear form for dd point linear convolution, then ApiApi, BpiBpi and CpiCpi can be taken to be

A p i = D φ ( p i ) B p i = E φ ( p i ) C p i = G p i F φ ( p i ) A p i = D φ ( p i ) B p i = E φ ( p i ) C p i = G p i F φ ( p i )
(12)

where GpiGpi represents the appropriate reduction operation and φ(·)φ(·) is the Euler totient function. Specifically, GpiGpi has the following form

G p i = I ( p - 1 ) p i - 1 - 1 -p - 1 ⊗ I p i - 1 I ( p - 2 ) p i - 1 - 1 0 p i - 1 + 1 , ( p - 2 ) p i - 1 - 1 G p i = I ( p - 1 ) p i - 1 - 1 -p - 1 ⊗ I p i - 1 I ( p - 2 ) p i - 1 - 1 0 p i - 1 + 1 , ( p - 2 ) p i - 1 - 1
(13)

if p≥3p≥3, while

G 2 i = I 2 i - 1 - I 2 i - 1 - 1 0 1 , 2 i - 1 - 1 . G 2 i = I 2 i - 1 - I 2 i - 1 - 1 0 1 , 2 i - 1 - 1 .
(14)

Note that the matrix RpeRpe block diagonalizes SpeSpe and each diagonal block represents a cyclotomic convolution. Correspondingly, the matrices AA, BB and CC of the bilinear form also have a block diagonal structure.

The Split Nesting Algorithm

We now describe the split-nesting algorithm for general length circular convolution [4]. Let n=p1e1⋯pkekn=p1e1⋯pkek where pipi are distinct primes. We have seen that

S n = P t R - 1 ⊕ d | n Ψ ( d ) R P S n = P t R - 1 ⊕ d | n Ψ ( d ) R P
(15)

where PP is the prime factor permutation P=Pp1e1,⋯,pkekP=Pp1e1,⋯,pkek and RR represents the reduction operations. For example, see Equation 46 in Preliminaries. RPRP block diagonalizes SnSn and each diagonal block represents a multi-dimensional cyclotomic convolution. To obtain a bilinear form for a multi-dimensional convolution, we can combine bilinear forms for one-dimensional convolutions. If (Apji,Bpji,Cpji)(Apji,Bpji,Cpji) describes a bilinear form for Φpji(s)Φpji(s) convolution and if

A = ⊕ d | n A d B = ⊕ d | n B d C = ⊕ d | n C d A = ⊕ d | n A d B = ⊕ d | n B d C = ⊕ d | n C d
(16)

with

A d = ⊗ p | d , p ∈ P A H d ( p ) B d = ⊗ p | d , p ∈ P B H d ( p ) C d = ⊗ p | d , p ∈ P C H d ( p ) A d = ⊗ p | d , p ∈ P A H d ( p ) B d = ⊗ p | d , p ∈ P B H d ( p ) C d = ⊗ p | d , p ∈ P C H d ( p )
(17)

where Hd(p)Hd(p) is the highest power of pp dividing dd, and PP is the set of primes, then (ARP,BRP,PtR-1C)(ARP,BRP,PtR-1C) describes a bilinear form for nn point circular convolution. That is

y = P t R - 1 C B R P h * A R P x y = P t R - 1 C B R P h * A R P x
(18)

computes the circular convolution of hh and xx.

As above (Apji,Bpji,Cpji)(Apji,Bpji,Cpji) can be taken to be (Dφ(pji),Eφ(pji),GpjiFφ(pji))(Dφ(pji),Eφ(pji),GpjiFφ(pji)) where (Dd,Ed,Fd)(Dd,Ed,Fd) describes a bilinear form for dd point linear convolution. This is one particular choice for (Apji,Bpji,Cpji)(Apji,Bpji,Cpji) - other bilinear forms for cyclotomic convolution that are not derived from linear convolution algorithms exist.

Example 2

A 45 point circular convolution algorithm:

y = P t R - 1 C B R P h * A R P x y = P t R - 1 C B R P h * A R P x
(19)

where

P = P 9 , 5 R = R 9 , 5 A = 1 ⊕ A 3 ⊕ A 9 ⊕ A 5 ⊕ ( A 3 ⊗ A 5 ) ⊕ ( A 9 ⊗ A 5 ) B = 1 ⊕ B 3 ⊕ B 9 ⊕ B 5 ⊕ ( B 3 ⊗ B 5 ) ⊕ ( B 9 ⊗ B 5 ) C = 1 ⊕ C 3 ⊕ C 9 ⊕ C 5 ⊕ ( C 3 ⊗ C 5 ) ⊕ ( C 9 ⊗ C 5 ) P = P 9 , 5 R = R 9 , 5 A = 1 ⊕ A 3 ⊕ A 9 ⊕ A 5 ⊕ ( A 3 ⊗ A 5 ) ⊕ ( A 9 ⊗ A 5 ) B = 1 ⊕ B 3 ⊕ B 9 ⊕ B 5 ⊕ ( B 3 ⊗ B 5 ) ⊕ ( B 9 ⊗ B 5 ) C = 1 ⊕ C 3 ⊕ C 9 ⊕ C 5 ⊕ ( C 3 ⊗ C 5 ) ⊕ ( C 9 ⊗ C 5 )
(20)

and where (Apji,Bpji,Cpji)(Apji,Bpji,Cpji) describes a bilinear form for Φpji(s)Φpji(s) convolution.

The Matrix Exchange Property

The matrix exchange property is a useful technique that, under certain circumstances, allows one to save computation in carrying out the action of bilinear forms [2]. Suppose

y = C A x * B h y = C A x * B h
(21)

as in Equation 18. When hh is known and fixed, BhBh can be pre-computed so that yy can be found using only the operations represented by CC and AA and the point by point multiplications denoted by **. The operation of BB is absorbed into the multiplicative constants. Note that in Equation 18, the matrix corresponding to CC is more complicated than is BB. It is therefore advantageous to absorb the work of CC instead of BB into the multiplicative constants if possible. This can be done when yy is the circular convolution of xx and hh by using the matrix exchange property.

To explain the matrix exchange property we draw from [2]. Note that y=Cdiag(Ax)Bhy=Cdiag(Ax)Bh, so that Cdiag(Ax)BCdiag(Ax)B must be the corresponding circulant matrix,

Cdiag ( A x ) B = x 0 x n - 1 ⋯ x 1 x 1 x 0 x 2 ⋮ x n - 1 x n - 2 x 0 . Cdiag ( A x ) B = x 0 x n - 1 ⋯ x 1 x 1 x 0 x 2 ⋮ x n - 1 x n - 2 x 0 .
(22)

Since Cdiag(Ax)B=JCdiag(Ax)BtJCdiag(Ax)B=JCdiag(Ax)BtJ where JJ is the reversal matrix, one gets

y = C A x * B h = Cdiag ( A x ) B h = J Cdiag ( A x ) B t J h = J B t d i a g ( A x ) C t J h = J B t A x * C t J h y = C A x * B h = Cdiag ( A x ) B h = J Cdiag ( A x ) B t J h = J B t d i a g ( A x ) C t J h = J B t A x * C t J h
(23)

As noted in [2], the matrix exchange property can be used whenever y=T(x)hy=T(x)h where T(x)T(x) satisfies T(x)=J1T(x)tJ2T(x)=J1T(x)tJ2 for some matrices J1J1 and J2J2. In that case one gets y=J1BtAx*CtJ2hy=J1BtAx*CtJ2h.

Applying the matrix exchange property to Equation 18 one gets

y = J P t R t B t C t R - t P J h * A R P x . y = J P t R t B t C t R - t P J h * A R P x .
(24)

Example 3

A 45 point circular convolution algorithm:

y = J P t R t B t u * A R P x y = J P t R t B t u * A R P x
(25)

where u=CtR-tPJhu=CtR-tPJh and

P = P 9 , 5 R = R 9 , 5 A = 1 ⊕ A 3 ⊕ A 9 ⊕ A 5 ⊕ ( A 3 ⊗ A 5 ) ⊕ ( A 9 ⊗ A 5 ) B t = 1 ⊕ B 3 t ⊕ B 9 t ⊕ B 5 t ⊕ ( B 3 t ⊗ B 5 t ) ⊕ ( B 9 t ⊗ B 5 t ) C t = 1 ⊕ C 3 t ⊕ C 9 t ⊕ C 5 t ⊕ ( C 3 t ⊗ C 5 t ) ⊕ ( C 9 t ⊗ C 5 t ) P = P 9 , 5 R = R 9 , 5 A = 1 ⊕ A 3 ⊕ A 9 ⊕ A 5 ⊕ ( A 3 ⊗ A 5 ) ⊕ ( A 9 ⊗ A 5 ) B t = 1 ⊕ B 3 t ⊕ B 9 t ⊕ B 5 t ⊕ ( B 3 t ⊗ B 5 t ) ⊕ ( B 9 t ⊗ B 5 t ) C t = 1 ⊕ C 3 t ⊕ C 9 t ⊕ C 5 t ⊕ ( C 3 t ⊗ C 5 t ) ⊕ ( C 9 t ⊗ C 5 t )
(26)

and where (Apji,Bpji,Cpji)(Apji,Bpji,Cpji) describes a bilinear form for Φpji(s)Φpji(s) convolution.

References

  1. Blahut, R. E. (1985). Fast Algorithms for Digital Signal Processing. Addison-Wesley.
  2. Johnson, H. W. and Burrus, C. S. (1985, February). On the Structure of Efficient DFT Algorithms. IEEE Trans. Acoust., Speech, Signal Proc., 33(1), 248-254.
  3. Myers, D. G. (1990). Digital Signal Processing: Efficient Convolution and Fourier Transform Techniques. Prentice Hall.
  4. Nussbaumer, H. J. (1982). Fast Fourier Transform and Convolution Algorithms. Sringer-Verlag.

Collection Navigation

Content actions

Download module as:

PDF | More downloads ...

Add:

Collection to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

Definition of a lens

Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

Who can create a lens?

Any individual member, a community, or a respected organization.

What are tags? tag icon

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks

Module to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

Definition of a lens

Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

Who can create a lens?

Any individual member, a community, or a respected organization.

What are tags? tag icon

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks