3
$\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
    Commented Feb 4, 2021 at 16:01
  • $\begingroup$ @TimWescott Do these URLs work? $\endgroup$ Commented Feb 4, 2021 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
    Commented Feb 4, 2021 at 16:12
  • $\begingroup$ @TimWescott Admittedly URLs over code aren't too visible; edited $\endgroup$ Commented Feb 4, 2021 at 16:15

2 Answers 2

3
$\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$
1
$\begingroup$

Periodization, or Dirac-combing (to coin a phrase), undoes the effect of domain-restriction and provides the link between a Fourier series and the corresponding Fourier transform. Adopting the convention $$1^x ≡ e^{2πix} = \cos(2πx) + i \sin(2πx)$$ (which I've termed the Ramanujan convention) and, defining forward and backward Fourier transforms by the following convention $$f(t) = ∫ \hat f(ν) 1^{+νt} dν ⇔ \hat f(ν) = ∫ f(t) 1^{-νt} dt,$$ periodization is done with the Dirac Comb $Ш(x)$, which satisfies the identities $$\sum_{k∈ℤ} δ(x - k) = Ш(x) = \sum_{n∈ℤ} 1^{nx},$$ and domain-restriction, or the restriction of a function to a domain $F$, is done with the characteristic function $χ_F(x)$ defined as 1 if $x∈F$ and 0 if $x∉F$, or more succinctly as $$χ_F(x) = ∫_{x'∈F} δ(x - x') dx'.$$

(Note that you can switch the signs to $k$ or $n$ on either side of the distributional identity, since the sum is over all integers $ℤ = \{⋯, -2, -1, 0, +1, +2, ⋯\}$.)

The Dirac Comb is its own Fourier transform, as can be seen fairly directly from the identity relations above and is ultimately what generates the discretization formulae we're familiar with. It yields - as both corollary and generalization - the Poisson summation formula $$\sum_{k∈ℤ} f(x - k) = \sum_{n∈ℤ} \hat f(n) 1^{nx}.$$

With our interest, in mind, of applying this to the time domain over equally-spaced sampled domains $T = \{ t₀ + k Δt: k ∈ ℤ \}$, with spacing $Δt > 0$, and noting the identity $δ(t - kΔt) Δt = δ(t/Δt - k)$, we can generalize the distributional identity to the following $$\sum_{t'∈T} δ(t - t') Δt = \sum_{k∈ℤ} δ(t - t₀ - kΔt) Δt = Ш\left(t - t₀ \over Δt\right) = \sum_{n∈ℤ} 1^{n{t - t₀ \over Δt}}.$$

The characteristic function is the kernel of domain restriction. We're interested to using it in the frequency domain $F = [ν₋,ν₊]$, with $ν_± = ν₀ ± ½Ω$ with finite bandwidth $Ω > 0$ and center frequency $ν₀$. So, we can write for its inverse Fourier transform, $$∫ χ_F(ν) 1^{+νt} dν = ∫_{ν∈F} 1^{+νt} dν = 1^{+ν_0t} Ω {\rm sinc}(Ωt),$$ with the cardinal example being that for the interval $F = [-½,+½]$: $${\rm sinc}(t) ≡ \int_{-½}^{+½} 1^{+νt} dν = {\sin πt \over πt}.$$

Domain-Restriction Versus Dirac-Comb
a.k.a. Fourier Series Versus Fourier Transform.

The question that arises is: given the following discrete Fourier transform relations $$f_F(t) = ∫_{ν∈F} \hat f(ν) 1^{+νt} dν ⇔ \hat f(ν) = \sum_{t∈T} f_F(t) 1^{-νt} Δt$$ where $Ω Δt ≤ 1$ for either tight sampling with $Ω Δt = 1$ or oversampling with a redundancy factor $1/(Ω Δt) > 1$, what relations do these have to the respective Fourier transforms $f(t)$ and $\widehat{f_F}(ν)$?

1: $\hat f(ν) ⇒ \widehat{f_F}(ν), f(t) ⇒ f_F(t)$ and Domain Restriction
First, given the previous discussion, we can immediately write $$\widehat{f_F}(ν) = χ_F(ν) \hat f(ν).$$ In the special case where the spectrum of $f$ is limited to $F$, i.e. $supp(\hat f) ⊆ F$, this gives us $\hat f = \widehat {f_F}$ (and, thus) $f = f_F$. Applying the inverse Fourier transform, this yields $$f_F(t) = ∫ f(t') 1^{ν₀(t - t')} Ω {\rm sinc}(Ω(t - t')) dt'.$$

2: $f_F(t) ⇒ f(t)$ and the generalized Shannon-Nyquist Sampling Theorem
Second, by directly substituting in the inverse Fourier transform of $\widehat{f_F}$ for the transform formula for $f(t)$, we find $$\begin{align} f(t) &= ∫ \left(\sum_{t'∈T} f_F(t') 1^{-νt'} Δt \right) 1^{νt} dν \\ &= \sum_{t'∈T} f_F(t') \left(∫ 1^{ν(t-t')} dν \right) Δt \\ &= \sum_{t'∈T} f_F(t') δ(t-t') Δt \\ &= f_F(t) \sum_{t'∈T} δ(t-t') Δt \\ &= f_F(t) Ш\left({t - t_0} \over Δt \right) \end{align}$$

Upon back-substitution, this yields the generalization of the Shannon-Nyquist sampling formula for general intervals $F$ and $T$: $$\begin{align} f_F(t) &= ∫ \left(\sum_{t'∈T} f_F(t') δ(t"-t') Δt\right) 1^{ν₀(t - t")} Ω {\rm sinc}(Ω(t - t")) dt" \\ &= \sum_{t'∈T} ∫ f_F(t') δ(t"-t') 1^{ν₀(t - t")} Ω {\rm sinc}(Ω(t - t")) dt" Δt \\ &= \sum_{t'∈T} f_F(t') 1^{ν₀(t - t')} Ω {\rm sinc}(Ω(t - t')) Δt. \end{align}$$

3: $\widehat{f_F}(ν) ⇒ \hat f(ν)$ and Periodization
Third, and last, since $Ш$ is its own forward and inverse Fourier transform, then it follows from the $f_F(t) ⇒ f(t)$ formula that $$\begin{align} \hat f(ν) &= ∫ \widehat{f_F}(ν') 1^{-(ν-ν')t_0} Ш((ν - ν')Δt) Δt dν' \\ &= ∫ \widehat{f_F}(ν') 1^{-(ν-ν')t_0} \sum_{n∈ℤ} δ\left(ν - ν' + {n \over Δt}\right) dν' \\ &= \sum_{n∈ℤ} ∫ \widehat{f_F}(ν') 1^{-(ν-ν')t_0} δ\left(ν - ν' + {n \over Δt}\right) dν' \\ &= \sum_{n∈ℤ} \widehat{f_F}\left(ν + {n \over Δt}\right) 1^{{n \over Δt}t_0} \end{align}$$

Note from the forward discrete transform formula for $\hat f$, that $\hat f$ is periodic up to phase, with $\hat f(ν + 1/Δt) = 1^{-t_0/Δt} \hat f(ν)$; i.e. $\hat f(ν) 1^{νt_0}$, itself, has period $1/Δt$. So, this is the generalization of periodization in the frequency domain that goes with a non-zero offset $t_0$ in the time domain.

This last formula, therefore, says that $\hat f(ν)$ is obtained by running $ν$ up or down by multiples of $1/Δt$, with appropriate adjustments in phase, until you get to whatever $ν'$ is as close to the domain $F$ as possible, or within it, and the result is $\widehat {f_F}(ν')$ if $ν'∈F$ or 0 otherwise.

Apparently, this can all be applied to undersampling as well, where $Ω Δt > 1$, except the sum in the $\widehat{f_F}(ν) ⇒ \hat f(ν)$ formula can yield two or more terms, rather than just 1 or 0. For tight sampling there is only one $ν'$ that lands in the interval $F$ ... except if it hits the boundary points $ν_±$. So, it's possible to fix that by adopting a convention, like $χ_F(ν_±) = ½$; and I think the technical formulation of the relevant theorems makes it mandatory.

$\endgroup$

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.