- Background
- Approach
- Symbol Definitions
- Core Equations
- Visualization
- Daubechies Coefficients
- Calculation of Scaling Function Initial Values
- Calculation of Scaling Function
- Calculation of Wavelet Function
- Source Code
- References
Background
Daubechies Wavelets are a family of orthogonal wavelets which are recursively-defined. They are discrete functions where each level of approximation fills in midpoints between points across calculated from the previous level and, in the infinite limit, becomes continuous. There are two equivalent nomenclatures: DX and dbX. The db4 (D8) wavelet contains 4 vanishing moments and 8 taps (coefficients) over the non-zero interval
Approach
A floating point approach is used for simplicity over the dyadic approach which uses fractional powers of two,
This article aims only to demonstrate the calculation of the wavelet and scaling function which are not well-explained in any of the many sources Iβve read. No theory, mathematical proofs, or description of the related wavelet transform will be given as those are covered elsewhere. The db4 wavelet will serve as the primary example here because the db2 wavelet is the dominant example in the literature. Values for wavelets other than db4 will occasionally be given to serve as an aid to general implementation.
Symbol Definitions
π(r) - Scaling function at point r
π(r) - Wavelet function at point r
N - number of Daubechies cofficients for the given wavelet
k - coefficient index
Core Equations
Initial Values for Scaling Function
Scaling Function
Wavelet Function
Visualization
The images below demonstrate db2 and db4 wavelets being recursively filled from 0 to 8 levels of approximation.
db2 wavelet and scaling functions
db4 wavelet and scaling functions
Number of points at various levels of approximation
0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | |
---|---|---|---|---|---|---|---|---|---|
db2 | 4 | 7 | 13 | 25 | 49 | 97 | 193 | 385 | 769 |
db3 | 6 | 11 | 21 | 41 | 81 | 161 | 321 | 641 | 1281 |
db4 | 8 | 15 | 29 | 57 | 113 | 225 | 449 | 897 | 1793 |
db5 | 10 | 19 | 37 | 73 | 145 | 289 | 577 | 1153 | 2305 |
db6 | 12 | 23 | 45 | 89 | 177 | 353 | 705 | 1409 | 2817 |
Evaluating via C#
int PointsAtLevel(int zeroLevelPoints, int level) =>
zeroLevelPoints * (1 << level) - ((1 << level) - 1);
Daubechies Coefficients
Daubechies coefficients for the scaling function are frequently given in texts and are given here as a table. The calculation of the scaling coefficients may be covered in a future blog. Briefly, the calculation involves a matrix solver with multiple constraints (example) and yields multiple solutions. The final constraint is the convention that the chosen solution is the one which has the largest values front-loaded. This is known as the extremal phase. The sum of coefficients for a given wavelet index is often normalized to 1 and this is the case here.
Scaling function coefficients
db2 | 0.6830127 | 1.1830127 | 0.3169873 | -0.1830127 | ||||
db3 | 0.47046721 | 1.14111692 | 0.650365 | -0.19093442 | -0.12083221 | 0.0498175 | ||
db4 | 0.32580343 | 1.01094572 | 0.89220014 | -0.03957503 | -0.26450717 | 0.0436163 | 0.0465036 | -0.01498699 |
Wavelet function coefficients
Wavelet function coefficients are directly derived from the scaling function coefficients via the relationship
The db4 wavelet coefficients calculations in detail:
db2 | -0.1830127 | -0.3169873 | 1.1830127 | -0.6830127 | ||||
db3 | 0.0498175 | 0.12083221 | -0.19093442 | -0.650365 | 1.14111692 | -0.4704672 | ||
db4 | -0.01498699 | -0.0465036 | 0.0436163 | 0.26450717 | -0.03957503 | -0.8922001 | 1.01094568 | -0.325803429 |
IEnumerable<float> ToWaveletCoefficients(IEnumerable<float> scalingCoefficients) =>
scalingCoefficients.Select((coef, index) => index % 2 == 0 ? -coef : coef).Reverse();
Calculation of Scaling Function Initial Values
The scaling function requires values at integer values of r within the waveletβs range of [0, N-1] with each end of the interval being zero (
And yields the following system of equations for db4:
This system can be greatly simplified, given that many values are zero. The endpoints and all other zero values will be removed:
Next, a matrix equation can be written of the form
Solving for Initial Values via Python
This matrix equation is solved with the constraint that
# Helper method for minimization: ||Ax||^2
def f(x):
y = np.dot(A, x)
return np.dot(y, y)
# maps 0 or coefficient values to matrix A
def map_index(a, i, j):
index = 2 * i - j
if index < 0 or index >= len(a):
return 0
else:
return a[index]
# db4 scaling coefficients
aVector = [ 0.32580343, 1.01094572, 0.89220014, -0.03957503, -0.26450717, 0.0436163, 0.0465036, -0.01498699 ]
N = len(aVector)
A_list = [[map_index(aVector, row, col) for col in range(1, N-1)] for row in range(1, N-1)]
A = np.array(A_list) - np.identity(N - 2)
print(f'A: {A}\n')
zeroVector = np.array(np.zeros(N - 2))
cons = ({'type': 'eq', 'fun': lambda x: x.sum() - 1})
res = optimize.minimize(f, zeroVector, method='SLSQP', constraints=cons, options={'disp' : False})
xbest = res['x']
print(f'db4 initial values: {xbest}\n')
print(f'sum of initial values: {xbest.sum()}')
Table of Initial Values
0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | |
---|---|---|---|---|---|---|---|---|
db2 | 0 | 1.36602544 | -0.36602544 | 0 | ||||
db3 | 0 | 1.28653824 | -0.38600335 | 0.09525722 | 0.00420788 | 0 | ||
db4 | 0 | 1.00716450 | -3.38492330e-2 | 3.96088947e-2 | -1.17337313e-2 | -1.25087227e-3 | 6.04419105e-5 | 0 |
Calculation of Scaling Function
The scaling function, π(r), is recursively-defined and requires the initial values from the previous step as well as the Daubechies coefficients. As mentioned in the intro, the function is not continuous until an infinite level of approximations are calculated, so we must settle for a discrete solution.
The calculation begins by initializing the initial level,
π(0) | π(1) | π(2) | π(3) | π(4) | π(5) | π(6) | π(7) |
---|---|---|---|---|---|---|---|
0 | 1.00716450 | -3.38492330e-2 | 3.96088947e-2 | -1.17337313e-2 | -1.25087227e-3 | 6.04419105e-5 | 0 |
Next, we calculate the values for
Simplify with the knowledge that for db4,
Note that the right-hand side contains known values of π(r) from the
Then simplify:
Once again, all of the values of π(r) required to evaluate the right-hand side were calculated during the previous step. The solutions for
π(0) | 0 | π(2.5) | -0.2420 | π(5) | -1.251e-3 | ||
π(0.25) | 0.1069 | π(2.75) | -0.1753 | π(5.25) | -4.077e-4 | ||
π(0.5) | 0.3281 | π(3) | 0.0396 | π(5.5) | 1.203e-4 | ||
π(0.75) | 0.6175 | π(3.25) | 0.1182 | π(5.75) | -2.691e-5 | ||
π(1) | 1.0072 | π(3.5) | 0.0343 | π(6) | 6.044e-5 | ||
π(1.25) | 1.1008 | π(3.75) | 0.0163 | π(6.25) | -1.845e-6 | ||
π(1.5) | 0.8773 | π(4) | -0.0117 | π(6.5) | -9.058e-7 | ||
π(1.75) | 0.5363 | π(4.25) | -0.0235 | π(6.75) | 1.358e-8 | ||
π(2) | -0.0338 | π(4.5) | 2.166e-3 | π(7) | 0 | ||
π(2.25) | -0.3020 | π(4.75) | 5.284e-3 |
Evaluating the Scaling Function via C#
An iterative approach is used here rather than a true recursive method. The following method calculates the scaling functionβs next level of approximation and is valid for
Vector2[] GetScalingLevel(int level, Vector2[] previous, float[] scalingCoefs)
{
var scaling = new Vector2[PointsAtLevel(scalingCoefs.Length, level)];
float maxRange = scalingCoefs.Length - 1f;
float distancePerIndex = maxRange / (scaling.Length - 1);
for (int n = 0; n < scaling.Length; n += 2)
scaling[n] = previous[n / 2];
for (int n = 1; n < scaling.Length; n += 2)
{
float sum = 0;
float r = n / (float)(1 << level);
for (int k = 0; k < scalingCoefs.Length; k++)
{
float rPrevious = 2 * r - k;
int rIndex = (int)MathF.Round(rPrevious / distancePerIndex);
if (rPrevious >= 0 && rPrevious <= maxRange)
sum += scalingCoefs[k] * scaling[rIndex].Y;
}
scaling[n] = new Vector2(r, sum);
}
return scaling;
}
Calculation of Wavelet Function
The wavelet function, π(r), is more straightforward once the scaling function has been calculated. Its equation is given by:
Given that
β¦
And simplifies to:
β¦
π(0) | 0 | π(2.5) | -0.0465 | π(5) | -0.0237 | ||
π(0.25) | -4.918e-3 | π(2.75) | -0.3901 | π(5.25) | -9.0904e-3 | ||
π(0.5) | -0.0151 | π(3) | -0.8872 | π(5.5) | 2.5044e-3 | ||
π(0.75) | -0.0284 | π(3.25) | -0.4322 | π(5.75) | 5.8323e-4 | ||
π(1) | -0.0463 | π(3.5) | 1.0437 | π(6) | 4.6864e-4 | ||
π(1.25) | -0.0229 | π(3.75) | 0.9951 | π(6.25) | -4.0116e-5 | ||
π(1.5) | 0.0449 | π(4) | -0.3976 | π(6.5) | -1.9692e-5 | ||
π(1.75) | 0.1358 | π(4.25) | -0.5611 | π(6.75) | -2.9513e-7 | ||
π(2) | 0.2633 | π(4.5) | 0.0616 | π(7) | 0 | ||
π(2.25) | 0.2069 | π(4.75) | 0.1116 |
Evaluating the Wavelet Function via C#
The wavelet function code is very similar to the scaling function code. The level is baked into scaling and isnβt required.
Vector2[] GetWaveletLevel(Vector2[] scaling, float[] waveletCoefs)
{
var wavelet = new Vector2[scaling.Length];
float maxRange = waveletCoefs.Length - 1f;
float distancePerPoint = maxRange / (wavelet.Length - 1);
for (int n = 0; n < wavelet.Length; n++)
{
float sum = 0;
float r = n * distancePerPoint;
for (int k = 0; k < waveletCoefs.Length; k++)
{
float rPrevious = 2 * r - k;
int rIndex = (int)MathF.Round(rPrevious / distancePerPoint);
if (rPrevious >= 0 && rPrevious <= maxRange)
sum += waveletCoefs[k] * scaling[rIndex].Y;
}
wavelet[n] = new Vector2(r, sum);
}
return wavelet;
}
Source Code
Main Project
Scaling and Wavelet Function (C#))
Initial Scaling Values (Python)
References
Wavelets Made Easy - Yves Nievergelt (2001)
Wavelet Browser - PyWavelets