# What is “filter periodization”?

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]

• 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. – TimWescott Feb 4 at 16:01
• @TimWescott Do these URLs work? – OverLordGoldDragon Feb 4 at 16:06
• You had the links in there and I didn't notice? Oh, yes! I'm paid to think! – TimWescott Feb 4 at 16:12
• @TimWescott Admittedly URLs over code aren't too visible; edited – OverLordGoldDragon Feb 4 at 16:15

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