Skip to content

The 2-D DFT

Compute a 2-D DFT as row then column DFTs, read an image's spectrum, blur it with a mask, and see why phase carries the shapes.

Before thisImage filtering (28.2)

5 more before it

Convergence and the Gibbs phenomenon (7.3), The DFT (13.2), Properties of the DFT (13.3), The FFT (14.1), Images as signals (28.1)

Before this28.2 · 5 more
Chapter 28 · Lesson 3 of 4

First, the picture

How much of each pattern of straight stripes does an image hold? Watch the 1-D DFT find out below, along every row and then every column: first five lines appear, then five dots.

Rows first, then columns

A constant plus two gratings, 64 × 64 pixels.

A constant, a grating of 4 cycles across, and a tilted grating of 3 across and 6 down.

stage
image
nonzero bins
—
0.00 / 12.00 s
Describe this picture

A constant plus two gratings, 64 × 64 pixels, in three square panels. The first is the image in grey levels. The second, after the rows, shades the size of each row’s DFT, divided by 64, from 0 to 0.5, with bins k1k_1 from −32 to 31 across and the rows running down as in the image. The third, after the columns, shades the size divided by 64264^2, from 0 to 0.5, over k1k_1 and k2k_2 from −32 to 31, with 0 at the centre. Each bin that is not zero also gets a ring, so that single bins show. The readouts are the stage and the number of nonzero bins. There is no control. The 12 s clip opens on the image: a constant, a grating of 4 cycles across, and a tilted grating of 3 across and 6 down. From 3 s the row DFTs sweep down the second panel: lines at k1k_1 = −4, −3, 0, 3 and 4, the same strength in every row, 320 nonzero bins. From 7.5 s the column DFTs sweep across the third panel: five dots, 0.500 at the centre and 0.125 at (±4,0)(\pm4,0), (3,6)(3,6) and (−3,−6)(-3,-6).

Rows first, then columns

In Images as signals (28.1), “A frequency with a direction” showed that a grating, a pattern of straight stripes, has a frequency across and a frequency down. Its spectrum was two dots. Real images are not one grating. So this page asks: how much of each grating does an image hold, and how do we find out?

The DFT (13.2) answered that question for a row of samples. In “Each bin: the signal times a probe, added up”, each bin multiplied the signal by a turning probe and added up the products. An image needs the same thing in two directions. You can do it with the 1-D DFT you already know: first along every row, then along every column.

Let’s watch that happen on an image I built from three parts. A constant, 0.5. A grating of 4 cycles across the image. And a tilted grating that makes 3 cycles across and 6 down:

x[n1,n2]=0.5+0.25cos⁡(2π⋅4n1/64)+0.25cos⁡(2π(3n1+6n2)/64).\begin{aligned} &x[n_1,n_2]=0.5\\ &\quad+0.25\cos(2\pi\cdot4n_1/64)\\ &\quad+0.25\cos\big(2\pi(3n_1+6n_2)/64\big). \end{aligned}

As in 28.1, n1n_1 is the column, counted across, and n2n_2 is the row, counted down. The image is 64 × 64 pixels, and its values run from 0.000 to 1.000.

First I take the DFT of each row and divide it by 64. Then I take the DFT of each column of that result and divide by 64 again. The divisions only set the scale: by 13.2, a cosine of height AA whose cycles fit the length exactly puts NA/2NA/2 in each of its two bins. After dividing by NN, each bin holds A/2A/2, and a constant cc holds cc in bin 0.

The picture at the top of the page does exactly this.

Watch the second panel at the columns k1=3k_1=3 and k1=−3k_1=-3. Each looks like a steady line, 0.125 in every row. But the tilted grating’s stripes move along as you go down, so the value in that column turns: its angle goes round 6 times in 64 rows. The column DFTs find that wave and gather it into the single dots at (3,6)(3,6) and (−3,−6)(-3,-6).

The 2-D DFT

Now the formula. For an N×NN\times N image, the 2-D DFT is

X[k1,k2]=∑n2=0N−1∑n1=0N−1x[n1,n2]×e−j2π(k1n1+k2n2)/N.\begin{aligned} &X[k_1,k_2]=\sum_{n_2=0}^{N-1}\sum_{n_1=0}^{N-1}x[n_1,n_2]\\ &\qquad\times e^{-j2\pi(k_1n_1+k_2n_2)/N}. \end{aligned}

Each bin is still the image times a probe, added up. The probe is now a grating: it turns k1k_1 times across the image and k2k_2 times down. In 28.1’s units that is k1/Nk_1/N cycles per pixel across and k2/Nk_2/N down. So X[k1,k2]X[k_1,k_2] measures how much of that one grating the image holds.

The probe splits into two factors, one for each direction:

e−j2π(k1n1+k2n2)/N=e−j2πk1n1/N e−j2πk2n2/N.\begin{aligned} &e^{-j2\pi(k_1n_1+k_2n_2)/N}\\ &\quad=e^{-j2\pi k_1n_1/N}\,e^{-j2\pi k_2n_2/N}. \end{aligned}

The second factor does not depend on n1n_1, so it comes out of the inner sum:

X[k1,k2]=∑n2=0N−1e−j2πk2n2/N×∑n1=0N−1x[n1,n2] e−j2πk1n1/N.\begin{aligned} &X[k_1,k_2]=\sum_{n_2=0}^{N-1}e^{-j2\pi k_2n_2/N}\\ &\qquad\times\sum_{n_1=0}^{N-1}x[n_1,n_2]\,e^{-j2\pi k_1n_1/N}. \end{aligned}

The inner sum is the 1-D DFT of row n2n_2, at bin k1k_1. The outer sum is the 1-D DFT, at bin k2k_2, of what the rows left in column k1k_1. That is the instrument, written down.

A transform that splits like this, one direction at a time, is called separable. The order does not matter: columns first, then rows, gives the same numbers.

The bins run from 0 to N−1N-1 in each direction. By “Bins in hertz, and the mirror” of 13.2, bin N−kN-k turns like bin −k-k. So I draw the bins from −N/2-N/2 to N/2−1N/2-1, with 0 in the centre, as the instrument did.

For a real image the mirror also holds in two dimensions: X[−k1,−k2]=X[k1,k2]∗X[-k_1,-k_2]=X[k_1,k_2]^*. That is why every grating gives a pair of dots, one on each side of the centre.

To get the image back, use the inverse DFT of 13.2 in both directions: e+je^{+j} in place of e−je^{-j}, and 1/N21/N^2 in front. It splits into rows and columns in the same way.

The cost

How much work is the formula? Each of the N2N^2 bins adds up N2N^2 products, so the direct sum takes N4N^4 multiplications. For a 64 × 64 image that is 644=16 777 21664^4=16\,777\,216.

Rows then columns is cheaper. There are NN rows and NN columns, so 2N2N 1-D DFTs of NN points. In The FFT (14.1), “The cost of the DFT” counted N2N^2 multiplications for each. Together that is 2N⋅N2=2N32N\cdot N^2=2N^3, which is 524 288524\,288 for 64 × 64: 32 times fewer.

Then each 1-D DFT can be an FFT. By “How the saving grows” of 14.1, a 64-point FFT costs 642log⁡264=192\tfrac{64}{2}\log_2 64=192 multiplications. The 128 of them cost 24 57624\,576, which is 683 times fewer than the direct sum.

That last number is no accident. A 64 × 64 image has 4096 pixels, and 14.1 found the same 24 57624\,576 for a 1-D FFT of 4096 points. The 2-D FFT costs what a 1-D FFT of all the pixels costs.

Reading a spectrum

Three things show up in the spectrum of an image.

  • The centre is bin (0,0)(0,0). Its probe does not turn at all, so it adds up the pixels. Divided by N2N^2 it is the average brightness, 0.500 in the instrument.
  • A grating is a pair of dots, at (k1,k2)(k_1,k_2) and (−k1,−k2)(-k_1,-k_2). Their distance from the centre is the grating’s frequency. Their direction is at right angles to its stripes, as in 28.1.
  • An edge is a line through the centre, at right angles to the edge.

The third needs a reason. Take an image whose rows are all the same, an edge running straight down from top to bottom.

After the row DFTs, every column holds one value repeated, and a constant’s DFT is a single bin at 0. So everything lands on the row k2=0k_2=0: a line along the across axis. A shorter edge spreads a little, but stays near that line.

Cut the spectrum, blur the image

Now let’s change a spectrum and see what happens to the image. I use the test card of Image filtering (28.2), 128 × 128 pixels. It has a gradient across, a bright disc of 0.9, a dark square of 0.1, and a patch of stripes 0.5+0.3cos⁡(2π(n1+n2)/6)0.5+0.3\cos(2\pi(n_1+n_2)/6). Its values run from 0.100 to 0.900.

Its spectrum shows all three things from the list. The square’s edges run down and across, so they draw a cross along the two axes. The disc’s edge runs in every direction, so it adds faint rings. The stripes repeat every 6 pixels along the diagonal, 1/61/6 cycle per pixel across and down, so they put two bright spots near (21,21)(21,21) and (−21,−21)(-21,-21).

There is one edge that is not drawn. The DFT treats the image as one tile of a pattern that repeats, as the DFT of a row did in Circular vs linear convolution (13.4). So the gradient’s 0.6 at the last column meets its 0.25 at the first column, and that edge adds to the cross too.

To filter, I multiply the spectrum by a mask H[k1,k2]H[k_1,k_2], one number for each bin, and take the inverse DFT:

Y[k1,k2]=H[k1,k2] X[k1,k2].Y[k_1,k_2]=H[k_1,k_2]\,X[k_1,k_2].

My mask keeps every bin within a circle around the centre and removes the rest. Its radius is the cutoff, Ωc\Omega_c, which I quote in cycles per pixel, like 28.1’s frequencies. A cutoff of 0.25 cycles per pixel (Ωc=π/2\Omega_c=\pi/2 rad/pixel) keeps the bins within 0.25⋅128=320.25\cdot128=32 bins of the centre:

H[k1,k2]={1,k12+k22≤32,0,otherwise.H[k_1,k_2]=\begin{cases}1, & \sqrt{k_1^2+k_2^2}\le 32,\\ 0, & \text{otherwise.}\end{cases}

The mask treats (k1,k2)(k_1,k_2) and (−k1,−k2)(-k_1,-k_2) alike, so the mirror survives, and the inverse DFT should come out real. In my NumPy computation its imaginary part is at most 2.5×10−162.5\times10^{-16}, for every cutoff on the slider below. That is rounding only.

Cut the spectrum, blur the image

The 128 × 128 test card, an ideal circular low-pass in the frequency plane.

Cutoff 0.25 cycles/pixel: 96.8 % of the energy kept; the stripes, at 0.236 cycles per pixel, stay, and only finer detail softens.

cutoff
0.25 cycles/pixel
energy kept
96.8 %
RMS change
0.039
0.00 / 13.00 s
Describe this picture

The 128 × 128 test card and an ideal circular low-pass in the frequency plane, in three square panels. The first is the card. The second, the spectrum, shades log⁡10\log_{10} of each bin’s size from −2 to 4, so each step of 1 is ten times larger, with zero frequency in the centre. The mask’s circle is drawn dashed, and the bins outside it are dimmed to 30 %. The third, filtered, is the inverse DFT, with values clipped to 0 to 1 for display. The readouts are the cutoff in cycles per pixel, the energy kept, the percent of the energy outside the centre bin that the circle keeps, and the RMS change, the root mean square of the filtered image minus the card. The 13 s clip opens at a cutoff of 0.25: 96.8 % of the energy kept and an RMS change of 0.039; the stripes, at 0.236 cycles per pixel, stay, and only finer detail softens. From 3.5 s the circle shrinks to 0.1: 85.2 % and 0.083, the stripes are gone and edges blur. From 8 s it shrinks to 0.05: 78.9 % and 0.099, a strong blur with rings beside the edges, and the image reaches 1.006, above the card’s brightest. When the clip ends, a full-width slider, “Cutoff”, sets 0.02 to 0.5 cycles per pixel in steps of 0.01 (arrow keys 0.01, Page Up and Page Down 0.05). The setting is kept in the link as mask.c, and the caption gives the energy kept and the RMS change.

Watch the filtered panel as the circle shrinks: first the stripes vanish, then the edges blur and rings appear beside them.

Try it: drag the cutoff from 0.25 down to 0.22. The energy kept falls from 96.8 % to 89.4 %, and the RMS change grows from 0.039 to 0.070. That energy was the stripes: as the circle passes inside their spots, they fade from the filtered image.

Why a sharp circle rings

Look again at the end frame. At 0.05 the image overshoots to 1.006, brighter than the card’s brightest, 0.900. It also dips to −0.008, below the dark square’s 0.1.

An average with positive weights never leaves the range of the values it averages. So this blur must have some negative weights.

Multiplying spectra is the same as convolving images, by the circular convolution of Properties of the DFT (13.3), “Circular convolution”, now in two directions. So the mask is a kernel, as in 28.2: the inverse DFT of HH. Its weights add up to H[0,0]=1H[0,0]=1, so by “What the weights add up to” in 28.2, flat areas stay flat.

And it does. Take a cutoff of 0.1, and walk along a row from the kernel’s centre. The weights are positive out to 6 pixels, negative from 7 to 11 and positive again from 12 to 16. Those rings of sign are what make rings beside every edge.

This is the Gibbs bump again. In Convergence and the Gibbs phenomenon (7.3), “The bump keeps its height” showed that stopping a Fourier series sharply leaves an overshoot beside a jump. The circle is a sharp stop in two directions.

A mask that falls gently to 0 instead does not ring. Try a Gaussian bell of the distance from the centre, with a spread of 0.05 cycles per pixel. Its kernel has no negative weights, and its blur stays between 0.1 and 0.9.

The rings are not the only side effect. Because the DFT sees the image as a repeating tile, the blur also mixes each border with the opposite one. At 0.05, the pixel in row 10 at the first column comes out 0.417, against 0.250 on the card.

Phase carries the shape

Every bin is a complex number, so it has a size and an angle:

X[k1,k2]=∣X[k1,k2]∣ ej∠X[k1,k2].X[k_1,k_2]=\lvert X[k_1,k_2]\rvert\,e^{j\angle X[k_1,k_2]}.

The size, the magnitude, says how much of that grating the image holds. The angle, the phase, says where its stripes sit.

You can see the second claim with “Shifting on a ring” of 13.3.

A circular shift of a row by n0n_0 samples turns each bin by −2πkn0/N-2\pi kn_0/N and leaves its size alone. The same holds in each direction of an image. So moving a shape around the tile changes only the phases. The magnitudes cannot know where anything is.

How much of the picture is in each? Let’s swap them. I need a second image, B, again from a formula. On a flat ground of 0.4 it has three things:

  • a ring of 0.95, between 14 and 22 pixels from column 88, row 36;
  • a triangle of 0.05, its tip at column 36, row 70, widening by one pixel on each side every two rows down to row 118;
  • three bars of 0.8, each 7 pixels wide, at columns 78 to 84, 92 to 98 and 106 to 112, from row 72 to row 120.

Its values run from 0.05 to 0.95.

Now I build two hybrids. The first takes the card’s magnitudes and B’s phases, and the second takes B’s magnitudes and the card’s phases. Each goes back through the inverse DFT. To say which image a hybrid looks like, I use the correlation coefficient from “A cloud that leans: correlation” in Random variables for signals (24.1), with the pixels as the draws.

Fig. Swap the phases: each hybrid looks like the image whose phase it took. Correlation with B: 0.71 (and with A, −0.03); the other hybrid with A: 0.71. The hybrids leave 0 to 1 (−0.240 to 1.422 and −0.184 to 1.222): each is shown mapped linearly from its own range to black and white.
Describe this picture

Four images of 128 × 128 pixels in grey levels, each named underneath. A is the card: a gradient, a bright disc, a dark square and a patch of diagonal stripes. B is a bright ring, a dark triangle and three bright bars on a grey ground. The third has the magnitudes of A and the phases of B, and in it the ring, the triangle and the bars of B appear. The fourth has the magnitudes of B and the phases of A, and in it the disc, the square and the stripes of the card appear. The hybrids leave 0 to 1 (−0.240 to 1.422 and −0.184 to 1.222), so each is shown mapped linearly from its own range to black and white.

Look at the two hybrids. The hybrid with B’s phases looks like B, and correlates with it at 0.71, though every magnitude in it is the card’s. With the card itself it correlates at −0.03, close to 0. The other hybrid correlates with the card at 0.71. For comparison, the card and B correlate at −0.11.

Why does the phase win? An edge is a place where many gratings line up: at the edge they all rise together. The sine waves that build a square wave do the same at its jump, in “Smooth arrows, flat tops” of Signals as sums of sinusoids (7.1).

The phases hold where each grating sits, so they hold where the gratings line up, and so where the edges are. The magnitudes hold only how strong each grating is, and the card and B have a similar spread of strengths: on average strong near the centre, weaker farther out.

Key idea

The magnitude says how much of each grating an image holds; the phase says where its stripes line up. Edges are places where many gratings line up, so the shapes live in the phase.

Worked example

1. A 2 × 2 image by hand. Take the image with rows 1, 2 and 3, 4: x[0,0]=1x[0,0]=1, x[1,0]=2x[1,0]=2, x[0,1]=3x[0,1]=3, x[1,1]=4x[1,1]=4. For N=2N=2 the probe is e−jπkn=(−1)kne^{-j\pi k n}=(-1)^{kn}, so a 2-point DFT gives the sum and the difference of two numbers.

The rows: 1, 2 give 3 and −1, and 3, 4 give 7 and −1. The columns of that: 3, 7 give 10 and −4, and −1, −1 give −2 and 0. So X[0,0]=10X[0,0]=10, X[0,1]=−4X[0,1]=-4, X[1,0]=−2X[1,0]=-2 and X[1,1]=0X[1,1]=0.

Read them. X[0,0]=10X[0,0]=10 is 222^2 times the average, 2.5. The brightness grows by 1 across and by 2 down, so the across bin is −2 and the down bin is −4, twice as large. X[1,1]=0X[1,1]=0 because nothing in the image changes along the diagonal alone.

2. Separability’s saving. For 64 × 64 pixels, the direct sum needs 644=16 777 21664^4=16\,777\,216 multiplications. Rows then columns needs 2⋅643=524 2882\cdot64^3=524\,288, which is 32 times fewer. With 64-point FFTs it is 128⋅192=24 576128\cdot192=24\,576.

3. Where a grating lands. Take cos⁡(2π(3n1+6n2)/64)\cos\big(2\pi(3n_1+6n_2)/64\big) with height 0.25. Write it as two turning arrows, half each:

0.125 ej2π(3n1+6n2)/64+0.125 e−j2π(3n1+6n2)/64.\begin{aligned} &0.125\,e^{j2\pi(3n_1+6n_2)/64}\\ &\quad+0.125\,e^{-j2\pi(3n_1+6n_2)/64}. \end{aligned}

The first arrow matches the probe (k1,k2)=(3,6)(k_1,k_2)=(3,6) exactly, so every product is 0.1250.125 and the sum is 642⋅0.12564^2\cdot0.125. Divided by 64264^2, that is 0.125 at (3,6)(3,6). The second gives 0.125 at (−3,−6)(-3,-6). Every other probe turns against the arrows and cancels, as in “Rows that cancel” of The DFT as a matrix (13.5).

Where you’ll meet this

Large blurs are done through the 2-D FFT. A 31 × 31 kernel takes 961 multiplications per pixel when it slides. Through the FFT, the cost per pixel grows only with the logarithm of the image’s size, whatever the kernel.

Phase correlation lines up two photos of the same scene. Divide one spectrum by the other, keep only the phase differences, and take the inverse DFT. A single sharp peak appears at the shift between them. It is used to align the frames of a shaky video and the tiles of a panorama.

An MRI scanner measures the 2-D spectrum of a slice of the body directly, one line of bins at a time. The picture you see is its inverse 2-D DFT.

JPEG uses a close cousin, the 2-D DCT, on blocks of 8 × 8 pixels. It is separable in the same way, and Image compression (28.4) builds it up.

The maths behind it · matrix products

Write the image as a matrix X\mathbf{X}, with row n2n_2 of the matrix holding row n2n_2 of the image. Its 2-D DFT is FXF\mathbf{F}\mathbf{X}\mathbf{F}, with the DFT matrix F\mathbf{F} of 13.5. Multiplying by F\mathbf{F} on the right transforms every row, and on the left every column: that is the separability.

The maths behind it · autocorrelation

The squared magnitudes of the 2-D DFT are the DFT of the image’s 2-D autocorrelation, as in Power spectral density (24.4), now in two directions. Texture statistics start from there.

Reference card

QuantityFormulaNotes
2-D DFT∑x[n1,n2]e−j2π(k1n1+k2n2)/N\sum x[n_1,n_2]e^{-j2\pi(k_1n_1+k_2n_2)/N}one bin for each grating
Separablerows, then columns2N2N 1-D DFTs of NN points
CostN4N^4 direct, 2N32N^3 by rows and columnsFFTs: N2log⁡2NN^2\log_2N
CentreX[0,0]/N2X[0,0]/N^2the average brightness
A gratingtwo dots, at ±(k1,k2)\pm(k_1,k_2)X[−k1,−k2]=X[k1,k2]∗X[-k_1,-k_2]=X[k_1,k_2]^*
An edgea line through the centreat right angles to the edge
MaskingY=H⋅XY=H\cdot Xsharp masks ring
Magnitudehow much of each gratingblind to where things are
Phasewhere the stripes line upcarries the shapes

End of lesson 28.3

Where to go next.

Phasorium
LibraryEvery lesson, in order

Parts

About Phasorium
Look