An Illustrated Introduction to the Truncated Fourier Transform

1 Introduction

In 1965 Cooley and Tukey developed the Discrete Fourier Transform (dft) to recover continuous functions from discrete samples. This was a landmark discovery because it allowed for the digital manipulation of analogue signals (like sound) by computers. Soon after a variant called the Fast Fourier Transform (fft) eclipsed the dft to the extent that fft is often mistakenly substituted for dft. As its name implies, the fft is a method for computing the dft faster.

ffts have an interesting application in Computer Algebra. Let $R$ be a ring with $2 \in R$ a unit. If $R$ has a primitive $n$ th root of unity $ω$ with $n = 2^{p}$ (i.e. $ω^{n / 2} = - 1$ ) then the fft computes the product of two polynomials $P, Q \in R [x]$ with $deg (P Q) < n$ in $O (n log n)$ operations in $R$ . Unfortunately, when $deg (P Q)$ is sufficiently far from a power of two many computations wasted. This deficiency was addressed by the signal processing community using a method called fft-pruning [3]. However, the difficult inversion of this method is due to van der Hoeven [4][5].

In Section 2 we outline the dft, including a method for its non-recursive implementation. In Section 3 we develop the “pruned” variant, called Truncated Fourier Transform (tft). Finally, in Section 4, we show how the tft can be inverted and outline the algorithm for doing so in Section 5.

2 The Discrete Fourier Transform

For this paper let $R$ be a ring with $2 \in R$ a unit and $ω \in R$ an $n$ th root of unity. The Discrete Fourier Transform^*^**In signal processing community this is called the “decimation-in-time” variant of the fft., with respect to $ω$ , of vector $a = (a_{0}, \dots, a_{n - 1}) \in R^{n}$ is the vector $^a = ({^a}_{0}, \dots, {^a}_{n - 1}) \in R^{n}$ with

{^a}_{i} = n - 1 \sum j = 0 a_{j} ω^{i j} .

Alternatively we can view these $n$ -tuples as encoding the coefficients of polynomials from $R [x]$ and define the dft with respect to $ω$ as the mapping

	${DFT}_{ω} : R [x]$	$\to R^{n}$
	$A (x) = a_{0} + \dots + a_{n - 1} x^{n - 1}$	$\mapsto (A (ω^{0}), \dots, A (ω^{n - 1})) .$

We let the relationship between $A$ and its coefficients be implicit and write

{DFT}_{ω} (a_{0}, \dots, a_{n - 1}) := (A (ω^{0}), \dots, A (ω^{n - 1}))

when $A = a_{0} + \dots + a_{n - 1} x^{n - 1}$ .

The dft can be computed efficiently using binary splitting. This method requires evaluation only at $ω^{2^{i}}$ for $i \in {0, \dots, p - 1}$ , rather than at all $ω^{0}, \dots, ω^{n - 1}$ . To compute the dft of $A$ with respect to $ω$ we write

(b_{0}, c_{0}, \dots, b_{n / 2 - 1}, c_{n / 2 - 1}) := (a_{0}, \dots, a_{n - 1})

and recursively compute the dft of $(b_{0}, \dots, b_{n / 2 - 1})$ and $(c_{0}, \dots, c_{n / 2 - 1})$ with respect to $ω^{2}$ :

	$({^b}_{0}, \dots, {^b}_{n / 2 - 1}) := {DFT}_{ω^{2}} (b_{0}, \dots, b_{n / 2 - 1}),$
	$({^c}_{0}, \dots, {^c}_{n / 2 - 1}) := {DFT}_{ω^{2}} (c_{0}, \dots, c_{n / 2 - 1}) .$

Finally, we construct $^a$ according to

	${DFT}_{ω} (a_{0}, \dots, a_{n - 1}) = ($	${^b}_{0} + {^c}_{0}, \dots, {^b}_{n / 2 - 1} + {^c}_{n / 2 - 1} ω^{n / 2 - 1},$
		${^b}_{0} - {^c}_{0}, \dots, {^b}_{n / 2 - 1} - {^c}_{n / 2 - 1} ω^{n / 2 - 1}) .$

This description has a natural implementation as a recursive algorithm, but in practice it is often more efficient to implement an in-place algorithm that eliminates the overhead of creating recursive stacks.

Definition.

Let $i$ and $p$ be a positive integers and let $i = i_{0} 2^{0} + \dots + i_{p} 2^{p}$ for $i_{0}, \dots, i_{p} \in {0, 1}$ . The length- $p$ bitwise reverse of $i$ is given by

[i]_{p} := i_{p} 2^{0} + \dots + i_{0} 2^{p} .

Example.

$[3]_{5} = 24$ and $[11]_{5} = 26$ because

	$[3]_{5}$	$= [1 \cdot 2^{0} + 1 \cdot 2^{1} + 0 \cdot 2^{2} + 0 \cdot 2^{3} + 0 \cdot 2^{4}]_{5}$
		$= 0 \cdot 2^{0} + 0 \cdot 2^{1} + 0 \cdot 2^{2} + 1 \cdot 2^{3} + 1 \cdot 2^{4}$
		$= 24$

and

	$[11]_{5}$	$= [1 \cdot 2^{0} + 1 \cdot 2^{1} + 0 \cdot 2^{2} + 1 \cdot 2^{3} + 0 \cdot 2^{4}]_{5}$
		$= 0 \cdot 2^{0} + 1 \cdot 2^{1} + 0 \cdot 2^{2} + 1 \cdot 2^{3} + 1 \cdot 2^{4}$
		$= 26.$

Notice if we were to write 3, 24, 11, and 26 as a binary numbers to five digits we have 00011 reverses to 11000 and 01011 reverses to 11010 — in fact this is the inspiration for the name “bitwise reverse.”

For the in-place non-recursive dft algorithm, we require only one vector of length $n$ . Initially, at step zero, this vector is

x_{0} = (x_{0, 0}, \dots, x_{0, n - 1}) := (a_{0}, \dots, a_{n - 1})

and is updated (incrementally) at steps $s \in {1, \dots, p}$ by the rule

[\begin{matrix} x_{s, i m_{s} + j} x_{s, (i + 1) m_{s} + j} \end{matrix}] := [\begin{matrix} 1 & ω^{[i]_{s} m_{s}} 1 & - ω^{[i]_{s} m_{s}} \end{matrix}] [\begin{matrix} x_{s - 1, i m_{s} + j} x_{s - 1, (i + 1) m_{s} + j} \end{matrix}]

(1)

where $m_{s} = 2^{p - s}$ and for all $i \in {0, 2, \dots, n / m_{s} - 2}$ , $j \in {0, \dots, m_{s} - 1}$ . Note that two additions and one multiplication are done in (1) as one product is merely the negation of the other.

We illustrate the dependencies of the $x_{s, i}$ values in (1) with

We call this a “butterfly” after the shape it forms and may say $m_{s}$ controls the width — the value of which decreases as $s$ increases. By placing these butterflies on a $s \times n$ grid (Figure 1) we can see which values of $x_{s}$ are required to compute particular entires of $x_{s + 1}$ (and vice-versa). For example

denotes that $x_{3, 9}$ and $x_{3, 11}$ are required to determine $x_{2, 9}$ and $x_{2, 11}$ (and vice-versa).

Figure 1: Schematic representation of Equation (1) using butterflies to illustrate the value dependencies at various steps $s$ . The grid has rows $s = 0, \dots, 4$ and columns $n = 0, \dots, 15$ for illustrating the $x_{s, i}$ values.

Using induction over $s$ ,

x_{s, i m_{s} + j} = {({DFT}_{ω^{m_{s}}} (a_{j}, a_{m_{s} + j}, \dots, a_{n - m_{s} + j}))}_{[i]_{s}},

for all $i \in {0, \dots, n / m_{s} - 1}$ and $j \in {0, \dots, m_{s} - 1}$ [4]. In particular, when $s = p$ and $j = 0$ , we have

x_{p, i} = {^a}_{[i]_{p}} and {^a}_{i} = x_{p, [i]_{p}}

for all $i \in {0, \dots, n - 1}$ . That is, $^a$ is a (specific) permutation of $x_{p}$ as illustrated in Figure 2.

The key property of the dft is that it is straightforward to invert, that is to recover $a$ from $^a$ :

{DFT}_{ω^{- 1}} (^a)_{i} = {DFT}_{ω^{- 1}} ({% DFT}_{ω} (a))_{i} = n - 1 \sum k = 0 n - 1 \sum j = 0 a_{i} ω^{(i - k) j} = n a_{i}

(2)

since $\sum_{j = 0}^{n - 1} ω^{(i - k) j} = 0$ whenever $i \neq k$ . This yields a polynomial multiplication algorithm that does $O (n log n)$ operations in $R$ (see [1, §4.7] for the outline of this algorithm).

Figure 2: The Discrete Fourier Transform for $n = 16$ . The top row, corresponding to $s = 0$ , are the initial values $a$ . The bottom row, corresponding to $s = 4$ , is a permutation of $^a$ (the result of the DFT on $a$ ).

3 The Truncated Fourier Transform

The motivation behind the Truncated Fourier Transform (tft) is the observation that many computations are wasted when the length of $a$ (the input) is not a power of two.^†^††The tft is exactly equivalent to a technique called “fft pruning” in the signal processing literature [3]. This is entirely the fault of the strategy where one “completes” the $ℓ$ -tuple $a = (a_{0}, \dots, a_{ℓ - 1})$ by setting $a_{i} = 0$ when $i \geq ℓ$ to artificially extend the length of $a$ to the nearest power of two (so the dft can be executed as usual).

However, despite the fact that we may only want $ℓ$ components of $^a$ , the dft will calculate all of them. Thus computation is wasted. We illustrate this in Figures 3 and 4. This type of wasted computation is relevant when using the dft to multiply polynomials — their products are rarely of degree one less some power of two.

The definition of the tft is similar to that of the dft with the exception that the input and output vector ( $a$ resp. $^a$ ) are not necessarily of length some power of two. More precisely the tft of an $ℓ$ -tuple $(a_{0}, \dots, a_{ℓ - 1}) \in R^{ℓ}$ is the $ℓ$ -tuple

(A (ω^{[0]_{p}}), \dots, A (ω^{[ℓ - 1]_{p}})) \in R^{ℓ} .

where $n = 2^{p}$ , $ℓ < n$ (usually $ℓ \geq n / 2$ ) and $ω$ a $n$ th root of unity.

Figure 3: The dft with “artificial” zero points (large black dots).

Figure 4: Removing all unnecessary computations from Figure 3 gives the schematic representation of the tft.

Remark.

van der Hoeven [4] gives a more general description of the tft where one can choose an initial vector $(x_{0, i_{0}}, \dots, x_{0, i_{n}})$ and target vector $(x_{p, j_{0}}, \dots, x_{p, j_{n}})$ . Provided the $i_{k}$ ’s are distinct one can carry out the tft by considering the full dft and removing all computations not required for the desired output. In this paper, we restrict our discussion to that of the scenario in Figure 4 (where the input and output are the same initial segments) because it can be used for polynomial multiplication, and because it yields the most improvement. $⋄$

If we only allow ourselves to operate in a size $n$ vector it is straightforward to modify the in-place dft algorithm from the previous section to execute the tft. (It should be emphasized that this only saves computation and not space. For a “true” in-place tft algorithm that operates in an array of size $ℓ$ , see Harvey and Roche’s [2].) At stage $s$ it suffices to compute $(x_{s, 0}, \dots, x_{s, j})$ with $j = ⌈ ℓ / m_{s} ⌉ m_{s} - 1$ where $m_{s} = 2^{p - s}$ .^‡^‡‡This is a correction to the bound given in [5] as pointed out in [4].

Theorem 1.

Let $n = 2^{p}$ , $1 \leq ℓ < n$ and $ω \in R$ be a primitive $n$ th root of unity in $R$ . The tft of an $ℓ$ -tuple $(a_{0}, \dots, a_{ℓ - 1})$ with respect to $ω$ can be computed using at most $ℓ p + n$ additions and $⌊ (ℓ p + n) / 2 ⌋$ multiplications of powers of $ω$ .

Proof.

Let $j = (⌈ ℓ_{s} / m_{s} ⌉) m_{s} - 1$ ; at stage $s$ we compute $(x_{s, 0}, \dots, x_{s, j})$ . So, in addition to $x_{s, 0}, \dots, x_{s, ℓ - 1}$ we compute

(⌈ ℓ / m_{s} ⌉) m_{s} - 1 - ℓ \leq m_{s}

more values. Therefore, in total, we compute at most

p ℓ + p \sum s = 1 m_{s}

= p ℓ + 2^{p - 1} + 2^{p - 2} + \dots + 1 = p ℓ + 2^{p} - 1 < p ℓ + n

values $x_{s, i}$ . The result follows. ∎

4 Inverting The Truncated Fourier Transform

Unfortunately, the tft cannot be inverted by merely doing another tft with $ω^{- 1}$ and adjusting the output by some constant factor like inverse of the dft. Simply put: we are missing information and must account for this.

Example.

Let $R = Z / 13 Z$ , $n = 2^{2} = 4$ , with $ω = 5$ a $n$ th primitive root of unity. Setting $A (x) = a_{0} + a_{1} x + a_{2} x^{2}$ , the tft of $a = (a_{0}, a_{1}, a_{2})$ at 5 is

⎡ ⎢ ⎢ ⎣ \begin{matrix} A (ω^{0}) A (ω^{2}) A (ω^{1}) \end{matrix} ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎣ \begin{matrix} A (1) A (- 1) A (5) \end{matrix} ⎤ ⎥ ⎦ = ⎡ ⎢ ⎣ \begin{matrix} a_{0} + a_{1} + a_{2} a_{0} - a_{1} + a_{2} a_{0} + 5 a_{1} - a_{2} \end{matrix} ⎤ ⎥ ⎦ .

Now, to show the tft of this with respect to $ω^{- 1}$ is not $a$ , define

b = ⎡ ⎢ ⎣ \begin{matrix} b_{0} b_{1} b_{2} \end{matrix} ⎤ ⎥ ⎦ = ⎡ ⎢ ⎣ \begin{matrix} a_{0} + a_{1} + a_{2} a_{0} - a_{1} + a_{2} a_{0} + 5 a_{1} - a_{2} \end{matrix} ⎤ ⎥ ⎦ .

The tft of $b$ with respect to $ω^{- 1} = - 5$ is

⎡ ⎢ ⎢ ⎣ \begin{matrix} B (ω^{0}) B (ω^{- 2}) B (ω^{- 1}) \end{matrix} ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎣ \begin{matrix} B (1) B (- 1) B (5) \end{matrix} ⎤ ⎥ ⎦ = ⎡ ⎢ ⎣ \begin{matrix} b_{0} + b_{1} + b_{2} b_{0} - b_{1} + b_{2} b_{0} - 5 b_{1} - b_{2} \end{matrix} ⎤ ⎥ ⎦ = ⎡ ⎢ ⎣ \begin{matrix} 3 a_{0} + 5 a_{1} + a_{2} a_{0} - 6 a_{1} - a_{2} - 5 a_{0} + a_{1} - 3 a_{2} \end{matrix} ⎤ ⎥ ⎦

which is not a constant multiple of ${TFT}_{ω} (a)$ .

This discrepancy is caused by the completion of $b$ to $(b_{0}, b_{1}, b_{2}, 0)$ — we should have instead completed $b$ to $(b_{0}, b_{1}, b_{2}, A (- 5))$ .

To invert the tft we use the fact that whenever two values among

x_{s, i m_{s} + j}, x_{s - 1, i m_{s} + j} and x_{s, (i + 1) m_{s} + j}, x_{s - 1, (i + 1) m_{s} + j}

are known, that the other values can be deduced. That is, if two values of some butterfly are known then the other two values can be calculated using (1) as the relevant matrix is invertible. Moreover, these relations only involve shifting (multiplication and division by two), additions, subtractions and multiplications by roots of unity — an ideal scenario for implementation.

As with dft, observe that $x_{p - k, 0}, \dots, x_{p - k, 2^{k} - 1}$ can be calculated from $x_{p, 0}, \dots, x_{p, 2^{k} - 1}$ . This is because all the butterfly relations necessary to move up like this never require $x_{s, 2^{k} + j}$ for any $s \in {p - k, \dots, p}$ and $j > 0$ . This is illustrated in Figure 5. More generally, we have that

x_{p, 2^{j} + 2^{k}}, \dots, x_{p, 2^{j} + 2^{k} - 1}

is sufficient information to compute

x_{p - k, 2^{j}}, \dots, x_{p - k, 2^{j} + 2^{k} - 1}

provided that $0 < k \leq j < p$ . (In Algorithm 1, that follows, we call this a “self-contained push up”.)

Figure 5: The computations in the boxes are self contained.

5 Inverse tft Algorithm

Finally, in this section, we present a simple recursive description of the inverse tft algorithm for the case we have restricted ourselves to (all zeroes packed at the end). The algorithm operates in a length $n$ array $x = (x_{0}, \dots, x_{n - 1})$ for which we assume access; here $n = 2^{p}$ corresponds to $ω$ , a $n$ th primitive root of unity. Initially, the content of the array is

x := (x_{p, 0}, \dots, x_{p, ℓ - 1}, 0, \dots, 0)

where $(x_{p, 0}, \dots, x_{p, ℓ - 1})$ is the result of the tft on $(x_{0, 0}, \dots, x_{0, ℓ - 1}, 0, \dots, 0)$ .

In keeping with our “Illustrated” description we use pictures, like Figure 6, to indicate what values are known (solid dots $∙$ ) and what value to calculate (empty dot $\circ$ ). For instance, “push down $x_{k}$ with Figure 6”, is shorthand for: use $x_{k} = x_{s - 1, i m_{s} + j}$ and $x_{k + m_{s} + j} = x_{s - 1, (i + 1) m_{s} + j}$ to determine $x_{s, i m_{s} + j}$ . We emphasize with an arrow that this new value should also overwrite the one at $x_{k}$ . This calculation is easily accomplished using (1) with a caveat: the values $i$ and $j$ are not explicitly known. What is known is $s$ , and therefore $m_{s}$ , and some array position $k$ . Observe that $i$ is recovered by $i = k quo m_{s}$ (the quotient of $k / m_{s}$ ).

Figure 6: Overwrite $x [i m_{s} + j]$ with $x [i m_{s} + j] + ω^{[i]_{s} m_{s}} x [(i + 1) m_{s} + j]$ .

The full description of the inverse tft follows in Algorithm 1; note that the initial call is InvTFT $(0, ℓ - 1, n - 1, 1)$ . A visual depiction of Algorithm 1 is given in Figure 9. A sketch of a proof of its correctness follows.

Algorithm 1 InvTFT $(head, tail, last, s)$

Theorem 2.

Algorithm 1, initially called with $I n v T F T (0, ℓ - 1, n - 1, 1)$ and given access to the zero-indexed length $n$ array

x = (x_{p, 0}, \dots, x_{p, ℓ - 1}, 0, \dots, 0)

(3)

will terminate with

x = (x_{0, 0}, \dots, x_{0, ℓ - 1}, 0, \dots, 0)

(4)

where (3) is the result of the tft on (4).

Termination.

Let ${head}_{i}$ , ${tail}_{i}$ , and ${last}_{i}$ be the values of head, tail, and last at the $i$ th recursive call. Consider the integer sequences given by

	$α_{i}$	$= {tail}_{i} - {head}_{i} \in Z,$
	$β_{i}$	$= {tail}_{i} - ⌊ \frac{{last}_{i} - {head}_{i}}{2} + head ⌋ \in Z .$

If ${head}_{i} > {tail}_{i}$ then we have termination. Otherwise, either branch (7) executes giving

	${head}_{i + 1}$
	${tail}_{i + 1}$	$= {tail}_{i}$
	${last}_{i + 1}$	$= {last}_{i}$

and thus $α_{i + 1} < α_{i}$ , or branch (13) executes, giving

	${head}_{i + 1}$	$= {head}_{i}$
	${tail}_{i + 1}$	$= {tail}_{i}$
	${last}_{i + 1}$

and thus $β_{i + 1} < β_{i}$ .

Neither branch can run forever since $α < 0$ causes termination and $β < 0$ means either $α$ strictly decreases or condition (13) fails, forcing termination. ∎

Sketch of correctness.

Figure 7 and Figure 8 demonstrate that self contained regions can be exploited to obtain the initial values required to complete the inversion. That is to say, for $n \in {0, \dots, p - 1}$ , that

x = (x_{p - n - 1, 0}, \dots, x_{p - n - 1, 2^{n + 1} - 1})

can always be calculated from

x = (x_{p, 0}, \dots, x_{p, ℓ - 1}, x_{p - n - 1, t a i l + 1}, \dots, x_{p - n - 1, 2^{n + 1} - 1}) .

∎

(a) Line (8): push up the self contained (dashed) region. This yields values sufficient to push down at line (9).

(d) Initially there is sufficient information to push down at line (14).

(a) Initial state of the algorithm. Grey dots are the result of the forward tft; larger grey dots are zeros.

6 Conclusions

The Truncated Fourier Transform is a novel and elegant way to reduce the number of computations of a dft-based computation by a possible factor of two (which may be significant). Additionally, with the advent of Harvey and Roche’s paper [2], it is possible to save as much space as computation. The hidden “cost” of working with the tft algorithm is the increased difficulty of determining the inverse tft. Although in most cases this is still less costly than the inverse dft, the algorithm is no doubt more difficult to implement.

Acknowledgements

The author wishes to thank Dr. Dan Roche and Dr. Éric Schost for reading a draft of this paper and offering suggestions.

References

[1] K. O. Geddes, S. R. Czapor, and G. Labahn. Algorithms for Computer Algebra. Kluwer Academic Publishers, 1992.
[2] David Harvey and Daniel S. Roche. An in-place truncated fourier transform and applications to polynomial multiplication. In Proceedings of the 2010 International Symposium on Symbolic and Algebraic Computation, ISSAC ’10, pages 325–329, New York, NY, USA, 2010. ACM.
[3] H.V. Sorensen and C.S. Burrus. Efficient computation of the dft with only a subset of input or output points. Signal Processing, IEEE Transactions on, 41(3):1184 –1200, mar 1993.
[4] J. van der Hoeven. Notes on the Truncated Fourier Transform. Technical report, Université Paris-Sud, Orsay, France, 2008.
[5] Joris van der Hoeven. The truncated fourier transform and applications. In ISSAC ’04: Proceedings of the 2004 international symposium on Symbolic and algebraic computation, pages 290–296, New York, NY, USA, 2004. ACM.

An Illustrated Introduction to the Truncated Fourier Transform

Authors

The inverse cyclotomic Discrete Fourier Transform algorithm

An in-place truncated Fourier transform

Interesting Open Problem Related to Complexity of Computing the Fourier Transform and Group Theory

A promotion for odd symmetric discrete Fourier transform

Deriving the Variance of the Discrete Fourier Transform Test Using Parseval's Theorem

midiVERTO: A Web Application to Visualize Tonality in Real Time

Numerical Algorithms for the Computation of Generalized Prolate Spheroidal Functions

1 Introduction

2 The Discrete Fourier Transform

Definition.

Example.

3 The Truncated Fourier Transform

Remark.

Theorem 1.

Proof.

4 Inverting The Truncated Fourier Transform

Example.

5 Inverse tft Algorithm

Theorem 2.

Termination.

Sketch of correctness.

6 Conclusions

Acknowledgements

References

An Illustrated Introduction to the Truncated Fourier Transform

Authors

Related Research

The inverse cyclotomic Discrete Fourier Transform algorithm

An in-place truncated Fourier transform

Interesting Open Problem Related to Complexity of Computing the Fourier Transform and Group Theory

A promotion for odd symmetric discrete Fourier transform

Deriving the Variance of the Discrete Fourier Transform Test Using Parseval's Theorem

midiVERTO: A Web Application to Visualize Tonality in Real Time

Numerical Algorithms for the Computation of Generalized Prolate Spheroidal Functions

This week in AI

1 Introduction

2 The Discrete Fourier Transform

Definition.

Example.

3 The Truncated Fourier Transform

Remark.

Theorem 1.

Proof.

4 Inverting The Truncated Fourier Transform

Example.

5 Inverse tft Algorithm

Theorem 2.

Termination.

Sketch of correctness.

6 Conclusions

Acknowledgements

References