1
$\begingroup$

A library defines periodize_filter_fourier, which is an equi-spaced averaging formulated by

$$ v_f[k] = \sum_{i=0}^{\text{n_periods}-1} h_f[i\cdot N + k], $$

where $v_f$ is periodization of $h_f$, $N=\text{len}(h_f)/\text{n_periods}=\text{len}(v_f)$, $k=[0,1,...,N-1]$, and $1/\text{n_periods}$ seems missing.

It's applied on the frequency-domain (DFT) Morlet wavelet. What's the purpose, and what's "periodizing" to do with it?


Example with h = [1, 2, 3, ..., 12]:

# n_periods == 2
v_f = [[1, 2, 3, 4, 5, 6],
       [7, 8, 9, 10, 11, 12]].mean(axis=0) 
    = [4, 5, 6, 7, 8, 9]

# n_periods == 3
v_f = [[1, 2, 3, 4],
       [5, 6, 7, 8],
       [9, 10, 11, 12]].mean(axis=0)
    = [5, 6, 7, 8]
$\endgroup$
4
  • $\begingroup$ Since this is edging toward the esoteric, please tell us what library this is, and if possible give us a link to the code or documentation on the library. I can think of reasons why you may want to do this, but wouldn't know if I'm right until I work through the math. $\endgroup$ – TimWescott Feb 4 at 16:01
  • $\begingroup$ @TimWescott Do these URLs work? $\endgroup$ – OverLordGoldDragon Feb 4 at 16:06
  • 1
    $\begingroup$ You had the links in there and I didn't notice? Oh, yes! I'm paid to think! $\endgroup$ – TimWescott Feb 4 at 16:12
  • $\begingroup$ @TimWescott Admittedly URLs over code aren't too visible; edited $\endgroup$ – OverLordGoldDragon Feb 4 at 16:15
1
$\begingroup$

TL;DR The idea is to modify the continuous-time representation of the filter, $f(t)$, such that it's suitable for discrete-time sampling - mainly by accounting for aliasing via "folding" in Fourier domain. The operation is also equivalent to (and used for) subsampling.

"Periodization" refers to the fact that sampling $f(t)$ at uniform intervals $s$ corresponds to making its Fourier transform $2\pi/s$-periodic.


Periodizing

A uniform sampling of $f$ with interval $s$ can be expressed as a weighted Dirac sum:

$$ f_d(t) = \sum_{n=-\infty}^{\infty} f(ns) \delta (t - ns) \tag{1} $$

Fourier transform of $\delta(t - ns)$ is $e^{-jns\omega }$, so the Fourier transform of $f_d$ is

$$ \hat f_d (\omega) = \sum_{n=-\infty}^{\infty} f(ns) e^{-ins\omega} \tag{2} $$

($f(ns)$ is the "weights" and thus constant in transform over $t$). It can be shown${}^1$ that this is same as

$$ \hat f_d (\omega) = \frac{1}{s}\sum_{k=-\infty}^{\infty} \hat f \left( \omega - \frac{2\pi}{s}k \right) \tag{3} $$

That is, sampling $f$ at intervals $s$ is equivalent to making its Fourier transform $2\pi/s$-periodic by summing all its translations $\hat f (\omega - (2\pi / s)k)$. This makes more sense with what follows.


Sampling

(Shannon-Whittaker sampling theorem) If the support of $\hat f$ is bounded in $[-\pi/s, \pi/s]$, then

$$ f(t) = \sum_{n=-\infty}^{\infty} f(ns) \phi_s(t - ns), \tag {4} $$

where

$$ \phi_s(t) = \frac{\sin{( (\pi/s)t )}}{(\pi / s) t}. \tag{5} $$

Note, if $n \neq 0$, the support of $\hat f(\omega - (\pi / s)n)$ does not intersect the support of $\hat f (\omega)$ because $\hat f(\omega) = 0 $ for $|\omega| > \pi / s$. So $(3)$ implies

$$ \hat f_d (\omega) = \frac{\hat f(\omega)}{s}\ \ \text{if}\ \ |\omega| \leq \frac{\pi}{s} \tag{6} $$

Fourier transform of $\phi_s$ is $\hat \phi_s = s \bf{1}_{[-\pi/s, \pi/s]}$ (indicator function, or "box"). Since the support of $\hat f$ is in $[-\pi/s, \pi/s]$, it follows from $(6)$ that $\hat f(\omega) = \hat \phi_s(\omega) \hat f_d(\omega)$. Its inverse Fourier transform gives

$$ \begin{align} f(t) = \phi_s \star f_d(t) &= \phi \star \sum_{n=-\infty}^{\infty} f(ns) \delta(t - ns) \\ &= \sum_{n=-\infty}^{\infty} f(ns) \phi_s(t - ns) \tag{7} \end{align} $$

This happens to be a proof of the theorem.


Sampling $\Leftrightarrow$ Periodizing

This can be seen as follows; suppose $\hat f$'s support is in $[-\pi/s, \pi/s]$. Then:

$(a)$ is simply $f(t)$ and its F.T. Sampling it at uniform interval $s$ has the effect of adding translations of its F.T. (i.e. $(3)$), which is what's shown - original $\hat f(\omega)$ translated and centered at multiples of $(2\pi /s)$. These "duplicates" extend out to infinity; we've thus "periodized" the Fourier transform of $f(t)$.

Per $(7)$, we recover $f(t)$ from $f_d(t)$ by convolving with $\phi_s$ (or multiplying in Freq domain):


Subsampling

Let's pose it as a problem: what effect does subsampling (i.e. taking every $k$-th sample) of a finite sequence $x[n] = [4, 7, ..., -1]$ have on its Fourier transform? Note this is equivalent to scaling the sampling interval: $s \rightarrow ks$. Suppose we wish to subsample by 2.

To answer this, we first visit aliasing.

Aliasing

In case $\hat f$ isn't bound to $[-\pi/s, \pi/s]$, the resulting $f_d(t)$ will alias. From $(3)$ and proof of $(4)$, we can predict this manifests as overlaps in Fourier domain:

The translations are still spaced by $(2\pi / s)$, while support of $\hat f$ exceeds it, yielding overlaps. The sampled $f_d(t)$ will then no longer resemble $f(t)$:

Subsampling (cont'd)

If $k=2$, we double our $s$, and thus halve the support $[-\pi/s, \pi/s]$. To subsample without aliasing, $\hat f$'s support must thus lie in $[-\pi/(2s), \pi/(2s)]$ (per original $s$):

Now suppose it isn't, and the support lies exactly in $[-\pi/s, \pi/s]$. Then:

What is the resulting $\hat f_d$? Well, it's simply the negative half of $\hat f(\omega)$ added to its positive half, and vice versa ("folding") - then divided by 2 (s *= 2 in eq $(3)$).

This is exactly what

$$ v_f[k] = \frac{1}{k} \sum_{i=0}^{\text{n_periods}-1} h_f[i\cdot N + k], \tag{8} $$

(expression in question) is doing, but instead with DFT bins, where the de-facto support w.r.t. continuous-time input is $[-\pi/s, \pi/s]$ (with or without aliasing).

Put one way, subsampling simulates the effects of aliasing. In code, the following holds (assume fold() implements $(8)$):

ifft(fold(fft(x), k)) == x[::k]

Periodizing a filter

Does $(8)$ serve any purpose beyond subsampling? Correctness: if we disregard folding and pretend aliasing isn't there, result's even worse than with aliasing.

More generally, given some arbitrary continuous-time function of a filter in Fourier domain, $\hat \psi(\omega)$, how do we sample it such that the entire support (assuming finite) is included, and that we use exactly $N$ samples? One approach is, "sample lots, then subsample" - which is exactly what Kymatio's morlet_1d does. (The alternative is to determine correct sampling bounds and get it right the first time.)

Once we have the correct $\hat \psi_d(\omega)$ that accounts for aliasing and captures at least one full period of the continuous-time function (i.e. $[-\pi/s, \pi/s]$), then $\hat \psi_d(\omega)$ is the correct continuous-time function, whether without or with aliasing, and taking IFFT we get closest to what we'd get with sampling the continuous-time $\psi (t)$ with rate $f_s=1/s$ (over one continuous-time period).

References

Figures and formulae (with much accompanying wording copied verbatim) taken from S. Mallat's Wavelet Tour, sections 3.1.0-3.1.2, where further explanations and proofs are given (including to footnote 1).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.