Wavelet CDF 9/7 Implementation
The MATLAB function waveletcdf97.m included in this package is a self-contained M-function for applying the Cohen–Daubechies–Feauveau 9/7 (CDF 9/7) wavelet transform. This wavelet is an especially effective biorthogonal wavelet, used by the FBI for fingerprint compression and selected for the JPEG2000 standard .
Y = waveletcdf97(X, L) decomposes
L stages of the CDF 9/7 wavelet. For the inverse transform,
waveletcdf97(X, -L) inverts
L stages. Filter boundary handling is half-sample symmetric.
X may be of any size; it need not have size divisible by
2L. For example, if
X has length 9, one stage of decomposition produces a lowpass subband of length 5 and a highpass subband of length 4. Transforms of any length have perfect reconstruction (exact inversion).
X is a matrix,
waveletcdf97 performs a (tensor) 2D wavelet transform. If
X has three dimensions, the 2D transform is applied along the first two dimensions.
This package includes a demo program,
waveletcdf97_demo, that uses
waveletcdf97 for signal approximation. Signal approximation is the problem of representing a signal with as few components as possible. This is similar to lossy image compression, but ignoring the problems of quantization and encoding.
Wavelets are particularly efficient for approximating piecewise-smooth signals. This first demo uses the CDF 9/7 wavelet to represent a piecewise-smooth signal.
NumComponents = 40; % Approximation with 40 components
%%% Construct the test signal %%%
N = 512; % Signal length
t = linspace(-1.7,1.7,N);
X = sign(t).*exp(-t.^4);
%%% Wavelet approximation %%%
Level = 9; % Use 9 levels of decomposition
Y = waveletcdf97(X,Level); % Transform the signal
Y = keep(Y,NumComponents); % Keep only 40 components
R = waveletcdf97(Y,-Level); % Invert to obtain the approximation
norm(X-R) % Compute error
%%% Fourier approximation %%%
Y = fft(X); % Transform
Y = keep(Y,NumComponents); % Keep only 40 components
R2 = real(ifft(Y)); % Invert
norm(X-R2) % Compute error
The figure below shows the resulting wavelet approximation using only 40 out of 512 components. Fourier approximation with 40 components is shown for comparison. The wavelet approximation has L2 error of 0.014 while the Fourier approximation has error of 2.244.
The second demo applies
waveletcdf97 to image approximation. First, the input image is converted from RGB to the JPEG Y'CbCr colorspace. The Y'CbCr image is transformed using
waveletcdf97, all but the largest transform coefficients are set to zero, and then inverse transformed.
X = double(imread('palm.jpg'))/255; % Load the demo image
X = RGBToYCbCr(X); % Convert to Y'CbCr
L = 6;
Y = waveletcdf97(X,L); % Transform the image
R = waveletcdf97(keep(Y,1/40),-L); % 40:1 approximation
R = waveletcdf97(keep(Y,1/80),-L); % 80:1 approximation
In the first test, a random signal X is transformed one stage, then inverse transformed, and the result is compared to the original. Mathematically, the transform is exactly inverted—the scheme is said to have perfect reconstruction. This test verifies that this property holds to machine precision.
Runs = 100; % Number of runs to average
N = ceil(logspace(... % Lengths to test
for k = 1:length(N)
% Create random input matrices
X = rand(N(k),1,Runs);
% Forward transform followed by inverse
R = waveletcdf97(waveletcdf97(X,1),-1);
% Compute the average error
AvgError(k) = mean(max(abs(permute(X - R,[1,3,2])),,1));
fprintf('%3d: Error = %.2e\n',N(k),AvgError(k));
15: Error = 3.34e-016
19: Error = 3.45e-016
24: Error = 4.25e-016
29: Error = 4.30e-016
36: Error = 4.63e-016
44: Error = 4.91e-016
55: Error = 5.00e-016
68: Error = 5.53e-016
84: Error = 5.55e-016
103: Error = 5.99e-016
128: Error = 5.90e-016
158: Error = 6.58e-016
196: Error = 6.90e-016
243: Error = 7.17e-016
300: Error = 7.00e-016
The CDF 9/7 wavelet is designed such that where the input signal is locally a polynomial of cubic degree or lower, the resulting detail (highpass) coefficients are equal to zero. A wavelet is said to have "N vanishing moments" if it has this property on polynomials up to degree N-1, so CDF 9/7 has 4 vanishing moments. This test transforms a piecewise polynomial signal and displays the largest detail coefficient magnitudes, verifying that the vanishing moments hold to reasonable accuracy.
N = 64;
t = (0:N-1)/N;
X = [t.^0,t.^1,t.^2,t.^3];
Y = waveletcdf97(X,1);
norm(Y(2*N+2:2.5*N-2),inf) % Largest detail coefficient from t^0
norm(Y(2.5*N+2:3*N-2),inf) % Largest detail coefficient from t^1
norm(Y(3*N+2:3.5*N-2),inf) % Largest detail coefficient from t^2
norm(Y(3.5*N+2:4*N-2),inf) % Largest detail coefficient from t^3
Locally constant Largest detail coefficient = 1.41e-012
Locally linear Largest detail coefficient = 1.30e-012
Locally quadratic Largest detail coefficient = 1.20e-012
Locally cubic Largest detail coefficient = 1.11e-012
The code is written specialized for the CDF 9/7 wavelet in lifting scheme implementation. Sweldens' lifting scheme  represents a wavelet transform as a sequence of predict and update steps. Let X=[X(1), X(2), ... X(2N)] be an array of length 2N. The lifting scheme begins with the "polyphase decomposition," splitting X into two subbands, each of length N:
Xo = [X(1), X(3), X(5), ... X(2N-1)], Xe = [X(2), X(4), X(6), ... X(2N)].
Since Xo and Xe can be merged to recover X, no information has been lost.
Next, the scheme performs lifting steps on the subbands Xo and Xe. Let p be a filter, then
Xe'= Xe + p∗Xo
is called a prediction step, where ∗ denotes convolution. Similarly, Xo' = Xo + u∗Xe is called an update step. Notice that Xe can always be recovered from Xe' with
Xe'= Xe - p∗Xo.
This simple relationship between a forward step and an inverse step is the key to the lifting scheme: any sequence of prediction and update steps can be undone to recover Xo and Xe.
Any FIR wavelet transform can be factored into a sequence of lifting steps . For the CDF 9/7 wavelet, the lifting scheme decomposition used in waveletcdf97 is
Xo = [X(1),X(3),X(5), ... X(2N-1)]
Xe = [X(2),X(4),X(6), ... X(2N)]
Xe1(n) = Xe(n) + α(Xo(n+1) + Xo(n))
Xo1(n) = Xo(n) + β(Xe1(n) + Xe1(n-1))
Xe2(n) = Xe1(n) + δ(Xo1(n+1) + Xo1(n))
Xo2(n) = Xo2(n) + γ(Xe2(n) + Xe2(n-1))
The subbands are then normalized with Xo3 = κ Xo2 and Xe3 = κ−1 Xe2. For a multi-level decomposition, the algorithm above is repeated with X = Xo3. The numbers α, β, δ, γ, κ are irrational values approximately equal to
α ≈ −1.58613432,
β ≈ −0.05298011854,
δ ≈ 0.8829110762,
γ ≈ 0.4435068522,
κ ≈ 1.149604398.
The inverse transform is done by performing the lifting steps in the reverse order and with α, β, δ, γ negated.
What if X has odd length 2N-1? The trick is to extrapolate one extra element X(2N)=x in such a way that transforming the augmented X has Xe3(N)=0. This zero element can then be thrown away without losing information. The result is a decomposition with N elements in Xo3 and N-1 elements in Xe3 for a total of 2N-1 elements; the decomposition is nonredundant. To invert an odd-length transform, append the zero element Xe3(N)=0 and proceed with the usual even-length inverse transform.
The formula for the extrapolated element x such that Xe3(N)=0 is a linear function of the elements of X that depends on the filter boundary handling. For the half-sample symmetric boundary handling used in
waveletcdf97, the extrapolation is
x = (−2/(1 + 2βδ)) [αβδ X(2N-3) + βδ X(2N-2) + (α + 3αβδ + δ) X(2N-1)].
- M. Unser and T. Blu. "Mathematical Properties of the JPEG2000 Wavelet Filters." IEEE Trans. on Image Proc., vol. 12, no. 9, Sep. 2003.
- W. Sweldens. "The Lifting Scheme: A Construction of Second Generation Wavelets." SIAM J. Mathematical Analysis, vol. 29, no. 2, pp. 511-546, 1997.
- I. Daubechies and W. Sweldens. "Factoring Wavelet Transforms into Lifting Steps." 1996.