An Algorithm for the 2D Radix-2 Sliding Window Fourier Transform
LEE RICHARDSON, Department of Statistics, Carnegie Mellon University
WILLIAM F. EDDY, Department of Statistics, Carnegie Mellon University
We present a new algorithm for the 2D Radix-2 Sliding Window Fourier Transform (SWFT). Our algorithm avoids

arXiv:1707.08213v1 [] 25 Jul 2017

repeating calculations in overlapping windows by using a tree representation of the Cooley-Tukey Fast Fourier Transform
(FFT). For an N 0 × N 1 array and n 0 = 2m0 × n 1 = 2m1 windows, our algorithm takes O(N 0 N 1n 0n 1 ) operations, which
is faster than taking a 2D FFT in each window. We provide a C implementation of the algorithm, compare ours with
existing algorithms, and show how the algorithm extends to higher dimensions.
Additional Key Words and Phrases: Fourier, Tree

1

INTRODUCTION

In the early 19t h century, Joseph Fourier discovered how to decompose a function into a sum of sine and
cosine waves with different frequencies, now called the Fourier transform. The Fourier transform burst into
the digital age when [Cooley and Tukey, 1965] published the Fast Fourier Transform (FFT) algorithm. Since
then, applications of the Fourier transform have soared ([Bracewell, 1986]).
The FFT reduces the number of operations for the discrete Fourier transform (DFT) of a 1D length n signal
from O(n 2 ) to O(n log n). The FFTs computational savings are so significant that the Society of Industrial
and Applied Mathematics (SIAM) named the FFT as a top-ten algorithm of the 20th century ([Cipra, 2000]).
That said, a downside of the FFT is that it operates on a global signal, meaning that local in time information
is lost. Desiring locality, [Gabor, 1946] introduced a transform balancing time and frequency, which we call
the Sliding Window Fourier Transform (SWFT).
After [Gabor, 1946], a plethora of researchers developed a variety of algorithms for the 1D Sliding Window
Fourier Transform (1D SWFT). Given the range of algorithms, we review the literature in Section 2, tying
together previous work and major developmental themes.
Authors’ addresses: Department of Statistics, Carnegie Mellon University, 5000 Forbes Ave, Pittsburgh, PA 15213.
Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that
copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first
page. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To
copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request
permissions from permissions@acm.org.
© 2017 Association for Computing Machinery.
Manuscript submitted to ACM
Manuscript submitted to ACM

1

2

Richardson, Eddy
Today, Fourier transform applications extend beyond 1D. Like 1D, 2D Fourier transforms operate globally,

but can capture local information using a 2D SWFT. This paper presents a new algorithm for the 2D SWFT.
Compared with existing algorithms, our algorithm is fast, stable, interpretable, and the only one with a
publicly available software implementation.
The rest of this paper is organized as follows. Section 2 reviews existing SWFT algorithms. Section 3
describes the 1D algorithm we extend in this paper. Section 4 derives our algorithm, discusses implementation,
compares ours with existing algorithms, and shows how our algorithm extends to higher dimensions. Section
5 concludes with a brief discussion.
2

PREVIOUS WORK

Since the 1960s, researchers have produced two different types of algorithms for the SWFT: recursive and
non-recursive. Recursive algorithms update DFT coefficients from previous windows using both new data
and the Fourier shift theorem, and non-recursive algorithms use repeated FFT calculations in overlapping
windows. This section briefly reviews the history of both algorithm types, starting with the recursive
algorithm, then the non-recursive algorithm, and concluding with recent developments in 2D.
The first recursive algorithm for the SWFT was introduced by [Halberstein, 1966], and has since been
rediscovered many times (e.g. [Amin, 1987, Aravena, 1990, Bitmead, 1982, Bongiovanni et al., 1976, Hostetter,
1980, Lilly, 1991]). Both [Sherlock and Monro, 1992] and [Unser, 1983] gave 2D versions of the recursive
algorithm, and [Sherlock, 1999] provided a 2D C implementation. [Sorensen and Burrus, 1988] generalized
the recursive algorithm to situations where the window moves more than one position, and [Park and
Ko, 2014] improved this generalization in the article, “The Hopping DFT”. [Macias and Exposito, 1998],
[Albrecht et al., 1997], and [Albrecht and Cumming, 1999] all proposed improvements to the recursive
algorithm. Finally, the most cited paper on the recursive algorithm is [Jacobsen and Lyons, 2003], likely
because this paper provides an excellent description.
One downside of the recursive algorithm is numerical error. In fact, [Covell and Richardson, 1991] proved
that the variance of the numerical error is unbounded. Researchers responded by proposing numerically
stable adaptations (e.g. [Douglas and Soh, 1997, Duda, 2010, Jacobsen and Lyons, 2003]). Although most of
these adapted algorithms substantially increase computational complexity, [Park, 2015b] recently proposed
the fastest, numerically stable recursive algorithm, called the “Guaranteed Stable SWFT”. [van der Byl and
Inggs, 2016] recently reviewed recursive SWFT algorithms.
A numerically stable alternative to the recursive algorithms are the non-recursive algorithms. Where
the recursive algorithm updates DFT coefficients from previous windows, the non-recursive algorithm
calculates an FFT in each window, and avoids repeating calculations already computed for previous windows.
[Covell and Richardson, 1991] first proposed the non-recursive algorithm, followed shortly by [Farhang and
Lim, 1992], who generalized the algorithm to arbitrary sized shifts ([Farhang-Boroujeny and Gazor, 1994]).
These authors illustrated their algorithms using a “Butterfly Diagram”, which shows calculations of the FFT
algorithm (see Figure 1). Specifically, both authors showed that calculations in the butterfly diagram of two
Manuscript submitted to ACM

An Algorithm for the 2D Radix-2 Sliding Window Fourier Transform

3

overlapping windows duplicated. Since then, [Exposito and Macias, 2000, 1999] improved the algorithm
for use in a digital relaying application. [Exposito and Macias, 1999] also pointed out that the butterfly
diagram forms a binary tree, which was known (e.g. [Van Loan, 1992]), but this was the first time the tree
was connected with the SWFT. Recently, [Montoya et al., 2012] and [Montoya et al., 2014] gave further
improvements to the non-recursive algorithm, including extension to the Radix-4 case.
Independently, [Wang et al., 2009] discovered a non-recursive algorithm during a magnetoencephalography experiment. These authors capitalized on the binary tree structure of the FFT algorithm, using it to
derive a 3D data-structure shaped like a long triangular prism. [Wang and Eddy, 2012] further developed a
parallel version of their algorithm. This paper extends the Wang-Eddy algorithm to 2D.
Recently, researchers have proposed algorithms for the 2D SWFT. [Park, 2015a] extended the 1D recursive
algorithm to 2D, and [Byun et al., 2016] proposed a 2D SWFT based on the 2 × 2 Vector-Radix algorithm
([Rivard, 1977] and [Harris et al., 1977]).

3

THE 1D FAST SLIDING WINDOW FOURIER TRANSFORM

This section describes the algorithm in [Wang et al., 2009], which we extend in Section 4. For convenience,
we call this the 1D Fast SWFT. After defining the 1D SWFT, we show a tree representation of the FFT, then
illustrate how the tree representation leads to the fast algorithm.

3.1

The 1D Sliding Window Fourier Transform

The 1D SWFT takes a sequence of DFTs within each position of a sliding window. Specifically, let x =
[x 0 , x 1 , . . . , x N −1 ] be a length N complex-valued signal, and let n ≤ N be the window size. Indexing the
window position by p = n − 1, n, . . . N − 1, the 1D SWFT is:
p

ak

=

n−1
1Õ
xp−n+1+j ωn−jk
n j=0

(1)

2π
2π
for k = 0, 1, . . . N − 1, where ωn = exp i2π
n = cos N + i sin N .

A straightforward calculation of Equation 1 takes Pn 2 operations, where P = N − n + 1, the number of

window positions. Replacing the DFT with an FFT for each window position reduces this to O(Pn log n)
operations. The fast algorithm described here further reduces this to O(Pn).
Before introducing the algorithm, we clarify a few points regarding the FFT algorithm used in this paper.
First, while many different FFT algorithms have been developed (e.g. [Rabiner et al., 1969] and [Good, 1958]),
we focus on the [Cooley and Tukey, 1965] algorithm. This means whenever we say FFT, we are referring to
the Cooley-Tukey algorithm. Next, since the FFT algorithm factorizes a length n signal, different algorithms
exist for different n. This paper only considers in detail when n is a power of two, called the Radix-2 case.
Manuscript submitted to ACM

4

Richardson, Eddy

Fig. 1. A butterfly diagram for n = 8. The squares on the left are the input data, and the circles on the right are the
DFT coefficients. Multiplication takes place at the beginning of each arrow, and addition takes place at the end. For
example, the first calculation on the bottom left means that x 7 is multiplied by ω 4 , then added to x 3 .

3.2

Butterflies, Overlapping Trees, and a Fast Algorithm

The 1D Fast SWFT algorithm takes a 1D FFT in each window, and avoids repeating calculations already
computed in previous windows. Understanding which calculations have already been computed requires
detailed understanding of the FFT. Figure 1 shows the famous butterfly diagram, giving FFT calculations for
a length n = 8 signal. The squares on the left of Figure 1 correspond to the input data, the circles on the
right are the output DFT coefficients, and the arrows in the middle are the calculations. Both [Covell and
Richardson, 1991] and [Farhang and Lim, 1992] derived their non-recursive algorithms by showing that
calculations in the butterfly diagram repeat in overlapping windows.
An equivalent tree representation of the FFT is shown in the top panel of Figure 2. It is important to
stress that the calculations in the butterfly diagram and this tree representation are the same. For sliding
windows, the tree representation has two major advantages. First, the input data and output coefficients are
ordered. Second, underneath x 7 is a binary tree with three (since 23 = 8) levels, and the final level of this
binary tree contains the DFT coefficients.
Figure 2 demonstrates the 1D Fast SWFT algorithm. The top panel shows calculations for the first window
position with input data [x 0 , x 1 , . . . , x 7 ], and the bottom panel shows the calculations for window position
two with input data [x 1 , x 2 , . . . , x 8 ]. In window two, solid arrows represent computations made in window
one, and dashed arrows represent new window two calculations. The number of window two calculations
Manuscript submitted to ACM

An Algorithm for the 2D Radix-2 Sliding Window Fourier Transform

5

Fig. 2. The calculations of the 1D Fast SWFT for the first two length n = 8 windows. The top shows a tree representation
of FFT calculations for input data [x 0 , x 1 , . . . , x 7 ]. In window one, the DFT coefficients are circles at the third level of
the binary tree underneath x 7 . The bottom shows FFT calculations in window two, where solid arrows indicate window
one calculations, and dashed arrows indicate window two calculations. In window two, the large square indicates the
new data-point (x 8 ) entering the window, and the small square indicates the data-point (x 0 ) leaving the window.

is exactly the size of the binary tree, and the difference between the number of solid and dashed arrows is
exactly the log factor speed up gained in the 1D Fast SWFT algorithm.
Implementing the 1D Fast SWFT requires calculating each node of each binary tree, corresponding to
three nested loops: 1. over N trees 2. over log2 (n) levels 3. over 2l nodes at level l of each tree. The only
restriction on loop order is that loop 2 (over levels) must precede loop 3 (over nodes), since nodes at lower
levels of the tree depend on nodes at higher levels.
4

THE 2D FAST SLIDING WINDOW FOURIER TRANSFORM

This section gives the 2D Fast SWFT algorithm. After defining the 2D SWFT, we derive the algorithm. We
then discuss implementation, compare ours with existing algorithms, and show how the algorithm extends
to higher dimensions.
4.1

The 2D Sliding Window Fourier Transform

The 2D SWFT of an N 0 × N 1 array calculates a 2D DFT for all n 0 × n 1 windows. Our algorithm requires
that n 0 = 2m0 and n 1 = 2m1 , the Radix-2 case. Let x be an N 0 × N 1 array. There are Pi = Ni − ni + 1; i = 0, 1
windows in each direction, making P0 P1 total window positions. Indexing the window position by (p0 , p1 ),
where pi = ni − 1, ni , . . . , Ni − 1, the 2D SWFT is:
p ,p

ak00,k11

=

n 0 −1 nÕ
1 −1
1 Õ
xp −n +1+j0,p1 −n1 +1+j1 ωn−j0 0 k0 ωn−j1 1 k1
√
n 0n 1 j =0 j =0 0 0
0

(2)

1

where ki = 0, 1, . . . , ni − 1; i = 0, 1. Equation 2 outputs a P0 × P1 × n 0 × n 1 array.
Manuscript submitted to ACM

6

Richardson, Eddy

Fig. 3. The 2D Fast SWFT for two windows, where window one has input data [x 0,0 , x 0,1 , x 0,2 , x 0,3 , x 1,0 , x 1,1 , x 1,2 , x 1,3 ]
and window two replaces column 0 (x 0,0 , x 1,0 ) with column 4 (x 0,4 , x 1,4 ). The circles on level three of the binary tree
underneath the bottom-right points of each window are the DFT coefficients. In the bottom panel, solid lines are
calculations from window one, and dashed lines are calculations from window two.

A straightforward calculation of Equation 2 takes P0 P 1n 20n21 operations. Replacing the 2D DFT with a

2D FFT (described next) reduces the number of operations to O(P 1 P0n 0n 1 log(n 0n 1 )). Our algorithm further
reduces this to O(P0 P1n 0n 1 ).
4.2

Derivation

Conceptually, the 2D Fast SWFT algorithm works like 1D. As in 1D, there is a tree representation of the
2D FFT (see top panel of Figure 3). Also like 1D, we take a 2D FFT in each window, and avoid repeating
calculations computed in previous windows. Figure 3 illustrates how the 2D Fast SWFT works. The solid
arrows in window two indicate 2D FFT calculations from window one, and the dashed arrows are the new
calculations required for window two. Once again, the number of dashed lines is exactly the size of the tree
underneath the bottom-right point (x 1,4 ) of window two.
Though conceptually similar, several details differentiate the 1D and 2D algorithms. First, the trees
underneath each data-point are two-dimensional, where the two dimensions are depths m 0 = log2 (n 0 ) and
m 1 = log2 (n 1 ). Second, we now require two twiddle-factor vectors; Twiddle factors (named by [Gentleman
and Sande, 1966]) are the trigonometric constants used for combining smaller DFTs during the FFT algorithm
(the ω’s in Figures 1 and 2 are twiddle factors). Finally, when the window slides by one position, we now
replace an entire row or column, opposed to a single data-point in 1D.
Manuscript submitted to ACM

An Algorithm for the 2D Radix-2 Sliding Window Fourier Transform

7

The key to resolving differences between the 1D and 2D algorithms is the “factorization” property of the
2D DFT. Factorization means that 2D DFTs can be computed using only 1D DFTs:
ak0,k1 =
where z j0,k1 =

Ín1 −1
j 1 =0

nÕ
0 −1 n
1 −1
Õ
j 0 =0 j 1 =0

x j0, j1 ωn−j0 0 k0 ωn−j1 1 k1 =

nÕ
0 −1
j 0 =0

z j0,k1 ωn−j0 0 k0

x j0, j1 ωn−j1 1 k1 . This implies the 2D FFT can be computed by first taking 1D FFTs of each

row, followed by 1D FFTs of the resulting columns, sometimes called the “row-column” FFT. Factorization
lets us update our algorithm when the 2D FFT switches from 1D row FFTs to 1D column FFTs. With this
preamble, we are now ready to derive the 2D Fast SWFT algorithm.
Following the above notation, let x be a N 0 × N 1 array, with window sizes n 0 × n 1 , where n 0 = 2m0 and
n 1 = 2m1 , the Radix-2 case. We do not require that n 0 = n 1 . Our trees have m 0 + m 1 levels: the first m 0 levels
correspond to row FFTs, the next m 1 levels correspond to column FFTs, and the final level (m 0 + m 1 ) is
i −1)
the 2D DFT coefficients. Let Ωi = [ωn0 i , ωn−1i , . . . , ωn−(n
] be a length ni twiddle-factor vector, for i = 0, 1.
i

Finally, let T be a tree data-structure, used to store the 2D FFT calculations for each window position. We
access T by window position, level, and node. For example, Tp0,p1,l,i 0,i 1 corresponds to node (i 0 , i 1 ) on the

l t h level of the tree at window position (p0 , p1 ).

We give the calculations for an arbitrary tree at position (p0 , p1 ). Each level l has 2l nodes. Level zero
of the tree is the data: Tp0,p1,0,0,0 = xp0,p1 . We use one additional piece of notation: let sld be the “shift” for
dimension d and level l. The “shift” (sld ) helps us identify which tree has the repeated calculation needed for
the current calculation. For example, if s 10 = 2, this means for level 1, the repeated calculation is located at

tree (p0 − 2, p1 ). For levels corresponding to row FFTs: sl1 = 2m1 −l , and for column FFT levels: sl0 = 2m0 +m1 −l .
Level one of tree (p0 , p1 ) has 21 = 2 nodes: (0, 0), (0, 1). The shift is s 11 = 2m1 −1 , meaning the repeated

calculation is at tree (p0 , p1 − s 11 ). For node (0, i 1 ), the calculation is a complex multiplication between node

(0, i 1 mod 21−1 ) = (0, 0) at level 1−1 = 0 of the current tree (p0 , p1 ), and the i 1 ·s 11 element from twiddle-factor
vector Ω1 . The complex multiplication output is then added to node (0, 0) from the shifted tree (p0 , p1 − s 11 ).
Using our tree notation:

Tp0,p1,1,0,i 1

= Tp0,p1 −s11,0,0,0 + Ωi1 ·s 1 · Tp0,p1,0,0,0
1

1

for i 1 = 0, 1, which is a length two DFT of xp0,p1 and xp0,p1 −s11 .

Calculating each node of T has the same form: TA = TB + ΩCi · TD . A, B, C and D are indices into either T,

Ω0 , or Ω 1 . The exact indices depend on whether the level of the tree corresponds to the row (1 to m 1 ) or

column (m 1 + 1 to m 1 + m 0 ) part of the 2D FFT algorithm. We give the exact indices for both situations next.
Level t ≤ m 1 of tree (p0 , p1 ) has 2t nodes: 1 in the row-direction, and 2t in the column-direction. The
repeated calculation comes from tree (p0 , p1 − st1 ), where st1 = 2m1 −t . For node (0, i 1 ); i 1 = 0, 1, . . . 2t − 1, the
calculation is:
Manuscript submitted to ACM

8

Richardson, Eddy

Tp0,p1,t,0,i 1

= Tp0,p1 −st1,t −1,0,i 1 mod 2t −1 + Ωi1 ·s 1 · Tp0,p1,t −1,0,i mod 2t −1
1

t

(3)

Level v > m 1 of tree (p0 , p1 ) has 2v nodes: 2v−m1 in the row-direction, and n 1 in the column-direction. The
repeated calculation comes from tree (p0 −sv0 , p1 ), where sv0 = 2m0 +m1 −l . For node (i 0 , i 1 ); i 0 = 0, 1, . . . n 1 , i 1 =
0, 1, . . . 2v−m1 − 1, the calculation is:
Tp0,p1,v,i 0,i 1

= Tp0 −sv0 ,p1,v−1,i 0 mod 2v −m1 −1,i 1 + Ωi0 ·s 0 · Tp0,p1,v−1,i 0 mod 2v −m1 −1,i 1
0

v

(4)

The final level (m 0 + m 1 ) contains the 2D DFT coefficients for window (p0 , p1 ). After calculating each
tree, all that remains is selecting the subset at level m 0 + m 1 of each tree, and the algorithm is complete.
4.3

Implementation

We have written a C function implementing the 2D algorithm. After creating the twiddle-factor vectors and
allocating memory, the core of our algorithm is six nested loops:
(1) over m 0 = log2 (n 0 ) levels corresponding to row FFTs
(2) over m 1 = log2 (n 1 ) levels corresponding to column FFTs
(3) over N 0 trees in the row-direction
(4) over N 1 trees in the column-direction
(5) over 2l0 nodes in the row-direction at level l 0
(6) over 2l1 nodes in the column-direction at level l 1
Inside the six loops is either Equation 3 or 4, depending on whether the level corresponds to the row or
column part of the 2D FFT algorithm.
Since the next level of the trees only depends on the previous level, our implementation allocates
2N 0 N 1n 0n 1 complex numbers of memory: N 0 N 0n 0n 1 for both the previous and current level. This is possible
because the two loops over levels are innermost. This also implies our algorithm is suitable for parallel
computing, following [Wang and Eddy, 2012].
Like 1D, we can swap the order of the loops. For our derivation, the only restrictions are that loop 1
must precede loop 2, and loops 5 and 6 must come after loops 1 and 2. A particularly interesting order is
when loops 3 and 4 (over trees) are innermost. This version can be tailored to a real-time task, opposed
to the “levels innermost” order, which requires all data to be available. For example, if the channels of a
multi-spectral scanner are considered as a two-dimensional array, and T is stored, we could calculate a new
window in O(n 0n 1 ) time.
The 2D SWFT is memory intensive, since the output is a P 0 ×P1 ×n 0 ×n 1 array. For example, if P0 = P 1 = 400,
and n 0 = n 1 = 32, then the output is 163, 840, 000 complex numbers. Because our implementation uses
Manuscript submitted to ACM

An Algorithm for the 2D Radix-2 Sliding Window Fourier Transform

9

16 byte double complex numbers, the output alone requires 16 · 163840000 = 2621440000 bytes, or ≈ 2.6
gigabytes.
These ideas extend directly to kD. For instance, if we had an N 0 ×N 1 ×. . .×Nk−1 array, with ni = 2mi ≤ Ni
windows, then the calculation for node (i 0 , i 1 , . . . , i k−1 ), level l, and dimension c is:
Tp0, ..,pk −1,l,i 0, ..,i k −1

= T
+

where slc = 2(

Ík −1
i =c

m i )−l ,

Ík −1

p0, ..,pc −slc , ..,pk −1,l −1,0, ..,i c mod 2l −1− i =c +1 mi , ..,i k −1
Ωicc ·s c · Tp0, ...,pk −1,l −1,0, ..,i mod 2l −1, ..,i k −1
l

and dimension c contains levels [

Ík −1

i=c+1 m i

+ 1,

Ík −1
i=c

(5)

mi ] of the trees. Extending

the implementation to higher dimensions is possible following Equation 5. However, it seems this would
require a separate program for each dimension.

4.4

Results and Comparisons

Like 1D, the run-time of the 2D Fast SWFT algorithm grows linearly in window (n 0n 1 ) and array (N 0 N 1 )
size. This is because calculating each node requires one complex multiplication and addition, defined as an
operation (the same definition was used in [Cooley and Tukey, 1965]). Since level l requires 2l operations,
the computational complexity per window is:
Ceach

=

mÕ
0 +m 1

2l = (2 + 4 + . . . + 2m0 +m1 ) = 2(2m0 +m1 − 1) = 2(n 0n 1 − 1)

l =1

The exact run-time is slightly more complicated, since windows with indices either less than n 0 in the row
direction or less than n 1 in the column direction do not require complete trees. These extra calculations are
negligible for large arrays, so the 2D algorithm takes O(P0 P1n 0n 1 ) operations. Table 1 compares the speed,
memory, and properties of the 2D Fast SWFT with existing algorithms. Our memory numbers for [Park,
2015a] and [Byun et al., 2016] come from the Table 2 in [Byun et al., 2016]. Out of existing O(P0 P1n 0n 1 )
algorithms, our 2D Fast SWFT is the only that is numerically stable, works for non-square windows, and
has an existing software implementation.

5

DISCUSSION

The goal of this paper is to describe our algorithm as clearly as we can, and instantiate it with a C program.
Our goal is not providing an optimized implementation of the algorithm, such as [Frigo and Johnson, 2005].
Finally, while our paper focused on the Radix-2 case, there is nothing conceptually preventing extension
to other factorizations. For instance, the Radix-3 implementation would simply replace binary with ternary
trees.
Manuscript submitted to ACM

10

Richardson, Eddy
Algorithm

Speed

Memory

2D DFT
O(P0 P1n 20n 21 )
P 0 P 1n 0n 1
2D FFT
O(P0 P1n 1n 0 log(n 0n 1 ))
P 0 P 1n 0n 1
[Park, 2015a]
O(P0 P1n 0n 1 )
P0 P1n 0n 1 + 2n + 1
[Byun et al., 2016]
O(P0 P1n 0n 1 )
P0 P1n 0n 1 + n 2 2 log2 (n − 1) + 1
This Paper
O(P0 P1n 0n 1 )
2N 0 N 1n 0n 1

Stable
Yes
Yes
No
Yes
Yes

Non-Square
Windows
Yes
Yes
Yes
No
Yes

Table 1. Speed, memory, and properties of existing 2D SWFT algorithms. The 2D DFT and FFTs take a 2D DFT or
2D FFT in each window. [Park, 2015a] extends the 1D recursive algorithm to 2D, and [Byun et al., 2016] gives a
Vector-Radix 2 × 2 algorithm. The memory results come from [Byun et al., 2016], and we use n instead of n 0n 1 , because
the algorithms in this paper were derived using non-square window sizes.

A

SOFTWARE

Our algorithm is available online at: http://stat.cmu.edu/~lrichard/fswft.tar.gz

Manuscript submitted to ACM

An Algorithm for the 2D Radix-2 Sliding Window Fourier Transform

11

REFERENCES
Sandor Albrecht and Ian Cumming. 1999. Application of momentary Fourier transform to SAR processing. IEE Proceedings-Radar,
Sonar and Navigation 146, 6 (1999), 285–297.
S Albrecht, I Cumming, and J Dudas. 1997. The momentary Fourier transformation derived from recursive matrix transformations. In
Digital Signal Processing Proceedings, 1997. DSP 97., 1997 13th International Conference on, Vol. 1. IEEE, 337–340.
Moeness G Amin. 1987. A new approach to recursive Fourier transform. Proc. IEEE 75, 11 (1987), 1537–1538.
Jorge L. Aravena. 1990. Recursive moving window DFT algorithm. IEEE Trans. Comput. 39, 1 (1990), 145–148.
R Bitmead. 1982. On recursive discrete Fourier transformation. IEEE Transactions on Acoustics, Speech, and Signal Processing 30, 2
(1982), 319–322.
Giancarlo Bongiovanni, Paolo Corsini, and Graziano Frosini. 1976. Procedures for computing the discrete Fourier transform on
staggered blocks. IEEE Transactions on Acoustics, Speech, and Signal Processing 24, 2 (1976), 132–137.
Ronald Newbold Bracewell. 1986. The Fourier transform and its applications. Vol. 31999. McGraw-Hill New York.
Keun-Yung Byun, Chun-Su Park, Jee-Young Sun, and Sung-Jea Ko. 2016. Vector Radix 2× 2 Sliding Fast Fourier Transform. Mathematical
Problems in Engineering 2016 (2016).
Barry A Cipra. 2000. The best of the 20th century: editors name top 10 algorithms. SIAM news 33, 4 (2000), 1–2.
James W Cooley and John W Tukey. 1965. An algorithm for the machine calculation of complex Fourier series. Mathematics of
computation 19, 90 (1965), 297–301.
M Covell and J Richardson. 1991. A new, efficient structure for the short-time Fourier transform, with an application in code-division
sonar imaging. In Acoustics, Speech, and Signal Processing, 1991. ICASSP-91., 1991 International Conference on. IEEE, 2041–2044.
SC Douglas and JK Soh. 1997. A numerically-stable sliding-window estimator and its application to adaptive filters. In Signals, Systems
& Computers, 1997. Conference Record of the Thirty-First Asilomar Conference on, Vol. 1. IEEE, 111–115.
Krzysztof Duda. 2010. Accurate, guaranteed stable, sliding discrete Fourier transform [DSP tips & tricks]. IEEE Signal Processing
Magazine 27, 6 (2010), 124–127.
A Gomez Exposito and JAR Macias. 2000. Fast non-recursive computation of individual running harmonics. IEEE Transactions on
Circuits and Systems II: Analog and Digital Signal Processing 47, 8 (2000), 779–782.
Antonio Gomez Exposito and Jose A Rosendo Macias. 1999. Fast harmonic computation for digital relaying. IEEE Transactions on
Power Delivery 14, 4 (1999), 1263–1268.
Behrouz Farhang and YC Lim. 1992. Comment on the computational complexity of sliding FFT. IEEE Transactions on Circuits and
Systems II: Analog and Digital Signal Processing 39, 12 (1992), 875–876.
Behrouz Farhang-Boroujeny and Saeed Gazor. 1994. Generalized sliding FFT and its application to implementation of block LMS
adaptive filters. IEEE Transactions on Signal Processing 42, 3 (1994), 532–538.
Matteo Frigo and Steven G. Johnson. 2005. The Design and Implementation of FFTW3. Proc. IEEE 93, 2 (2005), 216–231. Special issue
on “Program Generation, Optimization, and Platform Adaptation”.
Dennis Gabor. 1946. Theory of communication. Part 1: The analysis of information. Electrical Engineers-Part III: Radio and Communication Engineering, Journal of the Institution of 93, 26 (1946), 429–441.
W Morven Gentleman and Gordon Sande. 1966. Fast Fourier Transforms: for fun and profit. In Proceedings of the November 7-10, 1966,
fall joint computer conference. ACM, 563–578.
Irving John Good. 1958. The interaction algorithm and practical Fourier analysis. Journal of the Royal Statistical Society. Series B
(Methodological) (1958), 361–372.
JH Halberstein. 1966. Recursive, complex Fourier analysis for real-time applications. Proc. IEEE 54, 6 (1966), 903–903.
D Harris, J McClellan, D Chan, and H Schuessler. 1977. Vector radix fast Fourier transform. In Acoustics, Speech, and Signal Processing,
IEEE International Conference on ICASSP’77., Vol. 2. IEEE, 548–551.

Manuscript submitted to ACM

12

Richardson, Eddy

G Hostetter. 1980. Recursive discrete Fourier transformation. IEEE Transactions on Acoustics, Speech, and Signal Processing 28, 2 (1980),
184–190.
Eric Jacobsen and Richard Lyons. 2003. The sliding DFT. IEEE Signal Processing Magazine 20, 2 (2003), 74–80.
JH Lilly. 1991. Efficient DFT-based model reduction for continuous systems. IEEE transactions on automatic control 36, 10 (1991),
1188–1193.
JA Rosendo Macias and A Gomez Exposito. 1998. Efficient moving-window DFT algorithms. IEEE Transactions on Circuits and Systems
II: Analog and Digital Signal Processing 45, 2 (1998), 256–260.
Dan-El A Montoya, JA Rosendo Macias, and A Gomez-Exposito. 2012. Efficient computation of the short-time DFT based on a modified
radix-2 decimation-in-frequency algorithm. Signal Processing 92, 10 (2012), 2525–2531.
Dan-El A Montoya, JA Rosendo Macias, and A Gomez-Exposito. 2014. Short-time DFT computation by a modified radix-4 decimationin-frequency algorithm. Signal Processing 94 (2014), 81–89.
Chun-Su Park. 2015a. 2D discrete Fourier transform on sliding windows. IEEE Transactions on Image Processing 24, 3 (2015), 901–907.
Chun-Su Park. 2015b. Fast, Accurate, and Guaranteed Stable Sliding Discrete Fourier Transform [sp Tips&Tricks]. IEEE Signal
Processing Magazine 32, 4 (2015), 145–156.
Chun-Su Park and Sung-Jea Ko. 2014. The hopping discrete Fourier transform [sp tips&tricks]. IEEE Signal Processing Magazine 31, 2
(2014), 135–139.
L Rabiner, RW Schafer, and C Rader. 1969. The chirp z-transform algorithm. IEEE transactions on audio and electroacoustics 17, 2 (1969),
86–92.
G Rivard. 1977. Direct fast Fourier transform of bivariate functions. IEEE Transactions on Acoustics, Speech, and Signal Processing 25, 3
(1977), 250–252.
BG Sherlock and DM Monro. 1992. Moving discrete Fourier transform. In IEE Proceedings F-Radar and Signal Processing, Vol. 139. IET,
279–282.
Barry G Sherlock. 1999. Windowed discrete Fourier transform for shifting data. Signal Processing 74, 2 (1999), 169–177.
Henrik V Sorensen and C Sidney Burrus. 1988. Efficient computation of the short-time fast Fourier transform. In Acoustics, Speech, and
Signal Processing, 1988. ICASSP-88., 1988 International Conference on. IEEE, 1894–1897.
Michael Unser. 1983. Recursion in short-time signal analysis. Signal Processing 5, 3 (1983), 229–240.
Andrew van der Byl and Michael R Inggs. 2016. Constraining errorâĂŤA sliding discrete Fourier transform investigation. Digital
Signal Processing 51 (2016), 54–61.
Charles Van Loan. 1992. Computational frameworks for the fast Fourier transform. Vol. 10. Siam.
Jianming Wang, Bronwyn Woods, and William F Eddy. 2009. MEG, RFFTs, and the hunt for high frequency oscillations. In Image and
Signal Processing, 2009. CISP’09. 2nd International Congress on. IEEE, 1–5.
Jian-Ming Wang and William F Eddy. 2012. Parallel Algorithm for SWFFT Using 3D Data Structure. Circuits, Systems, and Signal
Processing 31, 2 (2012), 711–726.

Manuscript submitted to ACM

