Ejemplo sencillo de números de onda en 1D y 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': 16}) # tamaño de fuente
# graficos dentro de la notebook
#%matplotlib inline
Caso 1D¶
En 1D, el número de onda angular adimensional o digital es
$$\kappa_j=\frac{2\pi}{N}j.$$¶
Evaluamos para valores entre $[0,\pi) \cup [-\pi, 0)$ de manera de obtener correctamente $|\kappa_j|$ entre $[0,\pi]$ para $j=0,\cdots,N-1$. De esta manera podemos procesar adecuadamente una señal 1D de $N$ muestras para calcular filtros vinculados a derivadas.
nx = 4 # 3,4,5,10,100
kx = np.zeros(nx,dtype=float)
# a) Evaluando cada k y preguntando si k >= Nyquist == pi
dk = (2*np.pi/nx)
for j in np.arange(nx):
kx[j]=dk*j
if kx[j] >= np.pi:
kx[j] = kx[j]-2.0*np.pi
print(kx) # vemos a k
print(np.abs(kx)) # vemos a |k|
Visualización¶
Graficamos $\kappa$ vs. número de índice $n\,{\in}\,\mathbb{N}$.
index = np.arange(nx,dtype=int)
plt.figure(1,figsize=(10,5))
plt.ylabel('${\kappa}[n]$')
plt.xlabel('$n$')
plt.plot(index, kx, linestyle='', marker='o', markersize="10", color='tomato')
plt.xticks(index)
plt.yticks([-2*np.pi,-np.pi,0,np.pi,2*np.pi],
[r'$-2\pi$',r'$-\pi$',r'$0$',r'$+\pi$',r'$+2\pi$'])
plt.grid()
plt.show()
Comparamos con subrutina de numpy. La subrutina numpy.fftfreq calcula números de onda digitales entre $[0,+\frac{1}{2}){\cup}[-\frac{1}{2},0)$. Como recordarán de la Cátedra de Análisis de Señales en Geofísica, esto es equivalente (al multiplicar por $2\pi$) al número de onda angular digital entre $[0,+\pi){\cup}[-\pi,0)$, donde $\pi$ es la frecuencia de Nyquist.
kkx = (2*np.pi) * np.fft.fftfreq(nx)
print("check:", np.allclose(kx,kkx))
Caso 2D¶
Generamos el número de onda radial angular digital y tomamos su módulo. En 2D, el módulo del número de onda angular radial digital está dado por
$$|\kappa_r| = \sqrt{\kappa_x^2+\kappa_y^2}.$$Evaluamos para valores entre $[-\pi,\pi){\times}[-\pi,\pi)$ de manera de obtener correctamente $|\kappa_r|$ entre $[0,\pi]$. Para simplificar consideramos $\kappa_x=\kappa_y$.
ky,ny = kx,nx
kr = np.zeros(nx*ny).reshape(nx,ny)
# Utilizando dos bucles explícitos
for i in np.arange(nx):
for j in np.arange(ny):
kr[i,j] = np.sqrt(kx[i]**2+ky[j]**2)
# Utilizando un bucle explícito
# for i in np.arange(nx):
# kr[i,:] = np.sqrt(kx[i]**2+ky[:]**2)
Visualización¶
Graficamos $|\kappa_r|$ entre $[0,2\pi){\times}[0,2\pi)$.
cmap="inferno"
#cmap='gray'
plt.figure(2,figsize=(5,5))
plt.title(r"$|\kappa_r|$")
plt.xlabel(r"$\kappa_x$")
plt.ylabel(r"$\kappa_y$")
plt.xticks([0,np.pi,2*np.pi],[r'$0$',r'$\pi$',r'$2\pi$'])
plt.yticks([0,np.pi,2*np.pi],[r'$0$',r'$\pi$',r'$2\pi$'])
extent=[0,2*np.pi,0,2*np.pi]
plt.imshow(kr,origin="lower",cmap=cmap,extent=extent,interpolation='none')
plt.show()
Comparamos con subrutina de NumPy. En particular, numpy.meshgrid calcula matrices con las coordenadas del producto cartesiano $\kappa_x{\times}\kappa_y$.
kky = kkx
KX,KY = np.meshgrid(kkx,kky)
kkr = np.sqrt(KX**2+KY**2)
print("check:",np.allclose(kr,kkr))
Eso es todo por hoy.