How can I use FFT to find convolution function by python?

Question

I have output and input of a system, As I know f=conv(h,g) that h is convolution function in the FFT we can write F=H*G.

So Can I find H by : H=F/G or no?

I try to write a code like this:

    # -*- coding: utf-8 -*-
"""
Created on Wed Dec 15 09:46:52 2021

@author: 20210113-1
"""

import numpy as np
import scipy.signal
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq,ifft

# Define heaviside function
H = lambda x:np.sin(8*x)+np.cos(2*x)
#define gaussian
gauss = lambda x, sig: np.exp(-( (x)/float(sig))**2 )+3

X = np.linspace(-10, 30, num=1000)
X2 = np.linspace(-10,30, num=1000)

# convolute a heaviside with a gaussian
# H_c = np.convolve( H(X),  gauss(X, 1),  mode="same"  )
input_pulse=H(X)

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


output_pulse=H_c 
# deconvolute a the result
H_dc, er = scipy.signal.deconvolve(H_c, gauss(X2, 1) )



SAMPLE_RATE = 20 # Hertz
DURATION = 5  # Seconds
N =len(input_pulse)

yf = fft(input_pulse)
xf = fftfreq(N, 1 / SAMPLE_RATE)

SAMPLE_RATE = 20 # Hertz
DURATION = 5  # Seconds
N =len(output_pulse)

yff = fft(output_pulse)
xff = fftfreq(N, 1 / SAMPLE_RATE)



#### Plot #### 
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="Gauss filter", lw=3 )
ax[2].plot( H_c/H_c.max(), color="#325cab", label="convoluted" ,  lw=3 ) 
ax[3].plot( H_dc,          color="#ab4232", label="deconvoluted", lw=3 ) 
ax[4].plot( xf, np.abs(yf),         color="#ab4232", label="FFT input", lw=3 ) 
ax[5].plot(xff, np.abs(yff),         color="#ab4232", label="FFT output", lw=3 ) 

for i in range(len(ax)):
    # ax[i].set_xlim([0, len(X)])
    # ax[i].set_ylim([-0.07, 1.2])
    ax[i].legend(loc=4)
plt.show()

impulse=ifft(fft(output_pulse)/fft(input_pulse))

but I can not use it in particular way because return this error:

ValueError: operands could not be broadcast together with shapes (1999,) (1000,)

So I have two questions:

1- can I calculate convolution function by this method?

2- If I can , How can I implement in the code?

score 0 · Answer 1 · answered Dec 20 '21 at 09:42

You could pad the two signals to a common size, this can be achieved by specifying the size of the FFT to be applied.

impulse=ifft(fft(output_pulse, 2048)/fft(input_pulse, 2048))[:len(X)]

Then you remove the extra padding

For your particular example it worked well but this is not always the case. Deconvolution is not easy to be solved exactly unless you have more than one sample per symbol.

The scipy deconvolve method computes a convolution by polynomial division giving quotient and residual, the norm or the residual may be much larger than the norm of the quotient, if you are lucky the norm of the residual will be small and the quotient gives a good approximation of the deconvolution.

This fft deconvolution works if the signal is periodic, you concatenate the signal with itself infinitely many times. When we do this padding we are telling it to compute on a periodic signal with a "silence" between two consecutive copies.

How can I use FFT to find convolution function by python?

1 Answers1