Fourier transform of modulus of sum of sines

Question

$$ x(t) = |\cos(\omega_0 t) + \cos(\omega_1 t)| $$

with $\omega_0, \omega_1 > 0$.

Is there a known result for $\mathcal{F}\{x(t)\}$? Derivation not needed but is welcome. Of main interest is the DFT, but I don't reckon we can find it without $\mathcal{F}$.

---

Let $x(t) = |y(t)|$. Simplifying notation, we know $y$ is

$$ 2\cos(.5(A + B))\cos(.5(A - B)) $$

which is $2\pi$-periodic in both $A$ and $B$, so it seems one can reasonably find a periodic windowing function to isolate $d(t) = |y(t)| - y(t)$.

---

WA gives for $|\cos(A) + \cos(B)|$

$$ \frac{2}{\pi} - \frac{4}{\pi} \sum_{k=1}^{\infty} \frac{(-1)^k T_{2k}(\cos(A) + \cos(B))}{-1 + 4k^2} $$

where $T_k$ is the Chebyshev polynomial, 1st kind; one possible representation:

$$ T_n(x) = .5\left[ (x + \sqrt{x^2 - 1})^{-n} + (x + \sqrt{x^2 - 1})^n \right] $$

which means one of the things we need the $\mathcal{F}$ of, for all even $n$, is

$$ \frac{\cos(\omega_0 t) + \cos(\omega_1 t)} { \left((\cos(\omega_0 t) + \cos(\omega_1 t)) + \sqrt{(\cos(\omega_0 t) + \cos(\omega_1 t))^2 -1}\right)^n } $$

but that radicand isn't even $\geq 0$...

Answer 1

$$\begin{align*}\mathscr{F}\left\{x(t)\right\} &= \mathscr{F}\left\{\left|\cos\left(\omega_0t\right)+\cos\left(\omega_1t\right)\right|\right\}\\ \\ &= 2 \mathscr{F}\left\{\left|\cos\left(\dfrac{\omega_0-\omega_1}{2\pi}\pi t\right)\cos\left(\dfrac{\omega_0+\omega_1}{2\pi}\pi t\right)\right|\right\}\\ \\ &= 2 \cdot \mathscr{F}\left\{\left|\cos\left(\dfrac{\omega_0-\omega_1}{2\pi}\pi t\right)\right|\right\}*\mathscr{F}\left\{\left|\cos\left(\dfrac{\omega_0+\omega_1}{2\pi}\pi t\right)\right|\right\}\\ \\ &= 2 \cdot \dfrac{2\pi}{\left|\omega_0-\omega_1\right|}\left[\dfrac{1}{2}\mathrm{III}\left(\dfrac{2\pi}{\omega_0-\omega_1}s\right)\left(\mathrm{sinc}\left[\dfrac{2\pi}{\omega_0-\omega_1}s+\dfrac{1}{2}\right]+\mathrm{sinc}\left[\dfrac{2\pi}{\omega_0-\omega_1}s-\dfrac{1}{2}\right]\right)\right]\\ &\quad* \;\;\dfrac{2\pi}{\left|\omega_0+\omega_1\right|}\left[\dfrac{1}{2}\mathrm{III}\left(\dfrac{2\pi}{\omega_0+\omega_1}s\right)\left(\mathrm{sinc}\left[\dfrac{2\pi}{\omega_0+\omega_1}s+\dfrac{1}{2}\right]+\mathrm{sinc}\left[\dfrac{2\pi}{\omega_0+\omega_1}s-\dfrac{1}{2}\right]\right)\right]\\ \\ &= 2 \cdot \left[\dfrac{1}{2}\sum_{n=-\infty}^\infty\left(\mathrm{sinc}\left[n+\dfrac{1}{2}\right]+\mathrm{sinc}\left[n-\dfrac{1}{2}\right]\right)\delta\left(s-\dfrac{\omega_0-\omega_1}{2\pi}n\right)\right]\\ &\quad* \;\;\left[\dfrac{1}{2}\sum_{n=-\infty}^\infty\left(\mathrm{sinc}\left[n+\dfrac{1}{2}\right]+\mathrm{sinc}\left[n-\dfrac{1}{2}\right]\right)\delta\left(s-\dfrac{\omega_0+\omega_1}{2\pi}n\right)\right]\\ \\ &= 2 \cdot \left[\dfrac{1}{2\pi}\sum_{n=-\infty}^\infty(-1)^n\left(\dfrac{1}{n+\frac{1}{2}}-\dfrac{1}{n-\frac{1}{2}}\right)\delta\left(s-\dfrac{\omega_0-\omega_1}{2\pi}n\right)\right]\\ &\quad* \;\;\left[\dfrac{1}{2\pi}\sum_{n=-\infty}^\infty(-1)^n\left(\dfrac{1}{n+\frac{1}{2}}-\dfrac{1}{n-\frac{1}{2}}\right)\delta\left(s-\dfrac{\omega_0+\omega_1}{2\pi}n\right)\right]\\ \\ &= 2\cdot \dfrac{1}{(2\pi)^2} \int_{-\infty}^\infty \left[\sum_{n=-\infty}^\infty c_n\delta\left(s-\tau-\dfrac{\omega_0-\omega_1}{2\pi}n\right)\right]\left[\sum_{m=-\infty}^\infty c_m\delta\left(\tau-\dfrac{\omega_0+\omega_1}{2\pi}m\right)\right]d\tau\\ \\ &= 2\cdot \dfrac{1}{(2\pi)^2} \int_{-\infty}^\infty \sum_{n=-\infty}^\infty \sum_{m=-\infty}^\infty c_nc_m\delta\left(s-\tau-\dfrac{\omega_0-\omega_1}{2\pi}n\right)\delta\left(\tau-\dfrac{\omega_0+\omega_1}{2\pi}m\right)d\tau\\ \\ &= 2\cdot \dfrac{1}{(2\pi)^2} \sum_{n=-\infty}^\infty \sum_{m=-\infty}^\infty c_nc_m\delta\left(s-\dfrac{\omega_0+\omega_1}{2\pi}m-\dfrac{\omega_0-\omega_1}{2\pi}n\right)\\ \\ &= 2\cdot \dfrac{1}{(2\pi)^2} \sum_{n=-\infty}^\infty \sum_{m=-\infty}^\infty (-1)^{m+n}\dfrac{1}{\left(m^2-\frac{1}{4}\right)\left(n^2-\frac{1}{4}\right)}\delta\left(s-\left[\dfrac{\omega_0+\omega_1}{2\pi}m+\dfrac{\omega_0-\omega_1}{2\pi}n\right]\right)\\ \\ \end{align*}$$

Where

$$* \quad\text{denotes convolution}$$ $$\mathrm{III}(as) = \dfrac{1}{|a|}\sum_{n=-\infty}^\infty \delta\left(s-\dfrac{n}{a}\right)$$ $$\mathrm{sinc(s)} = \dfrac{\sin(\pi s)}{\pi s} $$ $$ \mathscr{F}\left\{x(t)\right\} = \int_{-\infty}^\infty x(t) e^{-i2\pi st} dt$$

So each of the two terms of that convolution are an infinite train of Dirac delta function spikes modulated by the sum of two sinc$()$ functions whose centers are equally spaced away from the origin in opposite directions.

Not terribly hard to understand conceptually: the end result is going to be an infinite number of Dirac delta functions at various spacings and strengths. The amplitudes and positions of the delta functions will be subject to beating of the sum and difference of the $\omega_{(0,1)}$.

Answer 2

I have implemented and validated @AndyWalls' result for discrete use:

Not yet performance optimized, and the optimal version is fairly quick, heaviest steps being 1) three $N$-sized FFTs, 2) one $N$-sized complex mult, and maybe 3) four spaced sinc samplings (but can pre-compute and reuse). The greater the f0 and f1, the more periodization periods (k - 1) are needed for sinc's decay.

Doesn't work for non-integer frequencies (need leaking sincs). Code, for now:

import numpy as np
from numpy.fft import fft, ifft
import matplotlib.pyplot as plt


def plot(x, title=None):
    plt.plot(x)
    if title is not None:
        plt.title(title, weight='bold', fontsize=18, loc='left')
        plt.gcf().set_size_inches(10, 7)
        plt.show()

def dtrain0(t, A):  # for reference
    def dirac(t):
        return (np.abs(t) < 1000*np.finfo(t.dtype).eps).astype('float64')

    n_terms = len(t)
    ns = np.arange(-n_terms//2, n_terms//2)
    return 1/abs(A) * np.sum([dirac(t - n/A) for n in ns], axis=0)

def dtrain1(t, A):
    M = len(t)
    o = np.zeros(M)
    A = abs(int(np.round(1/A)))
    scale = t.max() / len(t) * 2
    period = int(A / scale)
    o[:M//2 + 1][::period] = A
    o[M//2:-1][::-1][period - 2::period] = A
    return o

dtrain = (dtrain0, dtrain1)[1]

def _ssinc(pm, f):
    return np.sinc(2*np.pi / pm * f + .5) + np.sinc(2*np.pi / pm * f - .5)

def formula(w0, w1, f, f0, f1, k=0):
    p = w0 + w1
    m = w0 - w1

    ss0 = np.sinc(2*np.pi / m * f + .5) + np.sinc(2*np.pi / m * f - .5)
    ss1 = np.sinc(2*np.pi / p * f + .5) + np.sinc(2*np.pi / p * f - .5)
    o0 = dtrain(f, 2*np.pi / m) * ss0
    o1 = dtrain(f, 2*np.pi / p) * ss1
    # fold fourier <=> subsample time
    o0 = o0.reshape(2**k, -1).sum(axis=0)
    o1 = o1.reshape(2**k, -1).sum(axis=0)
    conv = ifft(fft(o0) * fft(o1)).real

    K0 = 2*np.pi / np.abs(m)
    K1 = 2*np.pi / np.abs(p)
    o = .5 * K0 * K1 * conv
    return o


fs, T, f0, f1 = 64, 2, 9, 23
k = 5  # number of periodization periods, minus 1
N = fs * T

w0, w1 = 2*np.pi * f0, 2*np.pi * f1
t = np.linspace(0, T, N, 0)
pos = np.arange(N//2 * 2**k + 1)
neg = -np.arange(1, N//2 * 2**k)[::-1]
f = np.hstack([pos, neg]) / T

x0 = np.cos(w0 * t)
x1 = np.cos(w1 * t)

o_exact = fft(np.abs(x0 + x1)) / N
assert np.allclose(o_exact.imag, 0)
o_exact = o_exact.real
o_formula = formula(w0, w1, f, f0, f1, k=k)

adiff = np.abs(o_exact - o_formula)
plot(o_exact)
plot(o_formula, title="MAE={:.2e} | f0, f1, fs, T, k = {}, {}, {}, {}, {}".format(
    adiff.mean(), f0, f1, fs, T, k))