Raghav K. Chhetri¶
- Discrete Fourier Transform (DFT) is a way to transform a discrete signal from time to frequency domain by summing over a sequence of sine and cosine waves. Fast Fourier Transform (FFT) is a series of algorithms that compute DFT by factorizing the matrix into a product of sparse matrices, thereby reducing the complexity and speed of computation enormously
There are multiple ways to compute FFT in Python:
numpy.fftscipy.fftpack-- Considered legacy. SciPy recommends usingscipy.fftscipy.fftpyfftwcupy.fft-- NVIDIA's library that evaluates FFTs on GPUstorch.FFT
- 2D FFT speed test shows
torch.fftas the fastest option scipy.fftis slightly faster thannumpy.fft, so using that in this notebook
import sys
print('Environment: {}'.format(sys.exec_prefix))
print('Executable: {}'.format(sys.executable))
print('Version: {}'.format(sys.version))
import numpy as np
import scipy.misc
from scipy import ndimage, signal
from scipy.fft import fft, ifft, fftn, ifftn, fftshift, ifftshift, fftfreq, rfft, irfft, rfftfreq
import skimage
from skimage.util import img_as_float
from skimage.data import camera, gravel, checkerboard
from skimage.filters import difference_of_gaussians, window
import matplotlib
from matplotlib import cm
import matplotlib.pyplot as plt
#plt.style.use('seaborn-poster')
%matplotlib inline
print('numpy: {}'.format(np.__version__))
print('scipy: {}'.format(scipy.__version__))
print('skimage: {}'.format(skimage.__version__))
print('matplotlib: {}'.format(matplotlib.__version__))
1D FFT (for time-series data)¶
General comments:
fft takes longer if npnts is not a multiple of 2 and even longer if it is a prime
The frequency axis has been defined in different ways in fft_powerspectrum_plot below:
Method 1 (Frequency mirroring above Nyquist):\ Manually define the frequency axis in this range,
[0, Sampling Rate]. The FFT Amplitude of the signal now has frequency peaks that are mirrored in the[Nyquist, Sampling Rate]too. Ignore the mirrored peaks and only consider this range,[0 to Nyquist frequency]. Recall that Nyquist frequency is half of the Sampling Rate.Method 2 (Frequency mirroring about zero):\ Define the frequency axis using
fftfreqwhich centers the output offft(). The frequency axis now extends from[-Nyquist, Nyquist]and the peaks are mirrored in the negative axis. Again, ignore the negative frequency and only consider this range,[0 to Nyquist frequency]Method 3 (No mirroring-- only positive frequencies):\ The frequency spectrum that
fft()outputted was reflected about the y-axis, which is caused by inputting real numbers instead of complex numbers tofft(). This symmetry can be taken advantage of to make Fourier transform faster by computing only the positive frequencies, which is whatrfft()does.
def fft_powerspectrum_plot(x, t, dt, npnts, freq_show, method=2, denoise=False, cutoff=0):
'''
- plot original signal
- compute FFT and power spectrum (power per frequency)
- plot power spectrum
- plot recovered signal after fft-> ifft
x: 1D signal
t: time vector
dt: sampling interval
npnts: number of time points
freq_show: upper range of frequency to plot
method: 1/2/3
denoise: filter out noise
cutoff: denoise below the cutoff amplitude
Call: plot_signal_fftamplitude(x, t, dt, npnts, 200, method=1, denoise=True, cutoff=3)
'''
if method == 1:
X = fft(x)
n = np.arange(npnts)
T = npnts*dt
freq = n/T
title = 'Frequency mirroring above Nyquist'
elif method == 2:
X = fft(x)
freq = fftfreq(npnts, dt)
title = 'Frequency mirroring about zero'
elif method == 3:
X = rfft(x)
freq = rfftfreq(npnts, dt)
title = 'No mirroring: only positive frequencies'
powerspect = 2*np.abs(X)/npnts
# Note: Normalized as 2*np.abs(X)/npnts instead of simply np.abs(X)
# returns the actual amplitude values of the sine and cosine functions
# instead of some arbitrarily-scaled values
if denoise:
powerspect = powerspect * (powerspect > cutoff) # Zero all frequencies with small power
X = X * (powerspect > cutoff) # Zero small Fourier coefficients
#To further zero out a peak at zero frequency -- occurs if noise is 'rand' instead of 'randn'
#powerspect = powerspect * (powerspect < 10)
#X = X * (powerspect < 10)
# PLOT
plt.figure(figsize = (9, 7))
plt.subplot(311)
plt.plot(t, x, 'k',label='original')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.subplot(312)
plt.title(title)
#plt.stem(freq, np.abs(X),'c', markerfmt=" ", basefmt="-b")
plt.stem(freq, powerspect,'c', markerfmt=" ", basefmt="-b")
plt.xlabel('Freq (Hz)')
plt.ylabel('FFT Amplitude')
if method == 2:
plt.xlim(-freq_show,freq_show)
else:
plt.xlim(0,freq_show)
plt.subplot(313)
if method == 3:
plt.plot(t, irfft(X), 'k--',label='recovered')
else:
plt.plot(t, ifft(X), 'k--',label='recovered')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.tight_layout()
plt.show()
print('Number of time points:', npnts, 'points')
print('Number of points in frequency range:', len(freq), 'points')
print('Frequency range:', min(freq), max(freq), 'Hz')
Let's generate a 1D signal and add some noise to it:
srate = 1000 # Sampling rate in Hz
dt = 1/srate
duration = 1 # in seconds
t = np.arange(0,duration,dt) # Time Vector in seconds
npnts = len(t)
# 1D signal
freq1 = 10.0
x = 2*np.sin(2*np.pi*freq1*t)
freq2 = 70.0
x += 4*np.cos(2*np.pi*freq2*t)
freq3 = 170.0
x += 6*np.cos(2*np.pi*freq3*t+np.pi/4)
noise = 6.0*np.random.randn(npnts)
x += noise
Let's plot to see what the signal looks like and what its power spectrum looks like (no denoising yet):
fft_powerspectrum_plot(x, t, dt, npnts, 200, 1)
Now, let's denoise frequencies with power spectrum below 1:
fft_powerspectrum_plot(x, t, dt, npnts, 200, 1, True, 1)
If we plot the frequency axis to show mirrored values in the negative axis too, this is what it looks like:
fft_powerspectrum_plot(x, t, dt, npnts, 200, 2, True, 1)
Or, using method 3, which uses rfft instead of fft:
fft_powerspectrum_plot(x, t, dt, npnts, 200, 3, True, 1)
Spurious frequency peaks in Amplitude Spectrum¶
Let's change the frequencies of the signal to non-integer values. When we do so, we see that the FFT amplitude shows up with spurious frequency peaks
srate = 1000 # Sampling rate in Hz
dt = 1/srate
duration = 1 # in seconds
t = np.arange(0,duration,dt) # Time Vector in seconds
npnts = len(t)
# 1D signal
freq1 = 10.1
x = 2*np.sin(2*np.pi*freq1*t)
freq2 = 70.5
x += 4*np.cos(2*np.pi*freq2*t)
freq3 = 170.7
x += 6* np.cos(2*np.pi*freq3*t+np.pi/4)
noise = 6.0*np.random.randn(npnts)
x += noise
fft_powerspectrum_plot(x, t, dt, npnts, 200)
One way to mitigate these additional frequency components is to use a longer signal, as shown below:
srate = 1000 # Sampling rate in Hz
dt = 1/srate
duration = 10 # in seconds
t = np.arange(0,duration,dt) # Time Vector in seconds
npnts = len(t)
# 1D signal
freq1 = 10.1
x = 2*np.sin(2*np.pi*freq1*t)
freq2 = 70.5
x += 4*np.cos(2*np.pi*freq2*t)
freq3 = 170.7
x += 6*np.cos(2*np.pi*freq3*t+np.pi/4)
noise = 6.0*np.random.randn(npnts)
x += noise
fft_powerspectrum_plot(x, t, dt, npnts, 200)
Another option is to denoise lower frequencies like we did earlier
2D FFT (working with images)¶
1. Blurring an image with 2D FFT¶
Adapted from the following numpy example
tmp = scipy.misc.face(gray=True) #Racoon face from SciPy
image = tmp[128:640, 256:768]
image.shape
plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.title('Original image')
plt.axis('off')
plt.imshow(image,cmap='gray');
# Take the 2-dimensional FFT and center the frequencies
ftimage = fftn(image)
ftimage = fftshift(ftimage)
plt.subplot(122)
#plt.imshow(np.abs(ftimage))
plt.imshow(np.log(np.abs(ftimage)), cmap='gray')
plt.title('FT(image)')
plt.axis('off')
plt.show()
# Build and apply a Gaussian filter
sigmax, sigmay = 7, 13
size = image.shape[0]
cy, cx = size/2, size/2
x = np.linspace(0, size, size)
y = np.linspace(0, size, size)
X, Y = np.meshgrid(x, y)
gmask = np.exp(-(((X-cx)/sigmax)**2 + ((Y-cy)/sigmay)**2))
#Note: mask is in Fourier space; same size as FT(image)
plt.figure(figsize = (12, 6))
plt.subplot(131)
plt.imshow(gmask, cmap='gray')
plt.title('Gaussian mask')
ftimagep = ftimage * gmask
plt.subplot(132)
plt.imshow(np.log(np.abs(ftimagep)),cmap='gray')
plt.title('FT(image)*mask')
plt.axis('off')
# Finally, take the inverse transform and show the blurred image
imagep = np.fft.ifftn(ftimagep)
plt.subplot(133)
plt.imshow(np.abs(imagep), cmap='gray')
plt.title('Blurred image')
plt.axis('off')
plt.show()
2. Gaussian blur implemented using FFT convolution: fftconvolve¶
#From https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.fftconvolve.html
window_x = signal.windows.gaussian(101, std=7)
window_y = signal.windows.gaussian(101, std=13)
plt.plot(window_x, label= "$\sigma_x$=7")
plt.plot(window_y, label= "$\sigma_y$=13")
plt.legend(loc="upper right")
plt.title("Gaussian window")
plt.ylabel("Amplitude")
plt.xlabel("Sample")
kernel = np.outer(window_x, window_y)
#Note: Kernel is in the same coordinate as the image; Its size is smaller compared to image
plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.imshow(kernel, cmap='gray')
plt.title('Gaussian Kernel')
blurred = signal.fftconvolve(image, kernel, mode='same')
plt.subplot(122)
plt.imshow(blurred, cmap='gray')
plt.title('Blurred image')
plt.axis('off')
plt.show()
Notice the dark borders around the image, due to zero-padding beyond its boundaries.
scipy.ndimage.gaussian_filter gets rid of this artifact.
# Image Denoising: Gaussian filter, Median filter
noisy = image + 0.4*image.std()*np.random.random(image.shape)
gauss_denoised = ndimage.gaussian_filter(noisy, sigma = 11)
med_denoised = ndimage.median_filter(noisy, size = 9)
plt.figure(figsize=(12,4))
plt.subplot(131)
plt.imshow(noisy, cmap=plt.cm.gray, vmin=40, vmax=220)
plt.axis('off')
plt.title('noisy', fontsize=20)
plt.subplot(132)
plt.imshow(gauss_denoised, cmap=plt.cm.gray, vmin=40, vmax=220)
plt.axis('off')
plt.title('Gaussian filter', fontsize=20)
plt.subplot(133)
plt.imshow(med_denoised, cmap=plt.cm.gray, vmin=40, vmax=220)
plt.axis('off')
plt.title('Median filter', fontsize=20)
plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, bottom=0, left=0, right=1)
plt.show()
Note that a Gaussian filter smooths out the noise and the edges. Median filter averages the noise and preserves the edges better
Band-pass filtering by Difference of Gaussians¶
Adapted from the following SciKit-Image documentation
Band-pass filters attenuate signal frequencies outside of a range (band) of interest. In image analysis, they can be used to denoise images while at the same time reducing low-frequency artifacts such a uneven illumination. Band-pass filters can be used to find image features such as blobs and edges.
One method for applying band-pass filters to images is to subtract an image blurred with a Gaussian kernel from a less-blurred image. This example shows two applications of the Difference of Gaussians approach for band-pass filtering.
More on DoG here
1. Denoise image and reduce shadows¶
image = gravel()
wimage = image * window('hann', image.shape) # window image to improve FFT
filtered_image = difference_of_gaussians(image, 1, 12)
filtered_wimage = filtered_image * window('hann', image.shape)
im_f_mag = fftshift(np.abs(fftn(wimage)))
fim_f_mag = fftshift(np.abs(fftn(filtered_wimage)))
fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(11, 12))
ax[0, 0].imshow(image, cmap='gray')
ax[0, 0].set_title('Original Image')
ax[0, 0].axis('off')
ax[0, 1].imshow(np.log(im_f_mag), cmap='magma')
ax[0, 1].set_title('FFT Magnitude (log)')
ax[0, 1].axis('off')
ax[1, 0].imshow(filtered_image, cmap='gray')
ax[1, 0].set_title('Filtered Image')
ax[1, 0].axis('off')
ax[1, 1].imshow(np.log(fim_f_mag), cmap='magma')
ax[1, 1].set_title('FFT Magnitude (log)')
ax[1, 1].axis('off')
plt.show()
2. Enhance edges in an image¶
image = camera()
#image = noisy
wimage = image * window('hann', image.shape) # window image to improve FFT
filtered_image = difference_of_gaussians(image, 1.5)
filtered_wimage = filtered_image * window('hann', image.shape)
im_f_mag = fftshift(np.abs(fftn(wimage)))
fim_f_mag = fftshift(np.abs(fftn(filtered_wimage)))
fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(11, 12))
ax[0, 0].imshow(image, cmap='gray')
ax[0, 0].set_title('Original Image')
ax[0, 0].axis('off')
ax[0, 1].imshow(np.log(im_f_mag), cmap='magma')
ax[0, 1].set_title('Original FFT Magnitude (log)')
ax[0, 1].axis('off')
ax[1, 0].imshow(filtered_image, cmap='gray')
ax[1, 0].set_title('Filtered Image')
ax[1, 0].axis('off')
ax[1, 1].imshow(np.log(fim_f_mag), cmap='magma')
ax[1, 1].set_title('Filtered FFT Magnitude (log)')
ax[1, 1].axis('off')
plt.show()
Diffraction-limited imaging simulation¶
- Diffraction-limited applies to a system when the effect of aberrations to the pupil wavefront is smaller than the effect of diffraction at the pupil boundary
Relationship between imaging properties:
cPSF(u,v) ⟸ FT ⟹ CTF(Fu,Fv) | | Mag. squared Autocorrelation ↓ 🡓 iPSF(u,v) ⟸ FT ⟹ OTF(Fu,Fv)where
(u,v)= image plane co-ordinatescPSF= coherent point spread function or coherent impulse response function,h(u,v)iPSF= intensity point spread function or incoherent impulse response function,PSF(u,v)CTF= coherent transfer function or amplitude transfer function,H(Fu,Fv)OTF= optical transfer function,OTF(Fu,Fv)Note that modulation transfer function,
MTF = |OTF|
- Examples below adapted from Computational Fourier Optics: A MATLAB Tutorial, Chapter 7, page 123
Coherent Imaging Simulation:¶
from google.colab import drive
drive.mount('/content/drive')
cd /content/drive/MyDrive/Colab Notebooks/
#Coherent Image Simulation
img = plt.imread('./USAF_Resolution_Chart.jpeg') # RGB image, 8-bit
img = img[:,:,0] # grayscale image, 8-bit
M, N = img.shape
L = 0.9e-3 # image length in m
du = L/M # sample interval (m)
print(f'Object is {L:.0e} m ({M} pixels), sampled every {du:.2e} m')
Nyquistf = 1/(2*du)
print(f'Thus, the highest frequency avaiable in an ideal image (Nyquist freqeuncy) is {Nyquistf:.2e} cycles/m')
# Image space co-ordinates
u = np.arange(-L/2, L/2, du)
v = np.arange(-L/2, L/2, du)
# binary Image, type: float
Ig = np.where(img > 128, 1.0, 0.0) #img>128, replace by 1.0, otherwise replace by 0.0
ug = np.sqrt(Ig) # an ideal geometric-optics predicted image field
plt.figure(figsize=(7,7))
plt.title('Ideal image')
plt.imshow(ug, extent=[-L/2, L/2, -L/2, L/2], cmap='gray')
plt.ticklabel_format(style="sci", scilimits=(0,0))
plt.xlabel('u (m)')
plt.ylabel('v (m)')
plt.show()
print('Note that u and v are image plane co-ordinates')
lamda = 0.5e-6 # wavelength in m
wxp = 6.25e-3 # exit pupil radius in m
zxp = 125e-3 # exit pupil distance in m
fnumber = zxp/(2*wxp)
print(f'f# of the lens is {fnumber}')
f0 = wxp/(lamda*zxp) # coherent cutoff frequency
print(f'Coherent cutoff frequency with this lens is {f0:.2e} cycles/m')
if Nyquistf >= 2*f0:
print('Imaging is in the diffraction-limited regime')
else:
print('Not in the diffraction-limited regime: object is undersampled')
#frequency coordinates
fu = np.arange(-1/(2*du), 1/(2*du), 1/L)
fv = np.arange(-1/(2*du), 1/(2*du), 1/L)
Fu,Fv = np.meshgrid(fu,fv)
H = np.sqrt(Fu**2 + Fv**2)/f0 < 1 # Coherent Transfer Function (Amplitude Transfer Function)
fig = plt.figure(figsize=(9,9))
fig.suptitle('Coherent Transfer Function')
ax = fig.add_subplot(projection='3d')
# Plot the surface
ax.plot_surface(Fu, Fv, H, cmap=cm.gray)
plt.ticklabel_format(style="sci", scilimits=(0,0))
fig.tight_layout()
ax.set_xlabel('fu (cycle/m)')
ax.set_ylabel('fv (cycle/m)')
plt.show()
Note that the ideal image field is blurred by the impulse response function of the imaging system, h(u,v)= FFT[H(Fu,Fv)]
i.e., Image(u,v) = Ideal_Image(u,v) $*$ h(u,v)
This is implemented in frequency-domain as follows:
Gi(Fu,Fv)= H(Fu,Fv).FFT(Ideal_Image)
So, actual image field, ui(u,v) = iFFT(Gi)
H = fftshift(H) #Coherent transfer function
Gg = fftn(fftshift(ug)) #FFT of an ideal geometric-optics predicted image field
Gi = np.multiply(Gg, H)
ui = ifftshift(ifftn(Gi))
Ii = (np.abs(ui))**2
#Simulated diffraction-limited coherent image
plt.figure(figsize=(7,7))
plt.title('Simulated image')
plt.imshow(Ii**0.5, extent=[-L/2, L/2, -L/2, L/2], cmap='gray')
plt.ticklabel_format(style="sci", scilimits=(0,0))
plt.xlabel('u (m)')
plt.ylabel('v (m)')
plt.show()
Next, let's account for surface roughness of the object by adding a random complex exponential phase term to the field:
ug1 = np.sqrt(Ig)*np.exp(1j*2*np.pi*np.random.randn(M))
H1 = H
Gg1 = fftn(fftshift(ug1))
Gi1 = np.multiply(Gg1, H1)
ui1 = ifftshift(ifftn(Gi1))
Ii1 = (np.abs(ui1))**2
#Simulated diffraction-limited coherent image
plt.figure(figsize=(7,7))
plt.title('Simulated image of an object w/ a rough surface')
plt.imshow(Ii1**0.5, extent=[-L/2, L/2, -L/2, L/2], cmap='gray')
plt.ticklabel_format(style="sci", scilimits=(0,0))
plt.xlabel('u (m)')
plt.ylabel('v (m)')
plt.show()
The irradiance image that we get is thus impeded further by the prescence of speckles (more apparent with imshow Ii**0.5 so plotting that above)
Incoherent Image simulation:¶
PSF = np.abs(fftn(fftshift(H)))**2
OTF = ifftn(PSF) #left in the shifted arrangement
OTF = OTF/OTF[0,0] #normalize wrt zero-frequency value
MTF = ifftshift(np.abs(OTF))
fig = plt.figure(figsize=(9,9))
fig.suptitle('Incoherent Transfer Function')
ax = fig.add_subplot(projection='3d')
# Plot the surface
ax.plot_surface(Fu, Fv, MTF, cmap=cm.gray)
plt.ticklabel_format(style="sci", scilimits=(0,0))
fig.tight_layout()
ax.set_xlabel('fu (cycle/m)')
ax.set_ylabel('fv (cycle/m)')
plt.show()
Gg2 = fftn(fftshift(Ig)) #FFT of an ideal geometric-optics predicted image intensity
Gi2 = np.multiply(Gg2, OTF)
Ii2 = ifftshift(ifftn(Gi2))
Ii2 = np.real(Ii2) #remove residual imaginary parts
Ii2 = np.where(Ii2 <= 0, 0, Ii2) #wherever Ii2 is <= 0, replace by 0, otherwise leave as is
#Simulated diffraction-limited coherent image
plt.figure(figsize=(7,7))
plt.title('Simulated image')
plt.imshow(Ii2, extent=[-L/2, L/2, -L/2, L/2], cmap='gray')
plt.ticklabel_format(style="sci", scilimits=(0,0))
plt.xlabel('u (m)')
plt.ylabel('v (m)')
plt.show()
Notice that the image quality in case of incoherent imaging is better than than with coherent imaging. Also, there are no ringing artifacts or speckle features in the image.
!pwd