What is wrong with my Meyer wavelet implementation?

Question

In general, given a mother wavelet $\psi \in L^2(\mathbb{R})$, a wavelet frame is constructed by $$ \psi_{j,k} (x) := 2^{j/2} \psi(2^jx - k), $$ for $j,k\in \mathbb{Z}$. Given such a wavelet frame, every element $x \in L^2(\mathbb{R})$ can be written as $$ \sum_{j,k\in \mathbb{Z}} \left\langle x, \psi_{j,k} \right\rangle \psi_{j,k}.$$ With $\Psi$ denoting the synthesis and $\Psi^*$ denoting the analysis operator, this an be written as $x = \Psi \Psi^* x$.

This also holds true in finite dimensions, where $x$ is a vector and $\Psi$ is a matrix. I want to implement a method that, given a wavelet function, returns the matrix $\Psi$, such that I can calculate the identity $x = \Psi \Psi^* x$, where $\Psi^*$ is the complex conjugate, i.e. the transpose if $\Psi$ is real.

I can easily implement this if $\psi$ is the Haar wavelet (see code below and left image). However, I can't make the same approach work for the Meyer wavelet (right image):

Can anyone give some insight on what might go wrong with the Meyer wavelet approach?

Thanks in advance!

import numpy as np
import matplotlib.pyplot as plt

def haar_wavelet(t):
    return np.where((t >= 0) & (t < 0.5), 1, np.where((t >= 0.5) & (t <= 1), -1, 0))

def haar_scaling(t):
    psi = np.zeros_like(t)
    psi[(0 <= abs(t)) & (abs(t) <= 1)] = 1
    return psi

def meyer_scaling(x):
    psi = (np.sin(2*np.pi/3*x) + 4/3*x*np.cos(4*np.pi/3*x))/(np.pi*x-16*np.pi/9*(x**3))
    psi[x==0] = 2/3 + 4/(3*np.pi)
    return psi

def meyer_wavelet(x):
    psi1 = (4/(3*np.pi)*(x-0.5)*np.cos(2*np.pi/3*(x-0.5))-1/np.pi*np.sin(4*np.pi/3*(x-0.5)))/((x-0.5)-16/9*(x-0.5)**3)
    psi2 = (8/(3*np.pi)*(x-0.5)*np.cos(8*np.pi/3*(x-0.5))+1/np.pi*np.sin(4*np.pi/3*(x-0.5)))/((x-0.5)-64/9*(x-0.5)**3)
    psi = psi1 + psi2
    return psi

N = 512
x = np.linspace(0, 1, N)
dx = x[1]-x[0]

# Max. number of scales
J = int(np.log2(N))

# Get the Haar scaling function
W_haar = [haar_scaling(x)[:,None]]
for j in range(0, J):
    # Number of shifts per scale
    ns = int(2**j) 
    for k in range(0, ns):
        # Get all translations k at scale j
        haar_jk = 2**(j/2)*haar_wavelet(2**j*x - k)
        W_haar.append(haar_jk[:,None])

# Synthesis operator for Haar wavelet
W_haar = np.concatenate(W_haar, axis=1) 

# Get the Meyer scaling function
W_meyer = [meyer_scaling(x)[:,None]]
for j in range(0, J):
    # Number of shifts per scale
    ns = int(2**j)
    for k in range(0, ns):
        # Get all translations k at scale j
        meyer_jk = 2**(j/2)*meyer_wavelet(2**j*x - k)
        W_meyer.append(meyer_jk[:,None])

# Synthesis operator for Meyer wavelet
W_meyer = np.concatenate(W_meyer, axis=1)

# Create example signal 
signal = np.zeros(n)
signal[N//2-N//4:N//2+N//4] = 1

# Calculate wavelet coefficients
coeff_haar = (W_haar.T @ signal) * dx
coeff_meyer = (W_meyer.T @ signal) * dx

# Calculate reconstruction
rec_haar = W_haar @ coeff_haar
rec_meyer = W_meyer @ coeff_meyer

# Plot The reconstructions
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
ax1.plot(x, signal, label="orignal")
ax1.plot(x, rec_haar, '--', label="reconstructed")
ax1.legend(loc="upper right")
ax2.plot(x, signal, label="original signal")
ax2.plot(x, rec_meyer, '--', label="reconstructed")
ax2.legend(loc="upper right")

Answer 1

There are problems in both the Haar and Meyer wavelet implementations.

Firstly, the Haar implementation. If you zoom in, you can see a difference between the signal and the reconstruction (it becomes quite obvious when N is small). This is caused by incorrectly normalising. One way to correct this is to replace your computation of coeff_haar with:

coeff_haar = (W_haar.T @ signal) / N

Secondly, the Meyer implementation. The difference between the signal and the reconstruction is significant. This is due to your Meyer matrix not forming an orthogonal basis set, nor does it even have full rank. The Meyer wavelet is a continuous wavelet, yet you are trying to treat it as if it is a discrete wavelet. If you want to pursue it, I found that Matlab has a discrete approximation of the Meyer wavelet, which references the book cited below, and that FIR approximations have been published.

P. Abry, Ondelettes et turbulence, Diderot ed., Paris, 1997, p. 268.