# How to get DFT spectral leakage from convolution theorem?

I have an issue, where my numerically calculated leakage from a DFT of a simple cosine does not match the theoretical prediction from the convolution theorem. I will try to present the example using technical units as opposed to the usual unit-less representation, so that the different Fourier transform can be more easily compared to each other. I am more or less closely following this blogpost, that has a nice explanation of the math involved.

From my understanding, the DFT of length $$N$$ of a some continuous, band-limited function $$f$$ sampled at a frequency $$\nu_\mathrm{s}$$ can be expressed as: $$\operatorname{DFT}_N\left[f\right](l)=\sum_{l=0}^{N-1}f_{l}e^{-i\,\frac{2\pi}{N}kl},$$ where $$f_l=f(\frac{l}{\nu_\mathrm{s}})$$. By using the convolution theorem the DFT can also be expressed using the three implicitly involved functions and distributions: $$\operatorname{DFT}_N\left[f\right](l)=\mathcal{F}\left[f \cdot s_{\nu_\mathrm{s}} \cdot W_N \right](\nu_l) = \left(\mathcal{F}\left[f\right] \ast \mathcal{F}\left[s_{\nu_\mathrm{s}} \cdot W_N\right]\right)(\nu_l) ,$$ where the continuous Fourier transform is defined as $$\mathcal{F}\left[f\right](\nu)=\int_{-\infty}^{\infty}f(t)\,e^{-\,2\pi\nu\,t}\,\mathrm{d}t,$$ the Dirac comb is defined as $$s_{\nu_\mathrm{s}}(t)=\sum_{k\in\mathbb{Z}}f_{k}\,\delta(t-\frac{k}{\nu_\mathrm{s}}),$$ and the rectangular window function $$W_N(t)$$ is nonzero only for the samples $$0$$ to $$N-1$$, and $$\nu_l = l \frac{\nu_\mathrm{s}}{N}.$$

The Fourier transform of the window function is the regular sinc $$\mathcal{F}\left[W_N\right](\nu) \propto \frac{\sin\left(\pi\frac{\nu}{\nu_\mathrm{s}} N\right)}{\pi\nu},$$ and its convolution with the Dirac comb gives the Dirichlet kernel (aliased periodic sinc function): $$\mathcal{F}\left[s_{\nu_\mathrm{s}} \cdot W_N\right](\nu)= \exp\left(- i \pi \frac{\nu}{\nu_\mathrm{s}} N\right) \frac{ \sin\left(\pi \frac{\nu}{\nu_\mathrm{s}} N \right) }{ \sin\left(\pi\frac{\nu}{\nu_\mathrm{s}}\right) }.$$

Now, let us assume that $$f=\cos\left( 2 \pi \mu_0 \right)$$. Then $$\mathcal{F}\left[f\right](\nu)= \frac{1}{2}\left(\delta(\nu-\mu_0) + \delta(\nu+\mu_0)\right).$$ Finally the complete convolution from above is then $$\mathcal{C}(\nu)= \frac{1}{2}\left(\exp\left(-i\pi\frac{\nu-\mu_{0}}{\nu_{\mathrm{s}}}N\right)\frac{\sin\left(\pi\frac{\nu-\mu_{0}}{\nu_{\mathrm{s}}}N\right)}{\sin\left(\pi\frac{\nu-\mu_{0}}{\nu_{\mathrm{s}}}\right)}+\exp\left(-i\pi\frac{\nu+\mu_{0}}{\nu_{\mathrm{s}}}N\right)\frac{\sin\left(\pi\frac{\nu+\mu_{0}}{\nu_{\mathrm{s}}}N\right)}{\sin\left(\pi\frac{\nu+\mu_{0}}{\nu_{\mathrm{s}}}\right)}\right),$$ that it, two Dirichlet kernels shifted to either the positive or negative frequency $$\mu_0$$.

Unfortunately, this convolution result does not in general match the result of the DFT, i.e., $$\mathcal{C}(\nu_l) \ne \operatorname{DFT}_N\left[f\right](l)$$ and this mismatch is particularly visible for small DFT sizes. Usually the values around the major peak and $$\pm$$ two bins match quite well, while all other bins do not match.

A "spectrum" that does match can be obtained by interpolating the DFT using the generalized DFT (or zero extension of the input): $$\mathcal{G}_{N}\left[f\right](\nu)=\sum_{l=0}^{N-1}f_{l}e^{-i\,\nu l}.$$

1. What am I missing? Did I make an error in the calculation?
2. Does this phenomenon have a name?

Thank you for reading such a long question :-)

I used the script below to generate the attached images. Uncomment either the "interpolation" or "convolution" in the elements set to switch between the two different spectra.

# %%
import numpy as np
import warnings

import matplotlib.pyplot as plt

warnings.resetwarnings()

def gdft(x, samples=None, fmin=None, fmax=None, fshift=None):
"""Calculate Generalized DFT."""
x = np.asanyarray(x)
n = x.size

if samples is None:
samples = n // 2 + 1
if fmin is None:
fmin = 0
if fmax is None:
fmax = n // 2
if fshift is None:
fshift = 0

f = np.linspace(fmin + fshift, fmax + fshift, samples)
j = np.arange(n)[:, np.newaxis]
z = 2 * np.pi * (f - fshift) * j / n
roots = np.exp(-1j * z)
return f, np.sum(x[:, np.newaxis] * roots, axis=0)

def dirichlet_dtft(f: float, fs: float, period: float = 1):
"""Calculate the DTFT of the rectangular window."""
f = np.asanyarray(f)
pif = np.pi * f
spif = np.sin(pif / fs)

def inner(pif):
spif = np.sin(pif / fs)
sinc = np.sin(pif * period) / spif / fs
return sinc

return np.piecewise(
pif,
np.isclose(spif, 0),
[
1.0 * period,
inner,
],
) * np.exp(-1j * pif * period)

def leaky_cosine_dtft(f: float, f0: float, fs: float, period: float = 1.0):
"""Calculate the sampled finite bandwidth spectrum of a clipped cosine."""
vw = dirichlet_dtft(f - f0, fs, period=period)
vw += dirichlet_dtft(f + f0, fs, period=period)
vw /= 2
return vw

# %%

elements = set(
[
# "interpolation",
"convolution",
]
)

# sample rate
fs = 1000 * 4
fft_size = 32
# DFT frequency spacing
df = fs / fft_size
# time
t = np.arange(0, fft_size) / fs
interpolate = 3001
# cosine frequency
f0 = 550

# sampled cosine
v = np.cos(2 * np.pi * f0 * t)

# sampled cosine DFT
vs = np.fft.fft(v) / fft_size
f = np.fft.fftfreq(fft_size) * fs

# interpolated cosine
ti = np.linspace(
-1 / fs * fft_size / 10,
1 / df + 1 / fs * fft_size / 10,
interpolate,
)
vi = np.cos(2 * np.pi * f0 * ti)

# interpolated "sinc"
fi, si = gdft(v, interpolate, fmin=-fft_size // 2, fmax=fft_size // 2)
fi *= df
si /= fft_size

period = fft_size / fs
fc = np.linspace(-fs / 2, fs / 2, interpolate)
sc = leaky_cosine_dtft(fc, f0, fs, period) / period

fig = plt.figure(figsize=(6, 10))
nrows=3,
ncols=1,
width_ratios=[1],
height_ratios=[0.68, 1, 1],
)

bx.tick_params("x", labelbottom=False)

ax.plot(t, v, color="black", ls="", marker=".")
ax.plot(ti, vi, color="black", ls="-", marker="", lw=2)
ax.axvspan(0, 1 / df, alpha=0.2)
ax.axvline(0, marker="", color="black", lw=1)
ax.axvline(1 / df, marker="", color="black", lw=1)
ax.axhline(0, marker="", color="black", lw=1)

bx.plot(
f,
np.real(vs),
color="blue",
ls="",
marker=".",
label=r"$$\Re[\operatorname{DFT}]$$",
)
bx.plot(
f,
np.imag(vs),
color="red",
ls="",
marker=".",
label=r"$$\Im[\operatorname{DFT}]$$",
)
if "interpolation" in elements:
bx.plot(fi, np.real(si), color="green", lw=1, label=r"$$\Re[\mathcal{G}]$$")
bx.plot(fi, np.imag(si), color="orange", lw=1, label=r"$$\Re[\mathcal{G}]$$")
if "convolution" in elements:
bx.plot(fc, np.real(sc), color="blue", lw=1, label=r"$$\Re[\mathcal{C}]$$")
bx.plot(fc, np.imag(sc), color="red", lw=1, label=r"$$\Im[\mathcal{C}]$$")
bx.axhline(0, marker="", color="black", lw=0.5, ls="dashed")

cx.plot(
f,
np.abs(vs),
marker=".",
color="black",
ls="",
label=r"$$\left|\operatorname{DFT}]\right|$$",
)
if "interpolation" in elements:
cx.plot(
fi,
np.abs(si),
marker="",
color="black",
ls="-",
label=r"$$\left|\mathcal{G}\right|$$",
)
if "convolution" in elements:
cx.plot(
fc,
np.abs(sc),
color="orange",
lw=1,
label=r"$$\left|\mathcal{C}\right|$$",
)
cx.axhline(0, marker="", color="black", lw=0.5, ls="dashed")

# cx.set_yscale("log")
# cx.set_ylim(ymin=1e-5, ymax=2)

for i in (0, 1, 2, 3, 4, -1, -2, -3):
dfx = fs / (fft_size + 0.3)
fzp = f0 + i * dfx
fzm = -f0 - i * dfx
for axi in (bx, cx):
for fzpos, color in ((fzp, "red"), (fzm, "blue")):
axi.axvline(
fzpos,
marker="",
color=color,
ls="dashed",
lw=0.3,
zorder=-100,
)

bx.legend()
cx.legend()

ax.set_xlabel("time $$t$$ (s)")
ax.set_ylabel("signal $$u$$ (V)")

bx.set_ylabel(r"amplitude $$\mathcal{F}$$ (V)")
cx.set_ylabel(r"magnitude $$\left|\mathcal{F}\right|$$ (V)")
cx.set_xlabel(r"frequency $$\nu$$ (Hz)")

fig.savefig(
f"fft-test.pdf",
bbox_inches="tight",
)

# %%



There is something wrong with your derivation of the Fourier transform for the rectangular window. A regular sinc corresponds to a rectangular window that is symmetric about $$t=0$$. We need a causal window function which corresponds to a phase-shifted sinc function.
Yet another simple way to derive the Fourier transform of the sampled rectangular window is to calculate its DTFT directly. In the following, I'll use a more common notation in the field of DSP, that is, $$n$$ denotes the discrete time index, $$\omega$$ is the normalized angular frequency.
\begin{aligned} \mathrm{DTFT}\{W_N[n]\} &= \sum_{n=-\infty}^{\infty} W_N[n] e^{-j\omega n} = \sum_{n=0}^{N-1} e^{-j\omega n} \\ & = \frac{1-e^{-j\omega N}}{1-e^{-j\omega}} = \frac{e^{-j\omega N/2}(e^{j\omega N/2} - e^{-j\omega N/2})}{e^{-j\omega/2}(e^{j\omega/2}-e^{-j\omega/2})} \\ & = e^{-j\omega\frac{N-1}{2}} \frac{\sin(\omega N /2)}{\sin(\omega/2)} \end{aligned}
• I think you are right - thank you very much, this was driving me insane. The only difference seems to be the $N-1$ instead of $N$ in the exponential. I just tested your formula, and it indeed works. In my derivation I assumed that the window can be arbitrary (as the DFT is periodic anyway). Why does it matter? Can you point me to some literature about it, please? Also, should I edit my question to reflect your answer? Jan 13 at 6:23
• @xaberus It seems that you calculate the Fourier transform of the window function by integral from $0$ to $NT$ where $T=1/f_s$ is the sampling period and $f_s$ is the sampling frequency, according to the definition of $W(t) = 1, 0\leq t < NT$. However after sampling $n=N$ is unreachable. The integral should be from $0$ to $(N-1)T$. Jan 13 at 7:14