Fast Wavelet Transforms via matrix decomposition. The matrices below (when applied efficiently to vectors of signal data) implement the Fast Wavelet Transform (FWT). > with(linalg): > with(plots): > First, we develop the FWT for the Haar wavelet > H8 := evalm( 1/2 * matrix(8,8,[ > 1, 0, 0, 0, 1, 0, 0, 0, > 1, 0, 0, 0, -1, 0, 0, 0, > 0, 1, 0, 0, 0, 1, 0, 0, > 0, 1, 0, 0, 0,-1, 0, 0, > 0, 0, 1, 0, 0, 0, 1, 0, > 0, 0, 1, 0, 0, 0,-1, 0, > 0, 0, 0, 1, 0, 0, 0, 1, > 0, 0, 0, 1, 0, 0, 0,-1])); > > H4 := diag( > evalm(1/2 * > matrix(4,4,[ 1, 0, 1, 0, > 1, 0, -1, 0, > 0, 1, 0, 1, > 0, 1, 0,-1])), > diag(1$4)); > > H2 := diag( > evalm(1/2 * > matrix(2,2,[1, 1, > 1, -1])), > diag(1$6)); > [1/2 , 0 , 0 , 0 , 1/2 , 0 , 0 , 0] [ ] [1/2 , 0 , 0 , 0 , -1/2 , 0 , 0 , 0] [ ] [0 , 1/2 , 0 , 0 , 0 , 1/2 , 0 , 0] [ ] [0 , 1/2 , 0 , 0 , 0 , -1/2 , 0 , 0] H8 := [ ] [0 , 0 , 1/2 , 0 , 0 , 0 , 1/2 , 0] [ ] [0 , 0 , 1/2 , 0 , 0 , 0 , -1/2 , 0] [ ] [0 , 0 , 0 , 1/2 , 0 , 0 , 0 , 1/2] [ ] [0 , 0 , 0 , 1/2 , 0 , 0 , 0 , -1/2] [1/2 0 1/2 0 0 0 0 0] [ ] [1/2 0 -1/2 0 0 0 0 0] [ ] [ 0 1/2 0 1/2 0 0 0 0] [ ] [ 0 1/2 0 -1/2 0 0 0 0] H4 := [ ] [ 0 0 0 0 1 0 0 0] [ ] [ 0 0 0 0 0 1 0 0] [ ] [ 0 0 0 0 0 0 1 0] [ ] [ 0 0 0 0 0 0 0 1] [1/2 1/2 0 0 0 0 0 0] [ ] [1/2 -1/2 0 0 0 0 0 0] [ ] [ 0 0 1 0 0 0 0 0] [ ] [ 0 0 0 1 0 0 0 0] H2 := [ ] [ 0 0 0 0 1 0 0 0] [ ] [ 0 0 0 0 0 1 0 0] [ ] [ 0 0 0 0 0 0 1 0] [ ] [ 0 0 0 0 0 0 0 1] > Haar_Wavelet_Transform := evalm( H8 &* H4 &* H2 ); Haar_Wavelet_Transform := [1/8 , 1/8 , 1/4 , 0 , 1/2 , 0 , 0 , 0] [ ] [1/8 , 1/8 , 1/4 , 0 , -1/2 , 0 , 0 , 0] [ ] [1/8 , 1/8 , -1/4 , 0 , 0 , 1/2 , 0 , 0] [ ] [1/8 , 1/8 , -1/4 , 0 , 0 , -1/2 , 0 , 0] [ ] [1/8 , -1/8 , 0 , 1/4 , 0 , 0 , 1/2 , 0] [ ] [1/8 , -1/8 , 0 , 1/4 , 0 , 0 , -1/2 , 0] [ ] [1/8 , -1/8 , 0 , -1/4 , 0 , 0 , 0 , 1/2] [ ] [1/8 , -1/8 , 0 , -1/4 , 0 , 0 , 0 , -1/2] > x := vector(8, i -> x.i ); > Hx := evalm( Haar_Wavelet_Transform &* x ): > matrix(8,1, (i,j) -> Hx[i] ); > x := [x1, x2, x3, x4, x5, x6, x7, x8] [1/8 x1 + 1/8 x2 + 1/4 x3 + 1/2 x5] [ ] [1/8 x1 + 1/8 x2 + 1/4 x3 - 1/2 x5] [ ] [1/8 x1 + 1/8 x2 - 1/4 x3 + 1/2 x6] [ ] [1/8 x1 + 1/8 x2 - 1/4 x3 - 1/2 x6] [ ] [1/8 x1 - 1/8 x2 + 1/4 x4 + 1/2 x7] [ ] [1/8 x1 - 1/8 x2 + 1/4 x4 - 1/2 x7] [ ] [1/8 x1 - 1/8 x2 - 1/4 x4 + 1/2 x8] [ ] [1/8 x1 - 1/8 x2 - 1/4 x4 - 1/2 x8] Notice that the Haar wavelet transform can be done in time O(N), since the H-transforms can be done quickly (each is almost the identity). For fun we look at the Haar wavelet transform graphically. > n := 8; > densityplot( Haar_Wavelet_Transform[(n-i+1),j], > j=1..n, i=1..n, grid=[n,n], axes=NONE); n := 8 Now: do the same analysis for the Daubechies D4 wavelet > alias( h0 = (1+sqrt(3))/4 ); > alias( h1 = (3+sqrt(3))/4 ); > alias( h2 = (3-sqrt(3))/4 ); > alias( h3 = (1-sqrt(3))/4 ); > > T8 := evalm( 1/sqrt(2) * matrix(8,8,[ > h0, h1, h2, h3, 0, 0, 0, 0, > 0, 0, h0, h1, h2, h3, 0, 0, > 0, 0, 0, 0, h0, h1, h2, h3, > h2, h3, 0, 0, 0, 0, h0, h1, > h3, -h2, h1, -h0, 0, 0, 0, 0, > 0, 0, h3, -h2, h1, -h0, 0, 0, > 0, 0, 0, 0, h3, -h2, h1, -h0, > h1, -h0, 0, 0, 0, 0, h3, -h2 > ])); > > T4 := diag( > evalm(1/sqrt(2) * > matrix(4,4,[ h0, h1, h2, h3, > h2, h3, h0, h1, > h3, -h2, h1, -h0, > h1, -h0, h3, -h2])), > diag(1$4)); > > T2 := diag( > evalm(1/sqrt(2) * > matrix(2,2,[h0+h2, h1+h3, > h3+h1,-h0-h2])), > diag(1$6)); > I, h0, h1, h2, h3 I, h0, h1, h2, h3 I, h0, h1, h2, h3 I, h0, h1, h2, h3 T8 := [1/2 sqrt(2) h0 , 1/2 sqrt(2) h1 , 1/2 sqrt(2) h2 , 1/2 sqrt(2) h3 , 0 , 0 , 0 , 0] [0 , 0 , 1/2 sqrt(2) h0 , 1/2 sqrt(2) h1 , 1/2 sqrt(2) h2 , 1/2 sqrt(2) h3 , 0 , 0] [0 , 0 , 0 , 0 , 1/2 sqrt(2) h0 , 1/2 sqrt(2) h1 , 1/2 sqrt(2) h2 , 1/2 sqrt(2) h3] [1/2 sqrt(2) h2 , 1/2 sqrt(2) h3 , 0 , 0 , 0 , 0 , 1/2 sqrt(2) h0 , 1/2 sqrt(2) h1] [1/2 sqrt(2) h3 , %1 , 1/2 sqrt(2) h1 , %2 , 0 , 0 , 0 , 0] [0 , 0 , 1/2 sqrt(2) h3 , %1 , 1/2 sqrt(2) h1 , %2 , 0 , 0] [0 , 0 , 0 , 0 , 1/2 sqrt(2) h3 , %1 , 1/2 sqrt(2) h1 , %2] [1/2 sqrt(2) h1 , %2 , 0 , 0 , 0 , 0 , 1/2 sqrt(2) h3 , %1] %1 := 1/2 sqrt(2) (-3/4 + 1/4 sqrt(3)) %2 := 1/2 sqrt(2) (-1/4 - 1/4 sqrt(3)) T4 := [1/2 sqrt(2) h0 , 1/2 sqrt(2) h1 , 1/2 sqrt(2) h2 , 1/2 sqrt(2) h3 , 0 , 0 , 0 , 0] [1/2 sqrt(2) h2 , 1/2 sqrt(2) h3 , 1/2 sqrt(2) h0 , 1/2 sqrt(2) h1 , 0 , 0 , 0 , 0] [1/2 sqrt(2) h3 , 1/2 sqrt(2) (-3/4 + 1/4 sqrt(3)) , 1/2 sqrt(2) h1 , 1/2 sqrt(2) (-1/4 - 1/4 sqrt(3)) , 0 , 0 , 0 , 0] [1/2 sqrt(2) h1 , 1/2 sqrt(2) (-1/4 - 1/4 sqrt(3)) , 1/2 sqrt(2) h3 , 1/2 sqrt(2) (-3/4 + 1/4 sqrt(3)) , 0 , 0 , 0 , 0] [0 , 0 , 0 , 0 , 1 , 0 , 0 , 0] [0 , 0 , 0 , 0 , 0 , 1 , 0 , 0] [0 , 0 , 0 , 0 , 0 , 0 , 1 , 0] [0 , 0 , 0 , 0 , 0 , 0 , 0 , 1] [1/2 sqrt(2) , 1/2 sqrt(2) , 0 , 0 , 0 , 0 , 0 , 0] [ ] [1/2 sqrt(2) , - 1/2 sqrt(2) , 0 , 0 , 0 , 0 , 0 , 0] [ ] [0 , 0 , 1 , 0 , 0 , 0 , 0 , 0] [ ] [0 , 0 , 0 , 1 , 0 , 0 , 0 , 0] T2 := [ ] [0 , 0 , 0 , 0 , 1 , 0 , 0 , 0] [ ] [0 , 0 , 0 , 0 , 0 , 1 , 0 , 0] [ ] [0 , 0 , 0 , 0 , 0 , 0 , 1 , 0] [ ] [0 , 0 , 0 , 0 , 0 , 0 , 0 , 1] > # D8 is the discrete Wavelet transform matrix > > D8 := evalm( T2 &* T4 &* T8 ); > D8 := [%5 , %4 , %5 , %4 , %5 , %4 , %5 , %4] [ [1/2 %3 sqrt(2) h0 + 1/2 %2 sqrt(2) h2 , 1/2 %3 sqrt(2) h1 + 1/2 %2 sqrt(2) h3 , 2 1/2 %3 sqrt(2) h2 + 1/2 h0 sqrt(2) , 1/2 %3 sqrt(2) h3 + 1/2 h0 sqrt(2) h1 , 1/2 h0 sqrt(2) h2 + 1/2 h3 sqrt(2) h0 , 1/2 h3 sqrt(2) h0 + 1/2 h3 sqrt(2) h1 , 1/2 h3 sqrt(2) h2 + 1/2 %2 sqrt(2) h0 , 2 ] 1/2 h3 sqrt(2) + 1/2 %2 sqrt(2) h1] [ [1/2 h3 h0 + 1/2 %2 h2 , 1/2 h3 h1 + 1/2 %2 h3 , 2 1/2 h3 h2 + 1/2 %1 h0 , 1/2 h3 + 1/2 %1 h1 , 2 1/2 %1 h2 + 1/2 h1 h0 , 1/2 %1 h3 + 1/2 h1 , ] 1/2 h1 h2 + 1/2 %2 h0 , 1/2 h3 h1 + 1/2 %2 h1] [ 2 [1/2 %1 h2 + 1/2 h1 h0 , 1/2 %1 h3 + 1/2 h1 , 1/2 h1 h2 + 1/2 %2 h0 , 1/2 h3 h1 + 1/2 %2 h1 , 1/2 h3 h0 + 1/2 %2 h2 , 1/2 h3 h1 + 1/2 %2 h3 , 2 ] 1/2 h3 h2 + 1/2 %1 h0 , 1/2 h3 + 1/2 %1 h1] [1/2 sqrt(2) h3 , 1/2 sqrt(2) %1 , 1/2 sqrt(2) h1 , 1/2 sqrt(2) %2 , 0 , 0 , 0 , 0] [0 , 0 , 1/2 sqrt(2) h3 , 1/2 sqrt(2) %1 , 1/2 sqrt(2) h1 , 1/2 sqrt(2) %2 , 0 , 0] [0 , 0 , 0 , 0 , 1/2 sqrt(2) h3 , 1/2 sqrt(2) %1 , 1/2 sqrt(2) h1 , 1/2 sqrt(2) %2] [1/2 sqrt(2) h1 , 1/2 sqrt(2) %2 , 0 , 0 , 0 , 0 , 1/2 sqrt(2) h3 , 1/2 sqrt(2) %1] %1 := -3/4 + 1/4 sqrt(3) %2 := -1/4 - 1/4 sqrt(3) %3 := -1/4 + 1/4 sqrt(3) %4 := 1/4 sqrt(2) h1 + 1/4 sqrt(2) h3 %5 := 1/4 sqrt(2) h0 + 1/4 sqrt(2) h2 Notice that: 1. D8 is an orthogonal matrix (its transpose is its inverse). 2. The wavelet transform D_N can be done in time O(N), since the T-transforms can be done quickly (each is almost the identity). > map(simplify,evalm( inverse(D8) - transpose(D8) )); > [0 0 0 0 0 0 0 0] [ ] [0 0 0 0 0 0 0 0] [ ] [0 0 0 0 0 0 0 0] [ ] [0 0 0 0 0 0 0 0] [ ] [0 0 0 0 0 0 0 0] [ ] [0 0 0 0 0 0 0 0] [ ] [0 0 0 0 0 0 0 0] [ ] [0 0 0 0 0 0 0 0] >