How to get DFT spectral leakage from convolution theorem?

Question

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?
3. Is there any literature I can consult about this?

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))
gs = fig.add_gridspec(
    nrows=3,
    ncols=1,
    width_ratios=[1],
    height_ratios=[0.68, 1, 1],
)

ax = fig.add_subplot(gs[0, 0])
bx = fig.add_subplot(gs[1, 0])
cx = fig.add_subplot(gs[2, 0], sharex=bx)
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",
)

# %%

Answer 1

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} $$