Physics 141

# Start writing code here... #load the necessary Python libraries import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats import random import math %matplotlib inline from tabulate import tabulate

2a) The Nyquist-Shannon sampling theorem stipulates the sampling rate required to capture the full information in a continuous signal (in time) with a finite bandgap B. If we have a signal f(t) with maximal frequency component given by A<B, then we can determine the complete signal by taking measurements spaced by 1/(2B)[s]. (With A and B in Hz). In practice a lowpass filter is employed to control noise by high frequency components, as any components with A/2 or any f/2 element f(t) will create aliasing when sampled.

# 2a ![sampling image](sampling.webp)

n = 256 f1 = 5 f2 = 20 t = np.arange(0,np.pi,np.pi/n) amp = np.sin(f1*2*np.pi*t) + np.sin(f2*2*np.pi*t) plt.figure(figsize=(10,2)) plt.plot(t, amp) plt.show() FT = np.fft.fft(amp) d = t[1] - t[0] freq = 1/interval f = np.linspace(0,freq, n) plt.figure(figsize=(5,3)) plt.plot(f[:n//2], np.abs(FT)[:n//2] / (d*n)) #where d/n is normalization according to documentation # we plot to n/2 as FT includes n-delta(x-f1) and n-delta(x-f2) plt.show()

n = 257 f1 = 5 f2 = 20 t = np.arange(0,np.pi,np.pi/n) amp = np.sin(f1*2*np.pi*t) + np.sin(f2*2*np.pi*t) plt.figure(figsize=(10,2)) plt.plot(t, amp) plt.show() FT = np.fft.fft(amp) d = t[1] - t[0] freq = 1/interval f = np.linspace(0,freq, n) plt.figure(figsize=(5,3)) plt.plot(f[:n//2], np.abs(FT)[:n//2] /(d*n)) #where d/n is normalization according to documentation # we plot to n/2 as FT includes n-delta(x-f1) and n-delta(x-f2) plt.show()

The shape of the np.fft.fft raw output is a 1D array of length n, where n is the number of samples from our signal. The significance of the i-th bin of this output is the value corresponding to the frequency bin centered at the i-th value in np.fft.fftfreq = [0,1, ... , n/2 - 1, -n/2, ... , -1]/(d*n) if n is even, or np.fft.fftfreq = [0,1, ... , (n-1)/2, -(n-1)/2, ..., -1]/(d*n) if n is odd, where d is the sample spacing. The frequency resolution is thus an array of length n as well, which covers values centered at zero with a magnitude less than half of the inverse of the sample spacing. In our code, we plot only the positive half.

#2c from PIL import Image im = Image.open('/work/Core.jpg','r') im.convert('1') pixels = list(im.getdata())

image = np.asarray(pixels) finalImage = image.reshape(256,256,3) plt.imshow(finalImage) gray = lambda pic : np.dot(pic[... , :3] , [0.299 , 0.587, 0.114]) test = gray(finalImage) plt.imshow(test, cmap = plt.get_cmap('gray'))

#2c X = np.reshape(np.asarray(test),(256,256)) Y = np.zeros([256,256]) noise = np.random.uniform(0, 1, np.shape(X)) fT1 = np.fft.fftn(X) plt.imshow(np.log(np.abs(np.fft.fftshift(fT1))**2)) plt.show() fT2 = np.fft.fftn(noise) plt.imshow(np.log(np.abs(np.fft.fftshift(fT2))**2)) plt.show() fT3 = np.fft.fftn(X + noise) plt.imshow(np.log(np.abs(np.fft.fftshift(fT3))**2)) plt.show()

We see a centered cross at F(120,120). With coloration dropping around this center value toward the outside in a circular manner. We see uniform coloration for the noise as expected. The centering of this cross is done somewhat arbitrarily by DC shifting and corresponds to the average brightness (normally done at F(0,0) by convention). The F(n-1,n-1) corresponds to the highest frequency, as we increase frequency when distancing from the center. As low frequencies contain more information than high frequencies, we see a difference in magnitude as explained above. The (i,j)th pixel output of the 2D fourier transform is thus the magnitude of frequency contribution by a frequency greater f0 given by r = sqrt(fx^2+fy^2). Thinking about the maximum limit of alternating black and white pixels, f_max = 1/(2pixels) so f = f_max / 2, where now we have a 2D expression for f_max which is proportional to r (sqrt(kx^2+ky^2)). Therefore it should not be possible to include bins representing to frequencies above a 2D nyquist frequency.

#2d n = 256 k_freq = np.fft.fftfreq(n)*n # 2D mesh k_freq_2D = np.meshgrid(k_freq, k_freq) # normalize mesh k_nrm = np.sqrt(k_freq_2D[0]**2 + k_freq_2D[1]**2) # flatten as we dont need wave vectors k_nrm = k_nrm.flatten() #frequency bins k_bins = np.arange(-0.5, n//2+1, 1.) # k_max is half the pixel size k_val = 0.5*(k_bins[1:]+k_bins[:-1]) def powerspectrum(image): fTn = np.fft.fftn(image) f_amp = np.abs(fTn)**2 f_amp = f_amp.flatten() A_bins, _, _ = stats.binned_statistic(k_nrm, f_amp, statistic = "mean", bins = k_bins) A_bins *= np.pi * (k_bins[1:]**2 - k_bins[:-1]**2) plt.loglog(k_val,A_bins) plt.xlabel('k') plt.ylabel('P(k)') plt.show()

#2e powerspectrum(test) powerspectrum(noise) powerspectrum(test+noise)