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
Reference¶
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:
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)
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.
## Data
def get_wavelength_adjusted_vec(dx,dy):
# 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)
_dx,_dy,wavelength = get_wavelength_adjusted_vec(dx,dy)
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,
color="red",head_width=0.5, head_length=0.3)
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.
# 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¶
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)
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)$.
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()
Inverse Fourier Transform¶
Using $F$ inferse Fourier function, we can transform back the FFT to get the original sinosoid image.
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
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 = fig.add_subplot(1,2,1)
ax.imshow(np.abs(fft),cmap="gray")
ax.set_title("FFT")
ax.set_xlabel("$u$",fontsize=30)
ax.set_ylabel("$v$",fontsize=30)
ax = fig.add_subplot(1,2,2)
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)
_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
_dx,_dy,wavelength = get_wavelength_adjusted_vec(dx,dy)
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))
ax = fig.add_subplot(1,2,1)
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,
color="red",head_width=0.5, head_length=0.3)
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 = fig.add_subplot(1,2,2)
im = ax.imshow(np.abs(fft_shift), cmap="gray",extent=extent)
fig.colorbar(im,ax=ax)
ax.set_title("Fourier Coefficient |f(u,v)|")
ax.set_xlabel("$u$",fontsize=30)
ax.set_ylabel("$v$",fontsize=30)
plot_nonzero_Fourier_coefficient_line(_y,_x,_N,_M)
plt.show()
