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", ) # %%