Ejemplo sencillo para calcular la TDF 2D¶
Métodos potenciales de prospección, FCAG, 2024.
import numpy as np # para calcular, numerical python.
import matplotlib.pyplot as plt # para graficar, matplotlib.
plt.xkcd()
plt.rcParams.update({'font.size': 14}) # tamaño de fuente
Repaso de señales en 1D ...¶
# Correr animación
from IPython.display import Image
Image("tdf.gif",width=400)
Ejemplo 1¶
El dato es una función constante, $f(x,y)=1$ y tomamos 16 muestras igualmente espaciadas en cada dirección.
# Subrutina para graficar:
def Plot(data,title="$f_{ij}$",cmap="coolwarm",edgecolors="gray"):
nx,ny = np.shape(data)
plt.figure(1,figsize=(5,5))
plt.pcolor(data, edgecolors=edgecolors, linewidths=2,cmap=cmap, shading="auto")
plt.title(title)
plt.xlabel("i")
plt.ylabel("j")
plt.xticks(range(nx))
plt.yticks(range(ny))
plt.show()
return True
cmap1 = "coolwarm"
cmap2 = "gray"
edgecolors = "gray"
nx,ny = 4,4
f = np.ones(nx*ny).reshape(nx,ny) # constante
Graficamos el dato 2D muestreado $f_{ij}$:
Plot(f);
Calculamos la transformada discreta bidimensional de Fourier (2D TDF) por medio de la TDF 1D. Calculamos la TDF de las filas y luego la TDF de las columnas. Utilizamos para ello numpy.fft.fft():
FF = np.zeros_like(f,dtype=complex)
# TDF por fila y por columna
for i in np.arange(nx):
FF[i,:] = np.fft.fft(f[i,:]) # filas
for j in np.arange(ny): # columnas (del arreglo anterior)
FF[:,j] = np.fft.fft(FF[:,j])
Graficamos el espectro de amplitud:
Plot(np.abs(FF),title="$|\mathcal{F}_{ij}|$",cmap="bone");
Lo mismo, pero en número de onda angular adimensional $0\leqslant \kappa_x,\kappa_y<2\pi$. Notamos que al no haber muestras sobre $2\pi$ se grafican bordes sin valores.
# En términos de los números de onda angulares adimensionales
kx = (2*np.pi)*np.arange(nx)/nx
ky = (2*np.pi)*np.arange(ny)/ny
edgecolors = "gray"
plt.figure(figsize=(5,5))
plt.pcolor(kx,ky,np.abs(FF), edgecolors=edgecolors, linewidths=2,cmap=cmap2, shading="auto")
plt.xticks([0,np.pi/2,np.pi,3*np.pi/2,2*np.pi],[r'$0$',r'$\pi$/2',r'$\pi$',r'$3\pi/2$',r'$2\pi$'])
plt.yticks([0,np.pi/2,np.pi,3*np.pi/2,2*np.pi],[r'$0$',r'$\pi$/2',r'$\pi$',r'$3\pi/2$',r'$2\pi$'])
plt.title("$|\mathcal{F}|_{}$")
plt.xlabel("$0 \leqslant \kappa_x < 2\pi$")
plt.ylabel("$0 \leqslant \kappa_y < 2\pi$")
plt.show()
Calculamos la TDF 2D con la subrutina numpy.fft.fft2() y comparamos:
F = np.fft.fft2(f)
print ("Check: ", np.allclose(F,FF))
Ejemplo 2¶
El dato es ahora una función cajón o boxcar 2D.
f = np.pad(f,(5,5),'constant')
F = np.fft.fft2(f)
plt.figure(3,figsize=(16,8))
plt.subplot(1,2,1)
plt.pcolor(f, edgecolors=edgecolors, linewidths=1,cmap=cmap1, shading="auto")
plt.title("$f_{ij}$")
plt.xlabel("i")
plt.ylabel("j")
plt.subplot(1,2,2)
plt.pcolor(np.abs(F), edgecolors=edgecolors, linewidths=1, cmap=cmap2, shading="auto")
plt.title("$|\mathcal{F}|_{ij}$")
plt.xlabel("i")
plt.ylabel("j")
plt.show()
Podemos hacer zero-padding para observar la TDF 2D en más puntos, con el objeto de reducir la incertidumbre al tomar muestras de la respuesta en frecuencia 2D del dato:
F = np.fft.fft2(f,s=(2**8,2**8)) # FFT 2D con padding dado por s=(,)
plt.figure(figsize=(15,10))
plt.subplot(1,2,2)
plt.imshow(np.abs(F), cmap=cmap2, origin="lower")
plt.title("$|\mathcal{F}|_{ij}$ con zero-padding")
plt.xlabel("i")
plt.ylabel("j")
plt.show()
Para terminar, si queremos graficar el resultado con la muestra en $(0,0)$ en el centro podemos utilizar np.fft.fftshift():
nx0,ny0 = np.shape(F)
plt.figure(4,figsize=(15,10))
plt.subplot(1,2,2)
plt.imshow(np.fft.fftshift(np.abs(F)), cmap=cmap2, origin="lower")
plt.title("$|\mathcal{F}|_{ij}$ con zero-padding")
plt.xlabel("i")
plt.ylabel("j")
plt.xticks([0,nx0//2,nx0],[-nx0//2,0,+nx0//2])
plt.yticks([0,ny0//2,ny0],[-ny0//2,0,+ny0//2])
plt.show()
Eso es todo por hoy.