In [79]:
def iacos(s, t, theta, phi):
r, c = np.indices(s)
tc = t / np.cos(theta)
tr = t / np.sin(theta)
f = np.cos(2*np.pi*(r/tr + c/tc) + phi)
return f
def iatile(f, new_size):
f = np.array(f)
if len(f.shape) == 1: f = f[newaxis,:]
aux = np.resize(f, (new_size[0], f.shape[1]))
aux = np.transpose(aux)
aux = np.resize(aux, (new_size[1], new_size[0]))
g = np.transpose(aux)
return g
def iaidft(F):
s = F.shape
if F.ndim == 1: F = F[np.newaxis,np.newaxis,:]
if F.ndim == 2: F = F[np.newaxis,:,:]
p,m,n = F.shape
A = ia.dftmatrix(m)
B = ia.dftmatrix(n)
C = ia.dftmatrix(p)
Faux = np.conjugate(A).dot(F)
Faux = Faux.dot(np.conjugate(B))
f = np.conjugate(C).dot(Faux)/(np.sqrt(p)*np.sqrt(m)*np.sqrt(n))
return f.reshape(s)
This demonstration illustrates the computation of the DFT - Discrete Fourier Transform of a particular discrete cosine wave. We realized that depending on the wave length and angle chosen, its DFT is not a pair of complex conjugate points. The demonstration shows the reason of this by calling the periodic nature of the DFT and when four discrete cosine waves are tiled together, it becomes evident that the wave is not an ideal discrete cosine wave. The demonstration goes on by computing the nearest normalized frequencies and and computes a new discrete cosine wave from the inverse DFT of a single pair of complex conjugate placed at the nearest normalized frequencies. In this case, when four images area tiled together, we see there are no abrupt changes in the boundaries of the images just tiled.
A discrete cosine wave of an image with size 256 columns and 128 rows, with wave length of 100 pixels at 45 degrees is computed and shown below.
In [80]:
import numpy as np
import sys,os
ia898path = os.path.abspath('/home/lotufo')
if ia898path not in sys.path:
sys.path.append(ia898path)
import ia898.src as ia
from math import pi, sin, cos, atan
%matplotlib inline
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
import matplotlib.colors as pltc
In [81]:
s = (128,256)
t = 100
theta = 45 * pi/180
f = iacos(s, t, theta, 0)
ia.adshow(ia.normalize(f), title='f - cosine wave')
The computation of its DFT - Discrete Fourier Transform reveals that the result is not a single pair of complex conjugate points, but also a number of points in the vertical and horizontal directions. This behavior has to do with abrupt changes in the image. As both the DFT and the image are periodic, we place four images as tiles. We can see that when these images are tiled together, there are many abrupt changes in the images as they do not fit together as a single cosine wave.
In [82]:
F = ia.dft(f)
ia.adshow(ia.dftview(F), title='F - DFT of f')
ia.adshow(iatile(ia.normalize(f),2 * np.array(np.shape(f))),title='f tiled 4 times')
From the wave length 100 and angle 45 degrees, we compute the value of the normalized horizontal and vertical frequencies and and noticed that they are not integers and so there is no correspondence to an ideal cosine wave in the discrete Fourier domain:
In [83]:
u = s[1] * cos(theta) / t
v = s[0] * sin(theta) / t
print('u,v = ', u,v)
The nearest integer u and v are 2 and 1, respectively. We create an synthetic Discrete Fourier image with a single pair of complex conjugate points at this nearest integer u and v.
In [84]:
ui, vi = round(u), round(v)
print('ui,vi = ', ui, vi)
FS = np.zeros(s)
FS[vi,ui] = 1
FS[s[0]-vi, s[1]-ui] = 1
print('FS is complex conjugate? ',ia.isccsym(FS))
ia.adshow(ia.dftview(FS),title='DFT of a single pair of complex conjugate points')
Computing the inverse DFT we can find the correspondent cosine wave. We observe now that tiling 4 images together, there are no abrupt changes in the boundaries. We compute the equivalent wave length and angle, and we notice that the nearest ideal cosine wave has now a wavelength of 90.5 pixels instead of 100, in the same direction (45 degrees).
In [86]:
fs = np.real(iaidft(FS))
ia.adshow(ia.normalize(fs), title='fs - ideal cosine wave')
ia.adshow(iatile(ia.normalize(fs), 2 * np.array(np.shape(fs))), title='fs tiled 4 times')
txi,tyi = s[1]/ui, s[0]/vi
thetai = atan(tyi/txi)
ti = txi * cos(thetai)
print('new period and angle = ', ti, thetai * 180/pi)
We place both the original cosine wave and the ideal cosine wave together for comparison.
In [88]:
ia.adshow(ia.normalize(f), title='f - cosine wave with period 100 pixels')
ia.adshow(ia.normalize(fs), title='fs - ideal cosine wave with period 90.5 pixels')