2
$\begingroup$

Assuming the convolution in the time domain produces a signal which spectrum is the product of individual spectrums, when the two convoluted signals are sinusoids with different frequencies, the DFT coefficients should be null. I'm trying to verify this property, but my spectrum (last column) is not null.

enter image description here

Assuming the code is correct, how to interpret the result?

import numpy as np
import scipy.fft as sf
import matplotlib.pyplot as plt

# Sampling
Fs = 200
T = 0.5
t = np.arange(0, T, 1/Fs)
N = len(t)

# Frequencies, Hz
f1 = 10
f2 = 50

# Functions
w = lambda f: 2 * np.pi * f
dft = lambda s: sf.fftshift(sf.fft(s))

# Time domain
x1 = 1 * np.cos(w(f1) * t)
x2 = 0.1 * np.cos(w(f2) * t)
ya = x1 + x2
yc = np.convolve(x1, x2, mode='same')

# Frequency domain
ff = sf.fftshift(sf.fftfreq(N, 1/Fs))
X1, X2, Ya, Yc = [dft(s) for s in (x1, x2, ya, yc)]

# Figure
mosaic = np.asarray([[f't{i}' for i in range(4)], [f'f{i}' for i in range(4)]])
nrows = len(mosaic)
ncols = len(mosaic[0])
fig_kw = dict(figsize=(ncols*2.6, nrows*1.5), layout='constrained')
fig, axs = plt.subplot_mosaic(mosaic=mosaic, **fig_kw)

for ax, title in zip(mosaic[0], ('x1(t)', 'x2(t)', 'Sum', 'Convolution')):
    axis = axs[ax]
    axis.set_title(title)
    axis.set_xlabel('time (s)')
for ax in mosaic[1]: axs[ax].set_xlabel('Hz')
fig.text(0.5, 1.04, ha='center', size='large', weight='bold', s='Sum vs. convolution')

for ax, x in zip(mosaic[0], (x1, x2, ya, yc)): axs[ax].plot(t, x)
for ax, X in zip(mosaic[1], (X1, X2, Ya, Yc)): axs[ax].stem(ff, abs(X))
$\endgroup$
2
  • 1
    $\begingroup$ Multiplication of DFTs corresponds to circular convolution. You've computed the linear convolution. $\endgroup$
    – Matt L.
    Commented May 25, 2023 at 11:17
  • 1
    $\begingroup$ You could answer your own question as a reference for future visitors, or otherwise just delete it. $\endgroup$
    – Matt L.
    Commented May 25, 2023 at 14:44

2 Answers 2

3
$\begingroup$

In short:

As corrected by @Matt, the DFT pair convolution-product refers to a circular convolution, not to a linear one.

For Fourier transform, the property convolution FT is the product of individual FT actually holds for continuous-time FT. For discrete-time FT this is not generally true, however this always holds for circular convolution, due to the fact DFT, contrary to FT, is periodic/cyclic (more).


Python

Numpy convolve doesn't provide this option, but scipy.ndimage.convolve1d can be used.

Adjusting the code this way:

import scipy.ndimage as nd
yc = nd.convolve1d(x1, x2, mode='wrap')

gives what is seemingly correct (scales for convolution result and related DFT are scaled by 1e-15 and 1e-13):

enter image description here


I believe using the (circular) convolution this way is not frequent. For long sequences, computing the DFT product and taking the inverse DFT is more efficient than performing the convolution. But I'm leaving this answer for those who, while learning DSP, also fall in the trap of mixing linear and circular convolutions.

$\endgroup$
0
0
$\begingroup$

Disjoint spectra do produce a zero result, but the convolution is circular. However, I encourage against thinking that FFT-based and "direct" (np.convolve) convolutions are incompatible. The DFT remains an invaluable tool for debugging convolutions and designing filters.

Instead of taking FFT of the result of np.convolve, set up FFT convolution in a way that achieves the same result - mode='full' can be replicated directly, and other modes with unpadding. Unpadding is aliasing, and arbitrary unpad amounts are difficult to understand spectrally, but if padding equals input's length (or integer multiple of), we can utilize "DFT unpad property" to understand its effects on spectrum. For a filter's frequency response, this most often isn't necessary and we can just look at its face value spectrum.

Related, there's a framework of thought in traditional signal processing that describes convolution in terms of "delay" and "time aliasing"; that's not all of SP. In time-frequency, and as I'd like to call "informative perspective", they can be called "temporal expansiveness" and "boundary effects", respectively. As such, padding becomes about imposing a "prior", and zero-padding an assumption of silence or acknowledgement of information deficiency - related. Put differently, it's not "wrong" to look at spectrum of a padding (it depends).

Example of extensive DFT-domain filter analysis here. In my line of work, a lot rests on the fact that disjoint spectra in convolution does mean zero output.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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