3
$\begingroup$

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

enter image description here

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. enter image description here

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

enter image description here

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

# %%

$\endgroup$
0

1 Answer 1

3
$\begingroup$

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

$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$
    – xaberus
    Jan 13, 2023 at 6:23
  • $\begingroup$ @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$. $\endgroup$
    – ZR Han
    Jan 13, 2023 at 7:14
  • $\begingroup$ Yes! I bungled up the windows size in the phase factor somewhere in my calculation. Thank you very much for finding this! This is my first question here, should I edit the question to reflect the solution & mention your answer or is it better to leave it as-is? $\endgroup$
    – xaberus
    Jan 13, 2023 at 17:32
  • $\begingroup$ @xaberus Both are fine. Glad to help you out. $\endgroup$
    – ZR Han
    Jan 16, 2023 at 1:50

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.