| SOCR ≫ | TCIU Website ≫ | TCIU GitHub ≫ |
1 Figure 3.1: 1D Waves and 2D Waves
Waves of constant wavelength propagate maintaining its shape but shifting its location like in the image below. However, this oversimplified representation does not work well for computing, as adding wave amplitudes does not work as expected. It simply shifts the wave amplitude up and doesn’t move the wave in the direction of motion. Similar problem arises with wave multiplication by a constant. Wave multiplication by a number increases the height of the peaks and valleys without moving them in space.
These problems are solved by an alternative wave representation where the wave has a complex amplitude, and complex multiplication by a number works well.
The wave phase represents the direction of the wave amplitude in a complex plane. Differences of phases yield interference phenomena. The interference of two waves of amplitude \(i\) will produce an increased amplitude of \(2i=i + i\). When the second wave is of amplitude \(-i\), then their interference destroys the aggregate wave, i.e., produces a zero amplitude wave \(i - i = 0\).
Let’s start with a wave of constant wavelength that represents a particle of definite momentum. For simplicity, we’ll consider just a one space dimension representing the horizontal axis, i.e., the spatial extension of the wave. The remaining two dimensions will represent the phase of the wave complex amplitude. The vertical axis represents the imaginary \(+i\) and \(-i\) directions, and depth represents the real \(+1\) and \(-1\) directions. In this representation, the wave magnitude is constant, however, the 2D phase determines the uniform distance advancement in space.
As the wave propagates, the amplitude cycles through all directions in the 2D complex plane perpendicular to the 1D space dimension (x-axis). The result is a wave function represented by a helical (“corkscrew”) shape. In the previous example, the illusion that the wave amplitude periodically vanishes is just an artifact of an incomplete representation of the wave.
The Complex wave representation clarifies this illusion arising by slicing through the 3D “corkscrew” wave shape and showing just the 1D values in the space of complex amplitudes.
This complex wave representation also addresses the above multiplication problem. Multiplying the 3D wave curve by \(i\) will relocate the amplitudes yielding spatial propagation of the wave, which explains the basic rule of time evolution of quantum theory. The effect of multiplying the wave amplitude by \(i\) at each point is effectively a rotation of the amplitude by one quarter turn around the space of values. Thus, the combination of these point-wise effects results in the wave advancing forward in space (propagation).
t <- seq(-6*pi, 6*pi, length.out=400)
r <- 3
x <- r * cos(t)
y <- r * sin(t)
z <- 2*t
c <- t%%(2*pi)
data1 <- data.frame(x, y, z)
p <- plot_ly(data1, x = ~x, y = ~y, z = ~z, type = 'scatter3d', mode = 'lines', showlegend = F,
line = list(width = 12, color = ~c, colorscale = list(c(0,'red'), c(1,'blue')))) %>%
# trace the z-axis
add_trace(data1, x = 0, y =0, z = ~z, type="scatter3d", mode="lines",
line = list(width = 10, color = ~c, colorscale = list(c(0,'red'), c(1,'blue'))),
name="Z", hoverinfo="none") %>%
# add projection on plane x=-3
add_trace(data1, x = -3, y = ~y, z = ~z, type="scatter3d", mode="lines",
line = list(width = 2, dash="solid", color = ~c,
colorscale = list(c(0,'gray'), c(1,'black'))),
name="Z", hoverinfo="none") %>%
# add projection on plane y=-3
add_trace(data1, x = ~x, y = -3, z = ~z, type="scatter3d", mode="lines",
line = list(width = 2, dash="solid", color = ~c,
colorscale = list(c(0,'gray'), c(1,'black'))),
name="Z", hoverinfo="none") %>%
# add a z-line at x=-3,y=0
add_trace(data1, x = -3, y = 0, z = ~z, type="scatter3d", mode="lines",
line = list(width = 5, dash="solid", color = "gray"),
name="Z", hoverinfo="none") %>%
# add a z-line at x=0,y=-3
add_trace(data1, x = 0, y = -3, z = ~z, type="scatter3d", mode="lines",
line = list(width = 5, dash="solid", color = "gray"),
name="Z", hoverinfo="none") %>%
# add a x-line at y=0, z=-40
add_trace(data1, x = ~x, y = 0, z = -40, type="scatter3d", mode="lines",
line = list(width = 15, dash="solid", color = "gray"),
name="Z", hoverinfo="none")
# add a few annotations/arrows
pRecall that for a single particle, the Hamiltonian operator corresponding to the total energy of the system is defined by:
\[\hat {H}={\underbrace{\frac{\hat {p}^{2}}{2m}}_{\text{particle kinetic energy}}}+ {\underbrace{\frac {1}{2} k \hat {x}^2}_{\text{particle potential energy}}}= {\frac {{\hat {p}}^{2}}{2m}}+{\frac {1}{2}}m\omega ^{2}{\hat {x}}^{2} ,\]
where \(m\) is the particle mass, \(k\) is the force constant, \(w =\sqrt {\frac {k}{m}}\) is the angular frequency of the oscillator, \(\hat {x}\) is the position operator (corresponding to the spatial location \(x\)), and \(\hat {p}\) is the momentum operator (corresponding to the momentum), which is expressed by \(\hat {p}= -i\hbar \frac{\partial}{\partial x}\). Then, the time-independent Schrodinger equation can be expressed as:
\[\hat {H}\left|\psi \right\rangle =E\left|\psi \right\rangle ,\]
where the eigenvalue \(E\) denotes a real number specifying a time-independent energy level and the wavefunction \(\psi\) represents the solution corresponding to that level’s energy eigenstate.
Another way of visualizing the propagation of a wavefunction over a 1D space in terms of its complex amplitude (\(\psi =Ae^{i(kx-wt)}\)), real and imaginary parts, where \(t\) is time, \(w = kc\), wave’s angular frequency which equals \(w=2\pi/T\), where \(T\) is the period of the wave, the position \(x\), and momentum \(p\) are conjugate variables. The wavefunction magnitude represents the probability of finding the particle at a given point \(x\in R^1\) is spread out like a waveform. There is no unique definite position of the particle. As the amplitude oscillates (increases and decreases), the result is a wave with an alternating amplitude.
\(k\) is the wave’s radian spatial frequency (angular wave number) measures how rapidly the disturbance changes over a given distance at a particular point in time, \(k=2\pi/\lambda\), where \(\lambda\) is the wavelength of the wave, and measured in radians per unit distance. Note that \((x,y,z)\), are not part of the equation because the wave’s magnitude and phase are the same at every point.
t <- seq(-6*pi, 6*pi, length.out=400)
r <- 3
x <- r * cos(t)
y <- r * sin(t)
z <- 2*t
c <- t%%(2*pi)
data1 <- data.frame(x, y, z)
scene = list(camera = list(eye = list(x = -1.25, y = 1.25, z = 1.25)),
aspectmode = "manual", aspectratio = list(x=1, y=1, z=5))
# define wave opacity as a parabola
const <- 1/(z[400]^2)
wave_opacity <- 0.3 -(0.7*const)*(z-z[1])*(z-z[400])
# plot(wave_opacity, type = "l")
plot_ly(x = ~z, y = ~wave_opacity, type = 'scatter', mode = 'lines')# bgcolor = toRGB("white", alpha = 0),
# cols = sample(grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)], 400)
# cols_func <- colorRampPalette(c("blue", "red")); cols <- cols_func(400)
# plot_ly(data1, x = ~x, y = ~y, z = ~z, type = 'scatter3d', mode = 'lines', showlegend = F,
# line = list(width = 40, color = toRGB(cols, alpha = wave_opacity)))
p <-
plot_ly(data1, x = ~x, y = ~y, z = ~z, type = 'scatter3d', mode = 'lines', showlegend = F,
line = list(width = 40, color = ~c, colorscale = list(c(0,'red'), c(1,'blue')),
alpha = wave_opacity)) %>%
# trace the z-axis
add_trace(data1, x = 0, y =0, z = ~z, type="scatter3d", mode="lines",
line = list(width = 10, color = ~c, colorscale = list(c(0,'red'), c(1,'blue'))),
name="Z", hoverinfo="none") %>%
# add projection on plane x=-3
add_trace(data1, x = -3, y = ~y, z = ~z, type="scatter3d", mode="lines",
line = list(width = 2, dash="solid", color = ~c,
colorscale = list(c(0,'gray'), c(1,'black'))),
name="Z", hoverinfo="none") %>%
# add projection on plane y=-3
add_trace(data1, x = ~x, y = -3, z = ~z, type="scatter3d", mode="lines",
line = list(width = 2, dash="solid", color = ~c,
colorscale = list(c(0,'gray'), c(1,'black'))),
name="Z", hoverinfo="none") %>%
# add a z-line at x=-3,y=0
add_trace(data1, x = -3, y = 0, z = ~z, type="scatter3d", mode="lines",
line = list(width = 5, dash="solid", color = "gray"),
name="Z", hoverinfo="none") %>%
# add a z-line at x=0,y=-3
add_trace(data1, x = 0, y = -3, z = ~z, type="scatter3d", mode="lines",
line = list(width = 5, dash="solid", color = "gray"),
name="Z", hoverinfo="none") %>%
# add a x-line at y=0, z=-40
add_trace(data1, x = ~x, y = 0, z = -40, type="scatter3d", mode="lines",
line = list(width = 15, dash="solid", color = "gray"),
name="Z", hoverinfo="none") %>%
layout(scene = scene)
# add a few annotations/arrows
p2 Figure 3.3: Fourier Transform of images
Let’s demonstrate the use of R’s fft() to calculate the 2D FT of images. Information on data acquisition frequency and block length (in sec) cannot be included into the fft() call.
R generates a single 2D vector of the same dimensions as the data containing a list of complex numbers. The function Mod(fft()) is used to extract the magnitudes of the Fourier coefficients, which are computed by: \[magnitude = \sqrt{(real \times real + imaginary \times imaginary)},\] where the real and imaginary components are extracted by functions Re(fft()) and Im(fft()), respectively.
The method fft() generates only meaningful frequency up to half the sampling frequency. The FT returns values of the discrete Fourier transform for both positive and negative frequencies. Although, as sampling a signal in discrete time intervals causes aliasing problems, R yields all frequencies up to the sampling frequency. For instance, sampling a 50 Hz sine wave and 950 Hz sine wave with 1000 Hz will generate identical results, as the FT cannot distinguish between the two frequencies. Hence, the sampling frequency must always be at least twice as high as the expected signal frequency. For each actual frequency in the signal, the FT will give 2 peaks (one at the “actual” frequency and one at sampling frequency minus “actual” frequency). This will make the second half of the magnitude vector a mirror image of the first half.
As long as the sampling frequency is at least twice as high as the expected signal frequency, all meaningful information will be contained in the first half of the magnitude vector. However, a peak in the low frequency range might be present when high “noise” frequency is present in the signal (or image). At this point, the vector of extracted magnitudes is only indexed by the frequencies but has no associated frequencies. To calculate the corresponding frequencies, the FT simply takes (or generates) the index vector (1, 2, 3, …, length(magnitude vector)) and divides it by the length of the data block (in sec).
In 1D, the phases would represent a vector of the same length as the magnitude vector with the phases (0 to \(2\pi\) or \(-\pi\) to \(+\pi\)) of each frequency. Phase shifts are translations in space (e.g., x-axis) for a given wave component that are measured in angles (radians). For instance, shifting a wave \(f(x)=0.5\sin(3wt)+0.25\sin(10wt)\) by \(\frac{\pi}{2}\) would produce the following Fourier series:
\[f(t)=0.5\sin\left (3wt+\frac{\pi}{2}\right )+0.25\sin\left (10wt+\frac{\pi}{2}\right ).\]
In 2D, The Fourier transform (FT/IFT) for images is defined by: \[\hat{f}(u,v)=F(u,v)=\int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}{f(x,y)e^{-i2\pi(ux+vy)}dxdy}},\] \[f(x,y)=\hat{\hat{f}}(x,y)=\hat{F}(x,y)=\int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}{F(u,v)e^{i2\pi(ux+vy)}dudv}},\] where \(u\) and \(v\) are the spatial frequencies, \(F(u,v)=F_R(u,v)+iF_I(u,v)\) is a complex number for each pair of arguments, \[|F(u,v)|=\sqrt{F_R^2(u,v)+F_I^2(u,v)}\] is the magnitude spectrum, and \[\arctan\left (\frac{F_I(u,v)}{F_R(u,v)}\right )\] is the phase angle spectrum.
The complex exponential \[e^{-i2\pi(ux+vy)}=\cos(2\pi(ux+vy)) +i\ \sin(2\pi(ux+vy))\] represents the real and imaginary (complex) sinusoidal terms in the 2D plane. The extrema of its real part (\(\cos(2\pi(ux+vy))\)) occurs at \(2\pi(ux+vy)=n\pi\). Using vector notation, \[2\pi(ux+vy)=2\pi \langle U, X \rangle =n\pi,\] where the extrema points \(U=(u,v)^T\) and \(X=(x,y)^T\) represent sets of equally spaced parallel lines with normal \(U\) and wavelength \(\frac{1}{\sqrt{u^2+v^2}}\).
Let’s define the index shifting paradigm associated with the discrete FT, which is simply used for convenience and better visualization. It has no other relevance to the actual calculation of the FT and its inverse, IFT.
When applying the forward or reverse generalized discrete FT it is possible to shift the transform sampling in time and frequency domain by some real offset values, \(a,b\). Symbolically,
\[\hat{f}(k) = \sum_{n=0}^{N-1} f(n) e^{-\frac{2 \pi i}{N} (k+b) (n+a)} \quad \quad k = 0, \dots, N-1.\]
Note: Remember that in R, the array indices start with 1, not 0, as in some other languages.
The function fftshift() is useful for visualizing the Fourier transform with the zero-frequency component in the middle of the spectrum. Its inverse counterpart, ifftshift(), is needed to rearrange again the indices appropriately after the IFT is employed, so that the image is correctly reconstructed in spacetime. The FT only computes half of the frequency spectrum corresponding to the non-negative (positive and zero if the length(f) is odd) frequencies in order to save computation time. To preserve the dimensions of the output \(\hat{f}=FT(f)\), the second half of the frequency spectrum (the complex conjugate of the first half) is just added at the end of this vector. In a 1D setting, the result of fft() is:
\(0\ 1\ 2\ 3\ ...\ (freq\ bins > 0)\ ... {\frac{Fs}{2}}\) and \(-{\frac{Fs}{2}}\ ... \ (freq\ bins < 0)\ ...\ -3\ -2\ -1.\)
where \(F_s\) is the frequency sample. The fftshift method sets the zero-frequency component in the center of the array, i.e., it just shifts (offsets) the second part with the negative frequency bins to the beginning and the first part to the end of the resulting FT vector, or matrix. Thus, the shifted discrete FT can be nicely plotted in the center covering the frequency spectrum from \(-{\frac{Fs}{2}}\) on the left to \({\frac{Fs}{2}}\) on the right. This is not necessary, but is used for better visualization aesthetics. To synthesize back the correct image, after using fftshift on the FT signal, we always have to undo that re-indexing by using ifftshift() on the inverse-FT.
# FFT SHIFT
fftshift <- function(img_ff, dim = -1) {
rows <- dim(img_ff)[1]
cols <- dim(img_ff)[2]
# planes <- dim(img_ff)[3]
swap_up_down <- function(img_ff) {
rows_half <- ceiling(rows/2)
return(rbind(img_ff[((rows_half+1):rows), (1:cols)], img_ff[(1:rows_half), (1:cols)]))
}
swap_left_right <- function(img_ff) {
cols_half <- ceiling(cols/2)
return(cbind(img_ff[1:rows, ((cols_half+1):cols)], img_ff[1:rows, 1:cols_half]))
}
#swap_side2side <- function(img_ff) {
# planes_half <- ceiling(planes/2)
# return(cbind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)], img_ff[1:rows, 1:cols, 1:planes_half]))
#}
if (dim == -1) {
img_ff <- swap_up_down(img_ff)
return(swap_left_right(img_ff))
}
else if (dim == 1) {
return(swap_up_down(img_ff))
}
else if (dim == 2) {
return(swap_left_right(img_ff))
}
else if (dim == 3) {
# Use the `abind` package to bind along any dimension a pair of multi-dimensional arrays
# install.packages("abind")
library(abind)
planes <- dim(img_ff)[3]
rows_half <- ceiling(rows/2)
cols_half <- ceiling(cols/2)
planes_half <- ceiling(planes/2)
img_ff <- abind(img_ff[((rows_half+1):rows), (1:cols), (1:planes)],
img_ff[(1:rows_half), (1:cols), (1:planes)], along=1)
img_ff <- abind(img_ff[1:rows, ((cols_half+1):cols), (1:planes)],
img_ff[1:rows, 1:cols_half, (1:planes)], along=2)
img_ff <- abind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)],
img_ff[1:rows, 1:cols, 1:planes_half], along=3)
return(img_ff)
}
else {
stop("Invalid dimension parameter")
}
}
ifftshift <- function(img_ff, dim = -1) {
rows <- dim(img_ff)[1]
cols <- dim(img_ff)[2]
swap_up_down <- function(img_ff) {
rows_half <- floor(rows/2)
return(rbind(img_ff[((rows_half+1):rows), (1:cols)], img_ff[(1:rows_half), (1:cols)]))
}
swap_left_right <- function(img_ff) {
cols_half <- floor(cols/2)
return(cbind(img_ff[1:rows, ((cols_half+1):cols)], img_ff[1:rows, 1:cols_half]))
}
if (dim == -1) {
img_ff <- swap_left_right(img_ff)
return(swap_up_down(img_ff))
}
else if (dim == 1) {
return(swap_up_down(img_ff))
}
else if (dim == 2) {
return(swap_left_right(img_ff))
}
else if (dim == 3) {
# Use the `abind` package to bind along any dimension a pair of multi-dimensional arrays
# install.packages("abind")
library(abind)
planes <- dim(img_ff)[3]
rows_half <- floor(rows/2)
cols_half <- floor(cols/2)
planes_half <- floor(planes/2)
img_ff <- abind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)],
img_ff[1:rows, 1:cols, 1:planes_half], along=3)
img_ff <- abind(img_ff[1:rows, ((cols_half+1):cols), (1:planes)],
img_ff[1:rows, 1:cols_half, (1:planes)], along=2)
img_ff <- abind(img_ff[((rows_half+1):rows), (1:cols), (1:planes)],
img_ff[(1:rows_half), (1:cols), (1:planes)], along=1)
return(img_ff)
}
else {
stop("Invalid dimension parameter")
}
}Now we will present the FT and IFT of a pair of synthetic 2D images (a square and a disk).
library(imager)
# install EBImage
# source("https://bioconductor.org/biocLite.R")
# biocLite("EBImage")
library(EBImage)# Define two synthetic images
square_arr <- matrix(nrow=256, ncol=256)
circle_arr <- matrix(nrow=256, ncol=256)
for (i in 1:256) {
for (j in 1:256) {
if ( abs(i-128) < 30 && abs(j-128) < 30)
square_arr[i,j]=1 # sqrt((i-128)^2+(j-128)^2)/30
else square_arr[i,j]=0
if ( sqrt((i-128)^2 + (j-128)^2)<30)
circle_arr[i,j]=1 # 1-sqrt((i-128)^2+(j-128)^2)/30
else circle_arr[i,j]=0
}
}
# image(square_arr); image(circle_arr)
# display(square_arr, method = "raster"); display(circle_arr, method = "raster")
# download.file("https://wiki.socr.umich.edu/images/e/ea/BrainCortex2.png",
# paste(getwd(),"./image.png", sep="/"), mode = 'wb')
# im <- load.image(paste(getwd(),"./image.png", sep="/"))
# plot(im)
# Grayscaled
# img_gray <- im[ , , 1]
# display(img_gray, title='Image')
# FFT
# img_ff <- fft(img_gray) #fftw2d
ft_square <- fft(square_arr) #fftw2d # Display Re(FT): display(fftshift(ft_square))
ft_circle <- fft(circle_arr) # display(fftshift(ft_circle))
# Magnitude and Phase
mag_ft_square <- sqrt(Re(ft_square)^2+Im(ft_square)^2)
mag_ft_circle <- sqrt(Re(ft_circle)^2+Im(ft_circle)^2)
# Phase <- atan(Im(img_ff)/Re(img_ff))
phase_ft_square <- atan2(Im(ft_square), Re(ft_square))
phase_ft_circle <- atan2(Im(ft_circle), Re(ft_circle))
# FFT SHIFT
shift_ft_square <- fftshift(mag_ft_square)
shift_ft_circle <- fftshift(mag_ft_circle)
# Display FT
display(log(shift_ft_square),title="FT Magnitude")
display(log(shift_ft_circle),title="FT Phase")
# Magnitude and Phase
mag_shift_ft_square <- sqrt(Re(shift_ft_square)^2+Im(shift_ft_square)^2)
# phase <- atan(Im(img_ff)/Re(img_ff))
phase_ft_square <- atan2(Im(ft_square), Re(ft_square))
display(fftshift(phase_ft_square))
phase_ft_circle <- atan2(Im(ft_circle), Re(ft_circle))
display(fftshift(phase_ft_circle))
# Implicitly invert the FT (IFT)
fftinv <- function( x ) { fft( x, inverse=TRUE ) / length( x ) }
display(Re(fftinv(fft(square_arr))),title="(IFT o FT) Magnitude")
display(Re(fftinv(fft(circle_arr))),title="(IFT o FT) Magnitude")
# FT of Square # Display the FT with a shift or not: display(fftshift(Re(X1))
X1 = fft(square_arr); display(fftshift(Re(X1)), method = "raster")
X1_mag <- sqrt(Re(X1)^2+Im(X1)^2); display(fftshift(X1_mag), method = "raster") # magnitude only
X1_phase <- atan2(Im(X1), Re(X1)); display(fftshift(X1_phase), method = "raster") # phase only
# Implicit Automated IFT
hat_X1 = Re(fft(X1, inverse = T)/length(square_arr))
display(hat_X1, method = "raster")
# Manually invert the FT (IFT) using the magnitudes and phases
Real1 = X1_mag * cos(X1_phase)
Imaginary1 = X1_mag * sin(X1_phase)
man_hat_X1 = Re(fft(Real1+1i*Imaginary1, inverse = T)/length(X1))
display(man_hat_X1, method = "raster")
# FT of Circle # No shift applied here (perhaps should be consistent or just show the difference?)
X2 = fft(circle_arr)
display(Re(X2), method = "raster")
X2_mag <- sqrt(Re(X2)^2+Im(X2)^2)
display(X2_mag, method = "raster") # magnitude only
X2_phase <- atan2(Im(X2), Re(X2))
display(X2_phase, method = "raster") # phase only
# Implicit Automated IFT
# hat_X2 = Re(fft(X2, inverse = T)/length(circle_arr)); display(hat_X2, method = "raster")
# Manually invert the FT (IFT) using the magnitudes and phases
Real2 = X2_mag * cos(X2_phase)
Imaginary2 = X2_mag * sin(X2_phase)
man_hat_X2 = Re(fft(Real2+1i*Imaginary2, inverse = T)/length(X1))
display(man_hat_X2, method = "raster")
# IFT Square-Magnitude and Circle-Phase
Real = X1_mag * cos(X2_phase)
Imaginary = X1_mag * sin(X2_phase)
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1mag_X2phase, method = "raster")
# IFT Circle-Magnitude and Square-Phase
Real = X2_mag * cos(X1_phase)
Imaginary = X2_mag * sin(X1_phase)
ift_X1phase_X2mag = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1phase_X2mag, method = "raster")
# IFT Circle-Magnitude and Nil-Phase
Real = X2_mag * cos(0)
Imaginary = X2_mag * sin(0)
ift_NilPhase_X2mag = Re(ifftshift(fft(Real+1i*Imaginary, inverse = T)/length(X2)))
display(ift_NilPhase_X2mag, method = "raster")
# IFT Square-Magnitude and Nil-Phase
Real = X1_mag * cos(0)
Imaginary = X1_mag * sin(0)
ift_NilPhase_X1mag = Re(ifftshift(fft(Real+1i*Imaginary, inverse = T)/length(X1)))
display(ift_NilPhase_X1mag, method = "raster")