0
$\begingroup$

It is well-known that convolution in real space corresponds to multiplication in the reciprocal space. Yet, I obtain different results using both methods in the case that I'm going to present. I am sure that I made a mistake, but I really can't see where.

The Python code:

%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt

# Create a signal: a gaussian
sig = np.exp(-(np.linspace(-10, 10, 1024))**2)
# Calculate its FFT transform
F_SIG = np.fft.fft(sig)

# Create an FIR filter from its transform
F_KER = np.empty_like(sig, dtype=np.complex128)
F_KER[:len(F_KER)//2] = np.exp(-1j*np.pi/3)
F_KER[len(F_KER)//2:] = np.exp(1j*np.pi/3)
F_KER[[0, len(F_KER)//2]] = np.cos(np.pi/3)
# Get its inverse FFT transform
kern = np.fft.ifft(F_KER) # is actually pure real value
print("Max real value:\t", np.abs(kern.real).max())
print("Max imaginary value:\t", np.abs(kern.imag).max())
print("Therefore 'kern' can be considered pure real.")
kern = kern.real

# From math books, 
# multiplication in Fourier/Frequency space <-> convolution in real space
multFFT_sig = np.fft.ifft(F_SIG*F_KER).real # Again, pure real, can be checked
conv_sig = np.convolve(sig, kern, "full")


# Plot both of them
plt.figure()
plt.xlim([0,1024])
plt.title("Convolution vs. FFT multiply")
plt.plot(multFFT_sig)
plt.plot(conv_sig)


# Padding problem? Circular convolution problem?
plt.figure()
plt.xlim([0,1024])
plt.title("Original signals")
plt.plot(sig)
plt.plot(kern)
plt.legend(["Orig. signal", "Conv. kernel"])

plt.figure()
plt.xlim([0,2048])
plt.title("Zero-padded signals")
p_sig = np.append([0]*512, np.append(sig, [0]*512))
p_kern = np.append(kern[:512], np.append([0]*1024, kern[512:]))
plt.plot(p_sig)
plt.plot(p_kern)
plt.legend(["pad. signal", "pad. conv. kernel"])

# Calculate convolution padded terms
plt.figure()
plt.xlim([0,2048])
plt.title("Padded convolution vs. FFT multiply")
p_F_SIG = np.fft.fft(p_sig)
p_F_KER = np.fft.fft(p_kern)
p_multFFT_sig = np.fft.ifft(p_F_SIG*p_F_KER).real # Again, pure real, can be checked
p_conv_sig = np.convolve(p_sig, p_kern, "full")
plt.plot(p_multFFT_sig)
plt.plot(p_conv_sig)

The output:

Max real value:  0.551327165647
Max imaginary value:     2.55232033425e-15
Therefore 'kern' can be considered pure real.

enter image description here

I also tested with scipy.signal.fftconvolve and compared with the analytical answer: the FFT multiplication gives the right answer (clearly asymmetrical peak in blue) while the convolution (in real space or using FFT) is wrong.

Why is convolution not working as intended?

$\endgroup$
2
$\begingroup$

Yes, zero-padding is the issue, however not as you would except...

@Soren is only partially right: you can add zeros at the beginning, middle or end of the signal just as you did, but absolutely not as you wish.

What's the general rule?

Look at how the signal repeats with and without zeros and choose the place that distorts the signal the least.

What does it mean?

  • In the case of a real, discrete and aperiodic time-space signal, there is a causal relationship between the value at n and the value at n+1. So all these values should be considered as a bloc, adding zeros in the middle would distort it. Here is an obvious reason from your code, you had a proper Gaussian, not something like this:

enter image description here

  • In the case of your convolution kernel, it's a bit more complicated. What do the last taps (second half) mean?
    • In the aperiodic DTFT (i.e. linear convolution, as opposed to circular), these mean that the convolved values depends of both [n, n+1, n+2, ...] (as expected) but also [len - 1, len - 2, len - 3, ...] (where len corresponds to the end of the array without padding). This "classical" linear convolution can be simulated by padding the end/beginning of the signal.
    • In the periodic DFT (circular convolution), the taps of the second half actually correspond to the values [n-1, n-2, n-3, ...] (not the values by the end of the signal).

np.convolve and bunch try to do a "classical" linear convolution while you seem to want a non-classical one (zeros added in the middle). Here is a comparison: enter image description here

As I understood, the actual kernel you want is the one described in "your kernel".

You can easily get the right kernel by fftshifting:

kern = np.fft.fftshift(kern)

This has the effect of centering your DC (using linear phase shifting) and then you can just pad zeros at the beginning/end as much as you wish and use the standard convolution/filtering tools.


In the end, this was an interesting question: people usually just pad blindly using convolution functions or options from the FFT transform without ever wondering what they are doing. Some even put fftshifts out of habit every time they deal with anything resembling a kernel.

$\endgroup$
0
$\begingroup$

When zero-padding you should insert zeros after the signal that you wish to pad.

Change

p_sig = np.append([0]*512, np.append(sig, [0]*512))
p_kern = np.append(kern[:512], np.append([0]*1024, kern[512:]))

to

p_sig = np.append(sig,np.zeros(1024))
p_kern = np.append(kern,np.zeros(1024))

and you'll get a correct result.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.