None convolucion_1D

Repaso de convolución en 1D

Métodos potenciales de prospección, FCAG, 2023.

Vemos ejemplos de como implementar la convolución lineal $\ast$ y la correlación $\star$ utilizando la TDF. Tomamos dos secuencias con las que realizaremos los cálculos:

In [24]:
import numpy as np
from numpy.fft import fft,ifft
import matplotlib.pyplot as plt
plt.xkcd()
plt.rcParams.update({'font.size': 16}) # tamaño de fuente
In [25]:
a=np.array([1, 2., 1, 0, 0])/4.
b=np.array([3., 3, 3, 3, 3])
M = len(a)
N = len(b)

Convolución lineal 1D $\ast$

In [26]:
conv     = np.convolve(a,b)
conv_fft = np.real(ifft(fft(a,N+M-1)*fft(b,N+M-1)))

#print(conv)
#print(conv_fft)
print(np.allclose(conv,conv_fft))

True

Si centramos el operador, obtenemos (por el teorema del corrimiento de fase) el resultado de la convolución lineal comenzando desde la muestra en el origen:

In [27]:
ac = np.array([0.5,0.25,0,0,0,0,0,0,0.25])
conv_fft_c = np.real(ifft(fft(ac)*fft(b,N+M-1)))

print("operador sin centrar:")
print(conv_fft)

print("operador centrado:")
print(conv_fft_c)

operador sin centrar:
[7.50000000e-01 2.25000000e+00 3.00000000e+00 3.00000000e+00
 3.00000000e+00 2.25000000e+00 7.50000000e-01 1.97372982e-16
 3.94745964e-16]
operador centrado:
[2.25000000e+00 3.00000000e+00 3.00000000e+00 3.00000000e+00
 2.25000000e+00 7.50000000e-01 0.00000000e+00 1.97372982e-16
 7.50000000e-01]
In [28]:
# Fig. 
cmap="inferno"

av = conv_fft.reshape(1,N+N-1) # para graficar
bv = conv_fft_c.reshape(1,N+M-1) # para graficar

vmin=-np.max(conv_fft)
vmax =np.max(conv_fft)
extent = [0,av.shape[1],0,av.shape[0]]
extent = [0,N+M-1,0,1]


plt.figure(figsize=(15,5))
plt.subplot(1,2,1)
plt.title("$\mathbf{a*b}$")
plt.imshow(av,interpolation="None",cmap=cmap,vmin=vmin,vmax=vmax,extent=extent);
plt.grid(color="black",linestyle='-',linewidth=1.5)
plt.xticks(np.arange(N+M-1),[]);plt.yticks([])
plt.subplot(1,2,2)
plt.title("$\mathbf{a_c*b}$")
plt.imshow(bv,interpolation="None",cmap=cmap,vmin=vmin,vmax=vmax,extent=extent);
#plt.xticks([]);plt.yticks([])
plt.grid(color="black",linestyle='-',linewidth=1.5)
plt.xticks(np.arange(N+M-1),[]);plt.yticks([])
plt.show()

Correlación cruzada lineal 1D $\star$

In [29]:
# Correlación cruzada lineal
corr           = np.convolve(np.flip(a),b)
corr_another   = np.correlate(b,a,"full")
corr_fft       = np.real(ifft(  fft(np.flip(a),N+M-1) * fft(b,N+M-1) ))

print(np.allclose(corr,corr_fft))

#print(corr)
#print(corr_another)
#print(corr_fft)

True

Convolución circular 1D $\circledast$

Tomemos el ejemplo del operador de promediación móvil pesado: $f=(1/4,\underline{1/2},1/4)$.

Si queremos hacer la convolución circular por medio de la TDF, debemos tomar el recaudo de calcular la fase correcta del operador. Para ello, como recordamos de la cátedra de Análisis de Señales en Geofísica, reescribimos $f_c=(\underline{1/2},1/4,1/4)$, ya que la TDF implementada en NumPy comienza desde el índice 0 y los vectores se consideran cíclicos.

De esta manera al utilizar numpy.fft.fft obtenemos la TDF correcta para el operador, que es $F=(\underline{1},1/4,1/4)$. Caso contrario, obtendremos la TDF con un corrimiento de fase, ya que calculariamos la TDF de la secuencia $(\underline{1/4},1/2,1/4)$.

In [30]:
s = np.array([1.,2,3,4,5])
f  = np.array([1,2,1])/ 4.  # Este arreglo se interpreta como si la muestra central estaría al principio
fc = np.array([2,1,1])/4.   # Ahora empieza en la muestra a "t=0", mandamos las muestras t<0 al final
F  = np.fft.fft(f)
FC = np.fft.fft(fc)
In [31]:
N  = len(s)
M  = len(f)
s  = s.reshape(1,N)
f  = f.reshape(1,M)
fc = fc.reshape(1,M)
F  = F.reshape(1,M)
FC = FC.reshape(1,M)
In [32]:
# Fig. operador con "stem"
plt.stem(np.arange(M),f.flatten(),use_line_collection=True,basefmt="")
plt.title("$\mathbf{f}$")
plt.xticks([0,1,2],[-1,0,+1])
plt.yticks([])
plt.xlabel("n")
plt.show()
In [33]:
# Fig. operador sin centrar

cmap="seismic"
cmap="inferno"

vmin=-np.max(f)
vmax =np.max(f)
vmins = -np.max(np.abs(F))
vmaxs = np.max(np.abs(F))
extent = [0,f.shape[1],0,f.shape[0]]


plt.figure(figsize=(15,5))
plt.subplot(1,3,1)
plt.title("$\mathbf{f}$")
plt.imshow(f,interpolation="None",cmap=cmap,vmin=vmin,vmax=vmax,extent=extent);
plt.grid(color="black",linestyle='-',linewidth=1.5)
plt.xticks(np.arange(M),[]);plt.yticks([])
plt.subplot(1,3,2)
plt.title("$\Re\{\mathbf{F}\}$")
plt.imshow(np.real(F),interpolation="None",cmap="bone",vmin=vmins,vmax=vmaxs,extent=extent);
#plt.xticks([]);plt.yticks([])
plt.grid(color="black",linestyle='-',linewidth=1.5)
plt.xticks(np.arange(M),[]);plt.yticks([])
plt.subplot(1,3,3)
plt.title("$\Im\{\mathbf{F}\}$")
plt.imshow(np.imag(F),interpolation="None",cmap="bone",vmin=vmins,vmax=vmaxs,extent=extent);
#plt.xticks([]);plt.yticks([])
plt.grid(color="black",linestyle='-',linewidth=1.5)
plt.xticks(np.arange(M),[]);plt.yticks([])
plt.show()

Como observamos, la TDF no se corresponde con la esperada para el operador $f$. Si ahora calculamos la TDF para el operador con la fase correcta, entonces observamos lo siguiente:

In [34]:
# Fig. operador centrado
vmins = 0#-2*np.max(np.abs(FC))
vmaxs = np.max(np.abs(FC))

plt.figure(figsize=(15,5))
plt.subplot(1,3,1)
plt.title("$\mathbf{f_c}$")
plt.imshow(fc,interpolation="None",cmap=cmap,vmin=vmin,vmax=vmax,extent=extent);
plt.grid(color="black",linestyle='-',linewidth=1.5)
plt.xticks(np.arange(M),[]);plt.yticks([])
plt.subplot(1,3,2)
plt.title("$\Re\{\mathbf{F}\}$")
plt.imshow(np.real(FC),interpolation="None",cmap="bone",vmin=vmins,vmax=vmaxs,extent=extent);
#plt.xticks([]);plt.yticks([])
plt.grid(color="black",linestyle='-',linewidth=1.5)
plt.xticks(np.arange(M),[]);plt.yticks([])
plt.subplot(1,3,3)
plt.title("$\Im\{\mathbf{F}\}$")
plt.imshow(np.imag(FC),interpolation="None",cmap="bone",vmin=vmins,vmax=vmaxs,extent=extent);
#plt.xticks([]);plt.yticks([])
plt.grid(color="black",linestyle='-',linewidth=1.5)
plt.xticks(np.arange(M),[]);plt.yticks([])
plt.show()

Ahora si, la fase es la deseada.

El mismo concepto lo aplicaremos cuando utilicemos operadores en dos dimensiones (2D) para procesar mapas gravimétricos o magnéticos.

Eso es todo por hoy.