Ejemplo sencillo para interpretar el espectro de amplitud de 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
Dato¶
Generamos una señal sintética sencilla para interpretar. En este caso, tomamos muestras de una onda plana
$$f(x,y)=A\cos\left(\mathbf{k}\cdot\mathbf{x}\right)=A\cos\left(k_{x}x+k_{y}y\right).$$Calculamos también su espectro de amplitud por medio de la transformada discreta de Fourier 2D.
El vector número de onda tiene magnitud $||\mathbf{k}||_2=k_r=\sqrt{k_x^2+k_y^2}$ y ángulo $\phi$, tal que $k_x= k_r \cos{\phi}$ y $k_y= k_r \sin{\phi}$ .
# definimos la señal y su espectro de amplitud
def f(x,y,amp=1.0):
f = amp*np.cos(2.0*np.pi*(x+y))
sf = np.abs(np.fft.fftshift(np.fft.fft2(f)))
return f,sf
# los puntos de observación como matriz
nx, ny = (2**7,2**7)
xx = np.linspace(-0.5, 0.5, nx)
yy = np.linspace(-0.5, 0.5, ny)
x, y = np.meshgrid(xx, yy)
# la señal discretizada y su espectro de amplitud para |k_r|= 2π y ϕ=45
kr = 1. # el 2π está dentro de la subrutina f()
angle = 45.
arg = angle*np.pi/180.
kx = kr*np.cos(arg)
ky = kr*np.sin(arg)
s,sf = f(kx*x,ky*y)
Gráficos¶
Graficamos la señal y su espectro de amplitud para $|\kappa_r|=2\pi$ y $\phi=45^\circ$. Utilizamos numpy.fft.fftshift para centrar la TDF 2D a los fines interpretativos.
cmap1='coolwarm'
cmap2='gray'
plt.figure(1,figsize=(16,8))
plt.subplot(121)
plt.title('$s(x,y)=A\,\cos(k_x\,x+k_y\,y)$')
plt.xticks([])
plt.yticks([])
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.imshow(s,origin='lower',cmap=cmap1,interpolation=None)
plt.subplot(122)
plt.title('$|S(k_x,k_y)|$')
plt.xticks([])
plt.yticks([])
plt.xlabel('$k_x=|k|\cos(\phi)$')
plt.ylabel('$k_y=|k|\sin(\phi)$')
plt.imshow(sf,origin='lower',cmap=cmap2,interpolation=None)
plt.show()
Visualización interactiva¶
Graficamos de forma interactiva para distintos valores de la amplitud $a$ y el vector número de onda $\mathbf{k}=(k_x,k_y)$. El vector número de onda tiene magnitud $||\mathbf{k}||_2=k_r=\sqrt{k_x^2+k_y^2}$ y ángulo $\phi$, tal que $k_x= k_r \cos{\phi}$ y $k_y= k_r \sin{\phi}$ .
from ipywidgets import interact,widgets # widgets para interacción
# interact dentro de la notebook
#%matplotlib notebook
def plotter(amp=1.0,kr=1.0,angle=0.0):
#from matplotlib.colors import LogNorm
# obtengo kx y ky de |k| y el ángulo
arg = angle*np.pi/180.
kx = kr*np.cos(arg)
ky = kr*np.sin(arg)
# calculo el dato y su espectro de amplitud
data,fdata = f(kx*x,ky*y,amp)
# grafico
plt.figure(2,figsize=(16,8))
plt.subplot(121)
plt.title('$s(x,y)=A\,\cos(k_x\,x+k_y\,y)$')
plt.xticks([])
plt.yticks([])
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.imshow(data,origin='lower',cmap=cmap1,interpolation=None)
plt.subplot(122)
plt.title('$|S(k_x,k_y)|$')
plt.xticks([])
plt.yticks([])
plt.xlabel('$k_x=|k|\cos(\phi)$')
plt.ylabel('$k_y=|k|\sin(\phi)$')
plt.imshow(fdata,origin='lower',cmap=cmap2,interpolation=None)#, norm=LogNorm(vmin=0.1))
plt.show()
# gráfico interactivo por defecto (etiquetas con el mismo nombre que la variable)
#interact(update,kr=(0.0,fN),angle=(0.0,180.0,10.0),amp=(1.0,3.0));
fN=0.5*np.min([nx,ny]) # nx == fS, nx/2 == fN
# parámetros explÃcitos para los sliders y descripción (las etiquetas pueden editarse)
amp_slider = widgets.FloatSlider(min=1.,max=3.,step=0.5,value=1.0,description='$A$')
angle_slider = widgets.FloatSlider(min=0.,max=180.,step=10.0,value=45.0,description="$\phi (\circ)$ ")
kr_slider = widgets.FloatSlider(min=0.,max=0.5*fN,step=fN/10.,value=1,description='$|k|$')
# gráfico interactivo con etiquetas
interact(plotter,kr=kr_slider,angle=angle_slider,amp=amp_slider);
Ejemplo con imagen satelital¶
Analizamos el espectro de amplitud de una imagen satelital contaminada por ruido periódico, la cual es analizada en el siguiente blog. El espectro de amplitud ha sido calculado en escala logarÃtmica y con zero padding para facilitar su análisis:
epsilon = 1e-15
IMG = np.fftshift(np.fft.fft2(img,(1024,1024)))
IMG = np.log(IMG+epsilon)
¿Notamos en el espectro de amplitud donde reside mayoritariamente el ruido periódico?
from matplotlib import image
img = image.imread("../../data/aerialpompeiiperiodic.jpg")
plt.figure(figsize=(20,8))
plt.imshow(img,cmap="gray")
plt.xticks([])
plt.yticks([])
plt.show()
plt.figure(figsize=(8,8))
plt.imshow(np.log(np.abs(np.fft.fftshift(np.fft.fft2(img,(1024,1024))))+1e-15),cmap="gray")
plt.xticks([])
plt.yticks([])
plt.show()
Eso es todo por hoy.