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}. $$
- What am I missing? Did I make an error in the calculation?
- Does this phenomenon have a name?
- 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",
)
# %%