Fourier Optics

Santiago R, Birge Ken Tok

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import cv2
import pandas as pd
import math
from IPython.display import clear_output

Amplitude Modulation

In [2]:
def I_A(A, a,b,c,d):
    I = a*np.cos(b*A+c)**2+d
    return I
popt = (2,2,1,1)
A = np.arange(0,10,0.1)
plt.plot(A, I_A(A, *popt))
Out[2]:
[<matplotlib.lines.Line2D at 0x7fef60750590>]

Fourier Transforms

In [3]:
img = cv2.imread('NasaLogo.png',0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20*np.log(np.abs(fshift))

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()
In [4]:
def fourier_image(input_path, Name):
    img = cv2.imread(input_path,0)
    f = np.fft.fft2(img)
    fshift = np.fft.fftshift(f)
    magnitude_spectrum = 20*np.log(np.abs(fshift))
    plt.imshow(magnitude_spectrum, cmap = 'gist_gray')
    plt.xticks([]), plt.yticks([])
    plt.savefig(Name, dpi=1200)
    plt.show()
In [5]:
fourier_image("NasaLogo.png", "NasaFft")
In [6]:
focal_lengths = [100,200] #mm

Pixel Pitch

$D = (278 \pm 2)mm$ Abstand Spiegel-Schirm
$\lambda = 532nm$
$M = (31 \pm 2)mm$ Abstand Spiegel-SLM

In [7]:
#in cm
d_low = [-4.2,-3.1,-2.1,-1.1]
d_high = [1.1,2.1,3.2,4.3]
d_right = [-4.2,-3.1,-2.1,-1]
d_left = [1, 2.1,3.1,4.2]
u_d = 0.2
u_m = 10
k = np.array([-4,-3,-2,-1,1,2,3,4])
In [8]:
def lin_fit(k,a,b):
    y = a*k+b
    return y
alpha_highlow = np.sin(np.arctan(np.array(d_low+d_high)*10/309))
u_alpha = np.sqrt((1/(1+(np.array(d_low+d_high)*10/309)**2)*1/309*u_d)**2+(1/(1+(np.array(d_low+d_high)*10/309)**2)*np.array(d_low+d_high)/(309**2)*u_m)**2)
u_sinalpha = np.cos(np.arctan(np.array(d_low+d_high)))*u_alpha
popt, pcov = curve_fit(lin_fit, k, np.sin(alpha_highlow), sigma = u_sinalpha)
plt.plot(k, lin_fit(k,*popt))
plt.scatter(k,alpha_highlow)
print("Der Pixelpitch beträgt", 532*1e-3/popt[0], "+/-", 532*1e-3/(popt[0]**2)*pcov[0,0]**0.5 ,"um")
print(popt)
print(pcov[0,0]**0.5)
print(pcov[1,1]**0.5)
Der Pixelpitch beträgt 15.687106781848366 +/- 0.049444311752608015 um
[0.0339132  0.00119024]
0.0001068912763243731
0.00035106050789269013
In [9]:
alpha_rightleft = np.sin(np.arctan(np.array(d_right+d_left)*10/309))
u_alpha = np.sqrt((1/(1+(np.array(d_right+d_left)*10/309)**2)*1/309*u_d)**2+(1/(1+(np.array(d_right+d_left)*10/309)**2)*np.array(d_right+d_left)/(309**2)*u_m)**2)
u_sinalpha = np.cos(np.arctan(np.array(d_right+d_left)))*u_alpha
popt, pcov = curve_fit(lin_fit, k, np.sin(alpha_rightleft), sigma = u_sinalpha)
plt.plot(k, lin_fit(k,*popt))
plt.scatter(k,alpha_rightleft)
print("Der Pixelpitch beträgt", 532*1e-3/popt[0], "+/-", 532*1e-3/(popt[0]**2)*pcov[0,0]**0.5 ,"um")
print(popt)
print(pcov[0,0]**0.5)
print(pcov[1,1]**0.5)
Der Pixelpitch beträgt 15.883818117196984 +/- 0.038896230709008585 um
[ 3.34932065e-02 -4.02404221e-09]
8.201803099780775e-05
0.00028893217718922
In [10]:
alpha_highlow = np.sin(np.arctan(np.array(d_low+d_high)*10/309))
u_alpha = np.sqrt((1/(1+(np.array(d_low+d_high)*10/309)**2)*1/309*u_d)**2+(1/(1+(np.array(d_low+d_high)*10/309)**2)*np.array(d_low+d_high)/(309**2)*u_m)**2)
u_sinalpha = np.cos(np.arctan(np.array(d_low+d_high)))*u_alpha
popt, pcov = curve_fit(lin_fit, k, np.sin(alpha_highlow), sigma = u_sinalpha)
plt.plot(k, lin_fit(k,*popt),label="Vertical Data Fit y = am+b")
plt.scatter(k,alpha_highlow,label="Data Vertical Diffraction")

alpha_rightleft = np.sin(np.arctan(np.array(d_right+d_left)*10/309))
u_alpha = np.sqrt((1/(1+(np.array(d_right+d_left)*10/309)**2)*1/309*u_d)**2+(1/(1+(np.array(d_right+d_left)*10/309)**2)*np.array(d_right+d_left)/(309**2)*u_m)**2)
u_sinalpha = np.cos(np.arctan(np.array(d_right+d_left)))*u_alpha
popt, pcov = curve_fit(lin_fit, k, np.sin(alpha_rightleft), sigma = u_sinalpha)
plt.plot(k, lin_fit(k,*popt), label="Horizontal DataFit y = am+b")
plt.scatter(k,alpha_rightleft,label="Data Horizontal Diffraction")
plt.ylabel("Angle "r'$sin(\theta _M)$')
plt.xlabel("Diffraction Order m")
plt.legend(loc="upper left")
plt.savefig("FitPixelPitch.pdf")

Amplitude Modulation

In [11]:
t = [0, 580,1440,1960,2880, 3460, 4100] #in schritten von pi in ms ab den 255-0 Sprung
u_t = 20 #in ms
#in sekunden
g_min = np.array([0.3])/4.14*256-1 
g_max = np.array([1])/4.14*256-1
p_t_min = np.array([1.4,1.46])
p_t_max = np.array([1.4,1.48])
print(g_min)
print(g_max)
p_t_min/4.14*256-1
[17.55072464]
[60.83574879]
Out[11]:
array([85.57004831, 89.28019324])
In [12]:
87*2
Out[12]:
174
In [13]:
data = pd.read_csv("Data_GreyValues.csv")
In [14]:
data
Out[14]:
Second Volt Unnamed: 2
0 -3.00 2.28
1 -2.99 2.36
2 -2.98 2.34
3 -2.97 2.48
4 -2.96 2.48
... ... ... ...
595 2.95 2.26
596 2.96 2.30
597 2.97 2.16
598 2.98 2.20
599 2.99 2.10

600 rows × 3 columns

In [15]:
t = np.array(data["Second"][60:440])
U = np.array(data["Volt"][60:440])
plt.plot(t,U)
Out[15]:
[<matplotlib.lines.Line2D at 0x7fef60214950>]
In [16]:
t = np.array(data["Second"])
U = np.array(data["Volt"])
plt.plot(t,U)
Out[16]:
[<matplotlib.lines.Line2D at 0x7fef602015d0>]
In [17]:
A = max(data["Volt"][60:440])-min(data["Volt"][60:440])
U_noise = min(data["Volt"][60:440])
Period = np.pi/1.5
Guesses = np.array([A,Period,0,U_noise])
popt, pcov = curve_fit(I_A,t[60:440],U[60:440], Guesses, sigma = 0.03*U[60:440])
t_array = np.arange(-3,3,0.02)
plt.plot(t_array, I_A(t_array, *popt),label="Fit")
plt.plot(t, U, label="Data")
plt.ylabel(r'$U in V \propto I$')
plt.xlabel("t in s")
plt.legend(loc="upper left")
print(popt)
print(U_noise)
[3.2343775  2.19075733 0.11077768 0.24562516]
0.18
In [18]:
#Fit mit A = U_max-U_min, d = U_min
def I_A(x,b,c):
    a = A
    d = U_noise
    I = a*np.cos(b*x+c)**2+d
    return I
Period = np.pi/1.5
Guesses = np.array([Period,0])
popt, pcov = curve_fit(I_A,t[60:440],U[60:440], Guesses, sigma = 0.03*U[60:440])
t_array = np.arange(-3,3,0.02)
plt.plot(t_array, I_A(t_array, *popt), label="Fit "r'$U(A) \propto I(A)$')
plt.plot(t, U, label="Data")
Gray_Values = np.array([1*np.pi/(4*popt[0])-2.68+0.12020481,3*np.pi/(4*popt[0])-2.68+0.12020481,2*np.pi/popt[0]-2.68+0.12020481])
plt.scatter(Gray_Values, I_A(Gray_Values, *popt), label="Gray Values", color = "black", s=50)
plt.ylabel(r'$U(t) \ in \ V \propto I(A)$')
plt.xlabel("t in s")
plt.legend(loc="upper right",fontsize=8)
plt.savefig("FitAmplitudeModulation.pdf")
print(popt)
print(pcov[0,0]**0.5)
print(pcov[1,1]**0.5)
print(A)
print(U_noise)
[2.19656736 0.12020481]
0.005516773591750394
0.007017437916589692
3.04
0.18
In [19]:
def residuals(f, popt, x, y, res_name):
    residuals = y - f(x,*popt)
    ss_res = np.sum(residuals**2)
    ss_tot = np.sum((y-np.mean(y))**2)
    R_2 = 1 - (ss_res / ss_tot)
    plt.scatter(x,residuals, label='Residuen, R^2 ='+str(np.round(R_2,7)))
    plt.xlabel("x-Werte")
    plt.ylabel("Differenz y-f(x)")
    plt.ylim(-np.abs(max(residuals))*1.5, np.abs(max(residuals))*1.5)
    plt.legend(loc="upper right", prop={'size': 8})
    plt.gca().set_aspect(aspect=int(1/(np.abs(max(residuals))*1.5)))
    #plt.rcParams["figure.figsize"] = (8,1)
    plt.savefig("Residuen"+res_name+".pdf", bbox_inches = "tight")
    print("R^2 =", R_2)
In [20]:
residuals(I_A, popt, t[60:440],U[60:440], "AmplitudenModulation")
R^2 = 0.9432242957177001
In [21]:
def A(phase):
    grey_value = phase/4.14*256-1
    return grey_value
A_2pi = A((2*np.pi)/popt[0])
A_max = A(3*np.pi/(4*popt[0]))
A_min = A(1*np.pi/(4*popt[0]))
print("A_2pi =",A_2pi)
print("A_max =",A_max)
print("A_min =",A_min)
A_2pi = 175.87846713845536
A_max = 65.32942517692076
A_min = 21.10980839230692

Gerchberg-Saxton Algorithm

In [22]:
from IPython.display import Image
Image(filename='Hologramme/3.jpg', width = 800, height = 300)
Out[22]:
In [23]:
#After this point, only images and their respective settings in the SLM setup are shown
#No further data processing is done, but all the recorded data can thus be saved in an
#orderly manner for potential future reference

Fourier Transforms

+ Comparison with FFT

In [24]:
path = "2021-04-19/Fourier/"
Image(filename=path+'abstract_black_white.jpg', width = 800, height = 300)
Out[24]:
In [25]:
path = "Hologramme/"
Image(filename=path+'abstract_black_white.jpg', width = 800, height = 300)
Out[25]: