Frequency analysis of images from scratch

This blog reviews frequency analysis on images. If you follow this blog, you will understand

• The two-dimensional discrete Fourier transform
• How to calculate wavelength of the Sinosoid
• What exactly np.fft.fft2 and np.fft.fftshift are doing.

This blog series on frequency analysis on images will continue Low and High pass filtering experiments

Sinosoid example¶

Troughout this blog, I will consider following Sinosoid as an example. $$g(v,u) = cos\Big(-2 \pi \left( u x + v y\right)\Big)$$

Here,

• $x$: frequency along horizontal axis. The larger it is, the shorter one period is. If $x=0$, then the cos wave has no change along horizontal direction.
• $y$: frequency along vertical axis.

What is this wavelength of this sinusoid?¶

The wavelength is the spatial period of a periodic wave—the distance over which the wave's shape repeats.

As g is a function of two variables, to talk about the periodicity, we need to have the direction along witch the periodicity shall be calculated. I will re define this function g in vector format. Let $\boldsymbol{u}=(u,v), \boldsymbol{x} = (x,y)$,

$$g(\boldsymbol{u}) = cos\Big(-2 \pi \boldsymbol{u}^T \boldsymbol{x} \Big) --(1)$$

Let $\boldsymbol{e}$ be a normalized vector along which the periodicity should be calculated. To calculate the periodicity, we need to find the wavelength $w$ that satisfy: $$g(\boldsymbol{u}+\boldsymbol{e}w) = g(\boldsymbol{u})$$

For example, if $\boldsymbol{e}=(1,0)$, horizontal axis,

\begin{array}{rl} g(\boldsymbol{u}+\boldsymbol{e}w) &= cos\Big( -2 \pi \big( \left(u+w\big) x + v y\right)\Big)\\ &= cos\Big(-\big( 2 \pi u x + 2 \pi v y + 2 \pi x w\big)\Big) -- (2) \end{array}

For which $w$ do (1) and (2) become the same? (1) and (2) are the same when $2 \pi x w = 2\pi n, n=1,2,3,...$. Especially, the wavelength is shortest when $n = 1$. This means that the wavelength along x-axis is: $$w = \frac{1}{x}$$ The wavelength along vertical axis can also be found and it is $w=1/y$. These observations give some suggestions to the choice of $x$ and $y$. When $x$ and

Now consider a normalized vector along vector $\boldsymbol{x}$, $\boldsymbol{e}=\boldsymbol{x}/||\boldsymbol{x}||$. \begin{array}{rl} g(\boldsymbol{u}+\boldsymbol{e}w) &= cos\Big( -2 \pi (\boldsymbol{u}+\boldsymbol{e}w)^T \boldsymbol{x}\Big)\\ &= cos\Big( -2 \pi \big(\boldsymbol{u}+\frac{\boldsymbol{x}}{||\boldsymbol{x}||}w\big)^T \boldsymbol{x}\Big)\\ &= cos\Big( -2 \pi \big(\boldsymbol{u}^T \boldsymbol{x} +\frac{w \boldsymbol{x}^T \boldsymbol{x} }{||\boldsymbol{x}||}\big) \Big)\\ &= cos\Big( -2 \pi \boldsymbol{u}^T \boldsymbol{x}-2 \pi w ||\boldsymbol{x}|| \Big)--(3)\\ \end{array}

(3) becomes the same as (1) when $2 \pi w ||\boldsymbol{x}||= 2\pi n, n=1,2,3,...$. Especially, the wavelength is shortest when $n = 1$. This means that the wavelength along $\boldsymbol{x}$ is: $$w = \frac{1}{||\boldsymbol{x}||} = \frac{1}{\sqrt{x^2 + y^2}}$$

OK. That is enough theory for now. Let's visualize this sinusoid with various $(x,y)$. First, define our sinusoid:

In [1]:
import numpy as np
import matplotlib.pyplot as plt
def myfun(v, u, dy, dx ):
'''
v  : (often) integer value
u  : (often) integer value
dy : frequency along vertical axis
dx : frequency along horizontal axis

cos(-2pi * (ux + vy)) = 1
when - 2 pi * (ux + vy) = 2 * pi *n  (n=0,1,..)
i.e.,          ux + vy = - n
this function taks its maximum value 1
'''
term = -2. * np.pi *( u * dx + v * dy )
return(np.cos(term))

def get_sinusoid_img(dy,dx,_M,_N):
orig_img = np.zeros((_N,_M))
startx = - _M/2
starty = - _N/2
countx = -1
for x in range(startx,startx + _M):
countx += 1
county = -1
for y in range(starty,starty + _N):
county += 1
orig_img[county,countx] = myfun(y,x,dy,dx)
return(orig_img)
/Users/yumikondo/Documents/blog/blog/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

Now visualize the sinusoid for (y,x):

• (0/18,10/42)
• (3/18,0/42)
• (7/18,1/42)
• (5/18,4/42)

Notice that the vector of length = wavelength go through one cycle.

In [2]:
## Data
# this calculation is wrong....
#if dx > 0.5:
#    dx_dist = 1 - dx
#else:
#    dx_dist = dx
#if dy > 0.5:
#    dy_dist = 1 - dy
#else:
#    dy_dist = dy
dx_dist = dx
dy_dist = dy
xynorm = np.sqrt(dx**2 + dy **2)
wavelength = 1 / np.sqrt(dx_dist**2 + dy_dist **2)
_dx, _dy = dx/xynorm *wavelength, dy/xynorm * wavelength
return(_dx,_dy,wavelength)

_N = 18 ## the number of pixels along y axis
_M = 42 ## the number of pixels along x axis
_ys = [0,3, 7,5]
_xs = [10,0,1,4]
for _y, _x in zip(_ys, _xs):
dy = _y / float(_N)
dx = _x / float(_M)
orig_img = get_sinusoid_img(dy,dx,_M,_N)

const = 0.5

fig = plt.figure(figsize=(_M*const,_N*const))
plt.imshow(orig_img,cmap="gray")
plt.arrow(x=3, y=3,
dx=_dx, dy=_dy,
plt.annotate(s="vector x (wavelength={:4.3f})".format(wavelength),
xy=(3,3),xytext=((3 + _dx)/2.0,(3 + _dy)/2.0),fontsize=20,color="red")
plt.title("$cos (-2\pi ( {}v/{} + {}u/{} ) )$".format(_y,_N,_x,_M),fontsize=30)
plt.xlabel("$u$",fontsize=30)
plt.ylabel("$v$",fontsize=30)
plt.colorbar()
plt.show()

Frequency analysis¶

Now, I will understand this sinosoid's frequency via the frequency analysis. Previously, I reviewed 1 dimensional discrete fourier transform. We can easily extent this to 2 dimensional setting. The two-dimensional discrete Fourier transform can be expressed with two equations:

\begin{array}{rl} F(v,u) &= \sum_{x= - M/2 }^{M/2 - 1} \sum_{y=-N/2}^{N/2 - 1} f(y,x) \exp\Big( 2 \pi j \left( \frac{x u }{M} + \frac{y v }{N}\right)\Big)\\ f(y,x) &= \frac{1}{M N}\sum_{u=- M/2 }^{M/2 - 1} \sum_{v=-N/2}^{N/2 - 1} F(v,u) \exp\Big( - 2 \pi j \left( \frac{x u }{M} + \frac{y v }{N}\right)\Big) \end{array}

Notice that the summation starts from $- M/2 + 1$ rather than 0, and run through $M$ times. You can also define the two-dimensional discrete fourier transform by summing from 0 to $M+1$. However, in image signal processing, it is more common to define the sum from $- M/2 + 1$. So that the 0 frequency lies at the center.

Now, let's focus on one senosoid: $$g(v,u) = cos\left(-2 \pi \left( x^* u + y^*v\right)\right)$$ where $$x^* = \frac{x'}{M}, y^* = \frac{y'}{N}$$ We will use these rational numbers with $M$ and $N$ as their denominators to ease the calculation.

In [3]:
# focus on one senosoid
_N = 18 ## the number of pixels along y axis
_M = 42 ## the number of pixels along x axis
_y = 5
_x = 14
dy = _y / float(_N)
dx = _x / float(_M)
orig_img = get_sinusoid_img(dy,dx,_M,_N)

Which Fourier coefficients $f(y,x)$ have nonzero values?¶

You can mathematically show that there are two sets of frequencies $(y,x)$ with nonzero Fourier coefficients within $y =- \frac{N}{2},- \frac{N}{2}+1,\cdots,0,1,\cdots,\frac{N}{2}-2,\frac{N}{2}-1$, and $x =- \frac{M}{2},- \frac{M}{2}+1,\cdots,0,1,\cdots,\frac{M}{2}-2,\frac{M}{2}-1$. The following lengthy math shows what these two sets exactly are.

First, our function $g$ can be expressed by two exponents:

\begin{array}{rl} cos\Big( - 2\pi \big(x^*u + y^*v\big) \Big) = \frac{1}{2} \Big( \exp\big( - 2\pi j \left(x^*u + y^*v\right) \big) + \exp\big( 2\pi j \left(x^*u + y^*v\right) \big) \Big) \end{array}

Therefore the Fouier coefficients can be calculated as, \begin{array}{rl} f(y,x) &= \frac{1}{M N}\sum_{u=- M/2 }^{M/2 - 1} \sum_{v=-N/2}^{N/2 -1} F(v,u) \exp\Big( - 2 \pi j \left( \frac{x u }{M} + \frac{y v }{N}\right)\Big)\\ &= \frac{1}{M N}\sum_{u=- M/2 }^{M/2 - 1} \sum_{v=-N/2}^{N/2 -1} \frac{1}{2} \Big( \exp\big( - 2\pi j \left(x^*u + y^*v\right) \big) + \exp\big( 2\pi j \left(x^*u + y^*v\right) \big) \Big)\exp\Big( - 2 \pi j \left( \frac{x u }{M} + \frac{y v }{N}\right)\Big)\\ &= \frac{1}{2}\frac{1}{M N}\sum_{u=- M/2 }^{M/2 - 1} \sum_{v=-N/2}^{N/2 -1} \exp\big( - 2\pi j \big( \left(x^*u + y^*v\right) + \left( \frac{x u }{M} + \frac{y v }{N}\right)\big)\big) + \exp\big( 2\pi j \big( \left(x^*u + y^*v\right) - \left( \frac{x u }{M} + \frac{y v }{N}\right)\big)\big)\\ &= \frac{1}{2}\frac{1}{M N}\sum_{u=- M/2 }^{M/2 - 1} \sum_{v=-N/2}^{N/2 -1} \exp\Big( - 2\pi j \Big( \left( u\left[ x^* + \frac{x }{M}\right] + v\left[y^* + \frac{y }{N}\right]\right)\Big) \Big)+ \exp\Big( 2\pi j \Big( \left( u\left[ x^* -\frac{x }{M}\right] + v\left[y^* -\frac{y }{N}\right]\right)\Big)\Big)\\ &= \frac{1}{2}\frac{1}{M N}\sum_{u=- M/2 }^{M/2 - 1} \sum_{v=-N/2}^{N/2 -1} \exp\Big( - 2\pi j \Big( \left( u\left[ \frac{x'}{M} + \frac{x }{M}\right] + v\left[ \frac{y'}{N}+ \frac{y }{N}\right]\right)\Big)\Big) --(4)\\ & \:+ \frac{1}{2}\frac{1}{M N}\sum_{u=- M/2 }^{M/2 - 1} \sum_{v=-N/2}^{N/2 -1} \exp\Big( 2\pi j \Big( \left( u\left[ \frac{x'}{M} -\frac{x }{M}\right] + v\left[ \frac{y'}{N}-\frac{y }{N}\right]\right)\Big)\Big)--(5) \\ \end{array}

Consider the term (4)¶

Consider the term within the summation of (4). $$\exp\Big( - 2\pi j \Big( u\left[ \frac{x'}{M} + \frac{x }{M}\right] + v\left[ \frac{y'}{N}+ \frac{y }{N}\right] \Big)\Big) = Cos\left(\theta\right) + j Sin\left( \theta \right)$$ where $$\theta = -2\pi\left( u\left[ \frac{x'}{M} + \frac{x }{M}\right] + v\left[ \frac{y'}{N}+ \frac{y }{N}\right]\right)$$ This term becomes 1 when $\theta=2\pi k, k =...,-3,-2,-1,0,1,2,3,...$. What values of $(y,x)$ is $\theta$ multiplicative of $2\pi$? Here, I will restrict the focus on $y =- \frac{N}{2},- \frac{N}{2}+1,\cdots,0,1,...,\frac{N}{2}-2,\frac{N}{2}-1$, and $x =- \frac{M}{2},- \frac{M}{2}+1,\cdots,0,1,...,\frac{M}{2}-2,\frac{M}{2}-1$.

• Clearly, $\theta = 0$ when $x=x'$ and $y=y'$
• $\theta$ is multiplicative of $2\pi$ when $x=M - x'$ and $y = M - y'$.

\begin{array}{rl} \theta &= -2\pi\left( u\left[ \frac{x'}{M} + \frac{x }{M}\right] + v\left[ \frac{y'}{N}+ \frac{y }{N}\right]\right)\\ &= -2\pi\left( u\frac{1}{M}\left[x' + (M-x')\right] + v\frac{1}{N}\left[ y'+ (N - y')\right]\right)\\ &= -2\pi\left( u + v\right)\\ \end{array} for any integer $u,v$.

Therefore, when $(y,x)=(-x',-y')$ or when $(y,x)=(M-x',N-y')$, (4)=1/2.

What would happen when $(y,x)\ne (-x',-y')$ or when $(y,x)\ne(M-x',N-y')$?¶

\begin{array}{rl} &\sum_{u=- M/2 }^{M/2 - 1} \sum_{v=-N/2}^{N/2 -1} \exp\Big( - 2\pi j\left( u\left[ \frac{x'}{M} + \frac{x }{M}\right] + v\left[ \frac{y'}{N}+ \frac{y }{N}\right]\right)\Big)\\ =& \sum_{u=- M/2 }^{M/2 - 1} \exp\Big( - \frac{2\pi j}{M} \left[ x' + x \right]\Big)^u \sum_{v=-N/2+1}^{N/2} \exp\Big(-\frac{2\pi j}{N}\left[ y' + y \right]\Big)^v\\ \end{array}

Notice that when $a < b$, $$S_{a,b} = r^a + r^{a+1} + r^{a+2} +... + r^{b-1} + r^b$$ Then when $r\ne 1$ $$S_{a,b} = \frac{r^{a} - r^{b+1}}{1 - r} = \frac{1 - r^{b+1-a}}{r^{-a}(1-r )}$$ Therefore, using this formula, setting $r=exp\big(-\frac{2\pi j}{M} [x'+x]\big)$, $a=-M/2$ and $b=M/2-1$, we can show that: $$\sum_{u=- M/2 }^{M/2 -1} \exp\Big( - \frac{2\pi j}{M} \left[ x' + x \right]\Big)^u= 0.$$ when $r \ne 1$ that is, when $x=x'$ and $x = M - x'$.

Also, using this formula, setting $r=exp\big(-\frac{2\pi j}{N}\left[ y' + y \right]\big)$, $a=-N/2$ and $b=N/2-1$, we can show that: $$\sum_{v=-N/2}^{N/2-1} \exp\Big(-\frac{2\pi j}{N}\left[ y' + y \right]\Big)^v=0$$ when $r \ne 1$ that is, when $y=y'$ and $y = N - y'$. So (4) = 0.

Consider the term (5)¶

Similar to the arguments used in term (4), $(x,y)=(x',y')$ and $(x,y)=(-M + x',-N + y')$ yield nonzero (5). And otherwise (5) = 0.

Therefore, we expect that $f(-x',-y')$,$f(M-x',N-y')$, $f(x',y')$ and $f(-M+x',-N+y')$ and $f(x',y')$ are the nonzero Fourier coefficients and otherwise Fourier coefficients are zero. There are two sets of Fourier coefficients that will lie within the range of $y =- \frac{N}{2},- \frac{N}{2}+1,\cdots,0,1,...,\frac{N}{2}-2,\frac{N}{2}-1$, and $x =- \frac{M}{2},- \frac{M}{2}+1,\cdots,0,1,...,\frac{M}{2}-2,\frac{M}{2}-1$.¶

Now, let's check if this is indeed the case!

Define $F$ and $f$ functions¶

In [4]:
def F(v,u,fs): ## y = v (N), x = u (M)
_F = 0
N, M = fs.shape
startu = - fs.shape[1]/2
startv = - fs.shape[0]/2
xcount = -1
for x in range(startu, startu + fs.shape[1]):
xcount += 1
ycount  = -1
for y in range(startv, startv + fs.shape[0]):
ycount += 1
term = x*u/float(M) + y*v/float(N)
_F += fs[ycount,xcount]*np.exp(2.0*np.pi*1j*term)
return(_F)
def f(y,x,Fs):
_f = 0
N, M = Fs.shape
startx = - Fs.shape[1]/2
starty = - Fs.shape[0]/2
ucount = - 1
for u in range(startx, startx + Fs.shape[1]):
ucount += 1
vcount  = -1
#print("u={:2.0f} count={:2.0f}".format(u,ucount))
for v in range(starty, starty + Fs.shape[0]):
vcount += 1
term = x*u/float(M) + y*v/float(N)
_f += Fs[vcount,ucount]*np.exp(-2.0*np.pi*1j*term)
_f /= np.prod(Fs.shape)
return(_f)
In [5]:
def plot_nonzero_Fourier_coefficient_line(_y,_x,_N,_M):
## horizontal and vertical lines along nonzero Fourier coefficient
if _x < _M/2: ## x is always positive
plt.axvline(_x)
plt.axvline(-_x)
else: ## x is not between - M/2,..., M/2-1
plt.axvline(_M - _x)
plt.axvline(-_M + _x)
if _y < _N/2:
plt.axhline(_y)
plt.axhline(-_y)
else: ## y is not between - N/2,..., N/2-1
plt.axhline(_N - _y)
plt.axhline(-_N + _y)

Discrete Fourier Transform on Sinusoidal Waves¶

You see that the discrete Fourier transform yields nonzero frequency coefficients at $f(5,4)$ and $f(-5,-4)$.

In [6]:
ySample, xSample = _N, _M
fs = np.zeros((ySample,xSample),dtype=complex)

startx = -orig_img.shape[1]/2
starty = -orig_img.shape[0]/2

countx = -1
xs = range(startx,startx + orig_img.shape[1])
ys = range(starty,starty + orig_img.shape[0])
for x in xs:
countx += 1
county = -1
for y in ys:
county += 1
fs[county,countx] = f(y,x,orig_img)
if np.abs(fs[county,countx]) > 0.0001:
print("x={:2.0f}, y={:2.0f}, f={:+5.3f}".format(x,y,fs[county,countx]))

## adjust the axis: fs[0,0]  corresponds to the frequency (starty, startx)
##                  fs[0,-1] corresponds to the frequency (starty, startx + orig_img.shape[1])
##                  fs[-1,0] corresponds to the frequency (starty + orig_img.shape[0], starty)
##                  fs[0,-1] corresponds to the frequency (starty + orig_img.shape[0], startx + orig_img.shape[1])
extent = (startx - 0.5, ## - 0.5 is necessary. See https://matplotlib.org/api/_as_gen/matplotlib.pyplot.imshow.html
startx - 0.5 + orig_img.shape[1],
starty - 0.5 + orig_img.shape[0],
starty - 0.5)

plt.figure(figsize=(20,10))
plt.imshow(np.abs(fs),cmap="gray",extent=extent) # left, right, bottom, top)
plot_nonzero_Fourier_coefficient_line(_y,_x,_N,_M)
plt.xticks(xs)
plt.yticks(ys)
plt.xlabel("$u$",fontsize=30)
plt.ylabel("$v$",fontsize=30)
plt.title("DFT shape={}".format(fs.shape),fontsize=20)
plt.colorbar()
plt.show()
x=-14, y=-5, f=+0.500+0.000j
x=14, y= 5, f=+0.500-0.000j

Inverse Fourier Transform¶

Using $F$ inferse Fourier function, we can transform back the FFT to get the original sinosoid image.

In [7]:
startx = - orig_img.shape[1]/2
starty = - orig_img.shape[0]/2
reco_img = np.zeros((ySample,xSample),dtype=complex)
countx = -1
for x in range(startx,startx + orig_img.shape[1]):
countx += 1
county  = -1
for y in range(starty, starty + orig_img.shape[0]):
county += 1
reco_img[county,countx] = F(y,x,fs)
plt.imshow(np.real(reco_img),cmap="gray")
plt.title("reconstructed image")

plt.show()

Understand np.fft.fftshift and np.fft.fft2¶

In practice you can do the calculation using Discrete Fourier Transform function in numpy

In [8]:
fft = np.fft.fft2(orig_img)
# fftshift
# If X is a matrix, then fftshift swaps the first quadrant of X with the third,
# and the second quadrant with the fourth.
fft_shift = np.fft.fftshift(fft)
fig = plt.figure(figsize=(20,10))
ax.imshow(np.abs(fft),cmap="gray")
ax.set_title("FFT")
ax.set_xlabel("$u$",fontsize=30)
ax.set_ylabel("$v$",fontsize=30)

ax.imshow(np.abs(fft_shift),cmap="gray",
extent=extent) # (left, right, bottom, top
ax.set_title("FFT f(0,0) at the center")
ax.set_xlabel("$u$",fontsize=30)
ax.set_ylabel("$v$",fontsize=30)
plot_nonzero_Fourier_coefficient_line(_y,_x,_N,_M)
plt.show()

Here I show you several example Sinosoid and its corresponding Frequency plot using np.fft.fft2. As we used simple Sinosoid, you can see where exactly the nonzero fourier coefficients are. (See blue lines)

In [9]:
_M = 32
_N = 32
xs = range(0,_M/2, 7)
ys = range(0,_N/2, 7)
startx = -_M/2.0
starty = -_N/2.0
extent = [startx - 0.5,      # Left
startx + _N - 0.5, # Right
starty +_M - 0.5,  # Bottom
starty - 0.5]      # Top
for _y in ys:
for _x in xs:
dy = _y/float(_M) ## v
dx = _x/float(_N) ## u
orig_img = get_sinusoid_img(dy,dx,_M,_N)
fft = np.fft.fft2(orig_img)
fft_shift = np.fft.fftshift(fft)

fig = plt.figure(figsize=(18,7))
im = ax.imshow(orig_img,cmap="gray")
fig.colorbar(im,ax=ax)
ax.set_title("$Cos(-2\pi(u{}/{} + v{}/{})$".format(_x,_N,_y,_M)) # -2. * np.pi *( u * dx + v * dy )
ax.set_xlabel("$u$",fontsize=30)
ax.set_ylabel("$v$",fontsize=30)
if ~np.isnan(_dx):
ax.arrow(x=3, y=3, dx=_dx, dy=_dy,
ax.annotate(s="vector x (wavelength={:4.3f})".format(wavelength),
xy=(3,3),xytext=((3 + _dx)/2.0,(3 + _dy)/2.0),fontsize=20,color="red")

ax.set_xlabel("$u$",fontsize=30)
ax.set_ylabel("$v$",fontsize=30)