Search code examples
pythonnumpyfftconvolution

fftshift or convolve from NumPy module incorrect operation

I'm trying to show in Python equivalence of two signals multiplication in time domain to convolution of them in frequency domain.

Thanks to Baddioes in this post, it was shown that such equivalence can be shown by using fftshift function, but the example was done in Matlab. I expected to get absolutely the same result in Python, but it is not the same, the code is below:

f = 2e6
t = np.linspace(0,1,64)

x = np.sin(2*np.pi*f*t)
y = np.sin(2*np.pi*0.5*f*t)

# multiplication in time domain
signal_xy = x*y

# Fourier of multiplied signals and individual signal
FFT_xy = np.fft.fft(signal_xy)
FFT_x = np.fft.fft(x)
FFT_y = np.fft.fft(y)

# plot the result
plt.plot(abs(FFT_xy))
plt.plot(np.fft.fftshift(abs(np.convolve(np.fft.fftshift(FFT_x),np.fft.fftshift(FFT_y), 'same')))/len(FFT_x), linestyle='--')

enter image description here

I also tried to use fftshift and convolve from scipy module, but the result is absolutely the same.

Does fftshift and convolve functions work correctly in Python?

Solution

  • (This is a heavily-edited answer, so I won't be offended if those who upticked it withdraw their vote.)

    Original Answer (more general one to follow)

    You can control the shifting and do it yourself using numpy.roll(). (Note the reversable shifting distances n//2-1 and n//2+1 in the below). Unfortunately, this is ad-hoc (and doesn't actually work for other values of n).

    n=len(x)
plt.plot(abs(FFT_xy))
plt.plot(np.roll(abs(np.convolve(np.roll(FFT_x,n//2-1),np.roll(FFT_y,n//2-1), 'same')),n//2+1)/n, linestyle='--')

    What I think is the real problem

    The convolution theorem ("FFT of convolution = product of FFTs") only holds if the x, y arrays have periodic extensions. However, numpy.convolve doesn't do that - it zero pads as below. For arrays of equal length this means that the convolution is only actually consistent with the periodic version at a single point.

    enter image description here

    A simple solution is to wrap the np.convolve() function with one which first doubles the length of the first array, then only takes the N to 2N-1 elements, which should be consistent with periodic arrays.

    enter image description here

    Full code:

    import numpy as np
import matplotlib.pyplot as plt

def conv( A, B ):
    N = len( A )
    AA = np.concatenate( ( A, A ) )
    return np.convolve( AA, B, 'full' )[N:2*N]

f = 2e6
t = np.linspace(0,1,64)

x = np.sin(2*np.pi*f*t)
y = np.sin(2*np.pi*0.5*f*t)

# multiplication in time domain
signal_xy = x*y

# Fourier of multiplied signals and individual signal
FFT_xy = np.fft.fft(signal_xy)
FFT_x = np.fft.fft(x)
FFT_y = np.fft.fft(y)

# plot the result
n = len( x )

# Forward plot ("Product of FFTs equals FFT of convolution")
plt.plot(abs(FFT_x*FFT_y))
plt.plot( abs(np.fft.fft( conv(x,y) ) ), linestyle='--')

# Backward plot ("FFT of product equals convolution of FFTs / N") - this was what was wanted by the OP
#plt.plot( abs(FFT_xy))
#plt.plot( abs( conv(FFT_x,FFT_y) ) / n, linestyle='--')

plt.show()

    In the forward direction ("FFT of convolution = product of FFTs") we get

    enter image description here

    In the backward direction (which actually what the OP wanted and can be found by switching the commented/uncommented lines in the code) this gives the following. Note that there is no need to shift the FFTs only to shift them back again.

    enter image description here

    It is possible that Matlab may already work with implied periodicity when it does the convolution, but I'm not in a position to check.