# Numerical evaluation of Fourier transform of a scaling function

Given a set of filter taps $$\{h_n\}_{n=0}^{m-1}$$, define a scaling function $$\phi$$ by $$\phi(x) = \sqrt{2}\sum_{n} h_n \phi(2x-n).$$ In keeping with the notation from Daubechies "Ten Lectures on Wavelets", define

$$m_0(\xi) := \frac{1}{\sqrt{2}} \sum_{n=0}^{m-1} h_{n}\exp(-in\xi).$$

The Fourier transform of the scaling function is

$$\mathcal{F}[\phi](\xi) = \frac{1}{\sqrt{2\pi}} \prod_{j=1}^{\infty} m_0(\xi/2^{j})$$

I wish to evaluate $$\mathcal{F}[\phi]$$ for any value of $$\xi$$, for the compactly supported Daubechies wavelets. Using Horner's method for evaluation of $$m_0$$, I have computed this function for various values of $$\xi$$ and timed the evaluation. Given a "sensible" finite set of filter taps an iterating over $$j$$ until $$m_0(\xi/2^j)$$ is within an ULP for 1, I have successfully got the values computed correctly. However, this method is extremely slow; to wit, I have measured it to take 163ns to recovery float precision and 352ns to recover double precision-making it about the slowest special function evaluation I know.

Is there a way to speed this calculation up?

Merge request: https://github.com/boostorg/math/pull/921/files

• What are typical values of $m$? The timing information in isolation is not useful without knowing the parameters of the platform. The latency seems unusually high. What is the processor architecture? What is the clock frequency? Does the processor have floating-point hardware, FMA support, SIMD support? What programming language and toolchain are used to compile the code? Have you profiled the code? FFTs in particular are often limited by memory bandwidth. Commented Jan 1, 2023 at 0:47
• Typical values of m around around (say) 8-16; see Table 1 here: services.math.duke.edu/~ingrid/publications/CPAM_Orth_bas.pdf Commented Jan 1, 2023 at 17:03
• Processor is ARM, 2.4GHz clock, FMA and SIMD available, C++ compiled with clang. Commented Jan 1, 2023 at 17:04
• @njuffa: Recompiled with -O3 -ffast-math; now I get 300ns evaluation in float precision and 950ns evaluation in double. I imagine this can still be sped up . . . but it is no longer just totally insane. Commented Jan 1, 2023 at 17:34
• Accuracy also does not seem great using this naive method; looks like I'm losing 4-5 bits off the mantissa in float precision. Commented Jan 1, 2023 at 18:15

Yes, an acceleration is possible, a basic acceleration for all admissible scaling sequences, and even faster convergence the smoother the resulting wavelet scaling function is. This is based on knowing a related scaling function with closed formula for the function itself and also its Fourier transform, more specifically the B-splines. The recursion equation for the difference has a smaller contraction constant, down to about $$2^{-A}$$ if $$A$$ is the number of vanishing moments (db2 has $$2$$, db8 with sequence length $$16$$ has $$A=8$$).

### My notation conventions

For my convenience, first switch notation to something more adapted to the task. Having a separate factor $$\sqrt2$$ is nice for exploring orthogonality of wavelets. For studying its convergence it is more helpful to have a factor $$2$$, which is also the reason to switch from $$h$$ to $$m_0$$ in the reference.

Shift operator Define $$(\tau f)(x)=f(x-1)$$. Then the scaling equation becomes $$\phi(x)=2(a(\tau)\phi)(2x)$$.

Dilation operator To remove the argument from the equation, additionally define a dilation or doubling operator $$(Df)(x)=f(2x)$$. Then the scaling equation becomes a fixed-point equation in function space $$\phi=D(a(\tau)\phi)$$. In a mash-up of these notations also use $$\hat \tau(\xi)=\exp(-i\xi)$$ for the basis shift.

Combining both operators results in $$\tau D=D\tau^2$$, and consequently the basic identity behind most of the following approach $$(1-\tau)D=D(1+\tau)(1-\tau).$$

Fourier and Z transforms It makes sense to introduce the formal (Laurent) polynomial $$a(z)$$. Denote the Fourier factor of the convolution operator $$a(\tau)$$ as $$\hat a(\xi)=a(\exp(-i\xi))$$. In consequence $$\hat\phi(\xi)=\hat a(\xi/2)\hat \phi(\xi/2)$$.

### Why slow convergence

The observed problem for the product is that the sequence $$\xi/2^j$$ converges fast toward $$0$$ and thus $$\hat m_0(\xi/2^j)$$ is ever closer to $$1$$, making for a very slow convergence.

### Vanishing moments and factorization

We know that $$a(1)=1$$ from the integrability of the scaling function and $$a(-1)=0$$ as necessary condition for a continuous scaling function. This implies a factorization $$a(z)=\frac{1+z}2p(z)$$ with $$p(1)=1$$. Thus $$p(z)=1+(1-z)q(z)$$.

### Comparing two scaling functions

Add a second pair $$(\tilde a,\tilde\phi)$$ of scaling sequence and function to the situation. Then with above conventions their difference can be represented as $$\phi-\tilde\phi=(1-\tau)\eta$$, and inserting the refinement equations \begin{align} (1-\tau)\eta &= 2D(a(\tau)\phi-\tilde a(\tau)\tilde\phi)\\ &=D(1+\tau)\Bigl(p(\tau)(\phi-\tilde\phi)+(p(\tau)-\tilde p(\tau))\tilde\phi)\Bigr)\\ &=(1-\tau)D\Bigl(p(\tau)\eta+(q(\tau)-\tilde q(\tau))\tilde\phi\Bigr) \end{align} This gives a recursion equation, one could say affine linear refinement equation, for $$\eta$$, $$\eta = D\Bigl(p(\tau)\eta+v(\tau)\tilde\phi\Bigr)$$ with $$v=q-\tilde q$$.

In the Fourier picture this gives the fixed-point iteration $$\hat\eta(\xi)=\frac12\left(\hat p(\xi/2)\hat\eta(\xi/2)+\hat v(\xi/2)\hat{\tilde\phi}(\xi/2)\right)$$ to be inserted in $$\hat\phi(\xi)=\hat{\tilde\phi}(\xi)+(1-\hat\tau(\xi))\hat \eta(\xi)\\$$ Applying the Banach fixed-point theorem, this gives uniform convergence if $$\sup|\hat p(\xi)|<2$$. This is now non-empty in the space of $$p$$ sequences, as the only a-priori restriction is $$p(1)=1$$.

Pointwise convergence is without restriction, as $$p$$ is polynomial with $$\hat p(0)=p(1)=1$$, so the sequence of factors $$\hat p(\xi/2^n)/2$$ converges to $$1/2$$.

A condition for convergence in $$L^2(\Bbb R)$$ is $$\sup|\hat p(\xi)|<\sqrt 2$$

### B-splines

A class of known solutions to the refinement equation are the basis or B-splines $$\beta_k$$ with $$\beta_0={\bf 1}_{[0,1]}$$ the box function and recursively $$\beta_{k+1}=\beta_0*\beta_k$$. $$\beta_1(x)=\max(0,1-|x-1|)$$ is the hat function. The refinement sequence for $$\beta_k$$ is given via $$a(z)=\left(\frac{1+z}{2}\right)^{k+1}$$. The Fourier transform gets computed via $$\hat\beta_0(\xi)=\int_0^1e^{-ix\xi}\,dx=\frac{1-e^{-i\xi}}{i\xi}=e^{-i\xi/2}\frac{\sin(\xi/2)}{\xi} \\ \implies\hat\beta_k=e^{-i(k+1)\xi/2}\sin^{k+1}(\xi/2)/\xi^{k+1}$$

### Examples with minimal support

Let $$q$$ be constant and take $$\tilde\phi=\beta_1$$, so that $$\tilde p=\frac{1+z}2=1-\frac12(1-z)$$. Then $$v(z)=q_0+\frac12$$. $$|\hat p(\xi)|=\sqrt{(1+q_0)^2-2(1+q_0)q_0\cos(\xi)+q_0^2}$$. The condition for uniform convergence is $$|1+q_0|+|q_0|=\max(1,|1+2q_0|)<2$$, that is, $$q_0\in(-1.5,0.5)$$.

The translation into an iteration for the Fourier transform follows directly.

The convergence is visually rapid in all cases, there is no visual difference between the 8th and (L=)50th iterate. How close $$q_0$$ is to $$-0.5$$ only determines how rapidly the amplitude falls, which is a measure for the smoothness of the scaling function itself.

The computation of a function table for a plot of the graph of the scaling function can be organized similarly. As $$p$$ has degree 1, length 2, the iteration function $$\eta$$ has support $$[0,1]$$. The value $$\eta(x)$$ is used in the iteration formula for the values at $$x/2$$ and $$(x+1)/2$$. $$\pmatrix{\eta(x/2)\\\eta(x/1+1/2)} = \pmatrix{1+q_0\\-q_0}\eta(x)+v_0\pmatrix{x\\1-x} \\ \pmatrix{\phi(x)\\\phi(x+1)} = \pmatrix{x\\1-x}+\pmatrix{1\\-1}\eta(x)$$ This form shows that $$\eta$$ converges uniformly and thus to a continuous function for $$\max(|1+q_0|,|q_0|)<1$$, that is, for $$q_0\in[-1,0]$$. However, only for $$q_0=-0.5$$ is the limit Lipschitz continuous, moving away one can show Hölder continuity, with the index going to $$0$$ at the boundaries of the interval.

## start with the most basic subdivision
x,eta = np.array(([0,1],[0,0])
etas = [[x,eta]]

## iterate the eta recursion
for k in range(10):
eta = np.concatenate([(1+q)*eta+v*x, -q*eta+v*(1-x)])
x = np.concatenate([x/2,(1+x)/2])
etas.append([x,eta])

## construct and evaluate the scaling function
for k in [1,3,7]:
x,eta = etas[k]
phi = np.concatenate([x+eta, 1-x-eta])
x = np.concatenate([x,1+x])
plot(x,phi)
# etc