¿Cómo puedo usar la transformada rápida de Fourier (FFT, por sus siglas en inglés) para encontrar la función de convolución utilizando Python?

21 Jul, 2023 Programación 0

Tengo la salida y la entrada de un sistema. Como sé que f=conv(h,g), donde h es una función de convolución en la FFT, podemos escribir F=H*G.

Entonces, ¿puedo encontrar H por medio de: H=F/G o no?

Intento escribir un código como este:

    # -<em>- coding: utf-8 -</em>-
"""
Created on Wed Dec 15 09:46:52 2021

<p>@author: 20210113-1
"""</p>

<p>import numpy as np
import scipy.signal
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq,ifft</p>

<h1>Define la función heaviside</h1>

<p>H = lambda x:np.sin(8<em>x)+np.cos(2</em>x)</p>

<h1>Define la función gaussiana</h1>

<p>gauss = lambda x, sig: np.exp(-( (x)/float(sig))**2 )+3</p>

<p>X = np.linspace(-10, 30, num=1000)
X2 = np.linspace(-10,30, num=1000)</p>

<h1>Convoluciona una función heaviside con una función gaussiana</h1>

<h1>H_c = np.convolve( H(X),  gauss(X, 1),  mode="same"  )</h1>

<p>input_pulse=H(X)</p>

<p>H_c = scipy.signal.convolve(H(X),  gauss(X2, 1), mode='full') #/ c.sum()</p>

<p>output<em>_pulse=H</em>_c </p>

<h1>Desconvoluciona el resultado</h1>

<p>H<em>_dc, er = scipy.signal.deconvolve(H</em>_c, gauss(X2, 1) )</p>

<p>SAMPLE<em>_RATE = 20 # Hertz
DURATION = 5  # Seconds
N =len(input</em>_pulse)</p>

<p>yf = fft(input<em>_pulse)
xf = fftfreq(N, 1 / SAMPLE</em>_RATE)</p>

<p>SAMPLE<em>_RATE = 20 # Hertz
DURATION = 5  # Seconds
N =len(output</em>_pulse)</p>

<p>yff = fft(output<em>_pulse)
xff = fftfreq(N, 1 / SAMPLE</em>_RATE)</p>

<h4>Gráfico</h4>

<p>fig , ax = plt.subplots(nrows=6, figsize=(6,7))
ax[0].plot( H(X),          color="#907700", label="Heaviside",    lw=3 ) 
ax[1].plot( gauss(X, 1),  color="#907700", label="Filtro de Gauss", lw=3 )
ax[2].plot( H<em>_c/H</em>_c.max(), color="#325cab", label="Convolucionado" ,  lw=3 ) 
ax[3].plot( H_dc,          color="#ab4232", label="Desconvolucionado", lw=3 ) 
ax[4].plot( xf, np.abs(yf),         color="#ab4232", label="FFT entrada", lw=3 ) 
ax[5].plot(xff, np.abs(yff),         color="#ab4232", label="FFT salida", lw=3 ) </p>

<p>for i in range(len(ax)):
    # ax[i].set<em>xlim([0, len(X)])
    # ax[i].set</em>ylim([-0.07, 1.2])
    ax[i].legend(loc=4)
plt.show()</p>

<p>impulse=ifft(fft(output<em>_pulse)/fft(input</em>_pulse))

pero no puedo usarlo de manera particular porque devuelve este error:

ValueError: las operaciones no se pueden transmitir juntas con formas (1999,) (1000,)

Así que tengo dos preguntas:

1- ¿Puedo calcular la función de convolución utilizando este método?

2- Si puedo, ¿cómo puedo implementarlo en el código?

Answer

davy.ai

July 21, 2023 at 7:38 pm
Respuesta:
1. ¿Puedo calcular la función de convolución mediante este método?
  No, no se puede calcular la función de convolución directamente dividiendo las FFT de las señales de salida e entrada. La división en el dominio de frecuencia solo es matemáticamente válida para la desconvolución.
2. Si puedo hacerlo, ¿cómo puedo implementarlo en el código?
  Si desea calcular la función de convolución ‘H’, puede utilizar la función scipy.signal.deconvolve, que devuelve tanto la señal desconvolucionada como la respuesta del sistema inversa (es decir, la función de convolución ‘H’). Aquí le mostramos cómo puede modificar su código para calcular ‘H’:
```
input_pulse = H(X)
filter_kernel = gauss(X2, 1)
H_c = scipy.signal.convolve(input_pulse, filter_kernel, mode='full')
output_pulse = H_c 
# para obtener H, desconvoluciona el resultado
_, H = scipy.signal.deconvolve(H_c, filter_kernel)
```
Tenga en cuenta que la salida de scipy.signal.deconvolve es una tupla que contiene la señal desconvolucionada y la respuesta del sistema inversa. Solo estamos interesados en el segundo valor devuelto, que es la función de convolución ‘H’.

Comments are closed.

es.davy.ai

¿Tienes una pregunta?

¿Cómo puedo usar la transformada rápida de Fourier (FFT, por sus siglas en inglés) para encontrar la función de convolución utilizando Python?

Answer

davy.ai