Two different value for the Wavelet Transform (Daubechies D4)

Question

I am testing this code

   protected final double sqrt_3 = Math.sqrt( 3 );
   protected final double denom = 4 * Math.sqrt( 2 );
   //
   // forward transform scaling (smoothing) coefficients
   //
   protected final double h0 = (1 + sqrt_3)/denom;
   protected final double h1 = (3 + sqrt_3)/denom; 
   protected final double h2 = (3 - sqrt_3)/denom; 
   protected final double h3 = (1 - sqrt_3)/denom;
   //
   // forward transform wavelet coefficients
   //
   protected final double g0 =  h3;
   protected final double g1 = -h2;
   protected final double g2 =  h1;
   protected final double g3 = -h0;

protected void transform( double a[], int n )
{
  if (n >= 4) {
     int i, j;
     int half = n >> 1;

 double tmp[] = new double[n];

 i = 0;
     for (j = 0; j < n-3; j = j + 2) {
        tmp[i]      = a[j]*h0 + a[j+1]*h1 + a[j+2]*h2 + a[j+3]*h3;
        tmp[i+half] = a[j]*g0 + a[j+1]*g1 + a[j+2]*g2 + a[j+3]*g3;
    i++;
     }

     tmp[i]      = a[n-2]*h0 + a[n-1]*h1 + a[0]*h2 + a[1]*h3;
     tmp[i+half] = a[n-2]*g0 + a[n-1]*g1 + a[0]*g2 + a[1]*g3;

     for (i = 0; i < n; i++) {
        a[i] = tmp[i];
     }
  }
 } // transform

to perform a Daubechies D4 wavelet transform on this discrete array:

[1,2,0,4,5,6,8,10]

the result is

  - 0 : 1.638357430415108
  - 1 : 3.6903274198537357
  - 2 : -2.6439375651698196
  - 3 : 79.01146993331695
  - 4 : 7.399237211089009
  - 5 : 0.3882285676537802
  - 6 : -39.6029588778518
  - 7 : -19.794010741818195
  - 8 : -2.1213203435596424
  - 9 : 0.0

but when I use python pywt.dwt on the same array, I get this:

import pywt
[cA, cD] = pywt.dwt([1,2,0,4,5,6,8,10], 'db4')


>>> >>> cA
array([ 7.14848277,  1.98754736,  1.9747116 ,  0.95510018,  4.90207373,
        8.72887094, 14.23995582])
>>> cD
array([-0.5373913 , -2.00492859,  0.01927609,  0.1615668 , -0.0823509 ,
       -0.32289939,  0.92816281])

Beyond the different values, one has 10 items and the other 7.

what am I missing?

Answer 1

I have never used any of these codes, also, not really sure about your question! But, maybe, this information might help you to get closer to an answer to your question:

Daubechies 4 Wiki

Daubechies Coefficients Wiki

• Before that, I think your input vector (signal) maybe too small to make Wavelet calculations come up right? Not sure though! Maybe, try something in 1x128 size.

• Maybe, Java code is Fast Wavelet Transform. Guessing based on the following methods:

Code

   /**
     Forward Daubechies D4 transform
    */
   public void daubTrans( double s[] )
   {
      final int N = s.length;
      int n;
      for (n = N; n >= 4; n >>= 1) {
         transform( s, n );
      }
   }


   /**
     Inverse Daubechies D4 transform
    */
   public void invDaubTrans( double coef[])
   {
      final int N = coef.length;
      int n;
      for (n = 4; n <= N; n <<= 1) {
         invTransform( coef, n );
      }
   }

Based on above methods, it seems this would be "Fast Wavelet Transform", which I'm also not so sure about their calculations, you might look into this link.

There are so many so-called, similar "terms" on Wavelet transforms that it might be best to go through their math to see things, and find out what the exact method might be (e.g., Discrete Wavelet Transform, Continuous Wavelet Transform, Discrete with Packet Decomposition). Every library has some terminologies and assumptions and make different calculations. You might print to see, if you would get anything close to D4 Wavelet = {−0.1830127, −0.3169873, 1.1830127, −0.6830127}; for DB4, first. Or, you may do other testing to see, if the calculations are correct.

Methods of Decomposition in Wavelets

It looks like cA and cD are coefficients of "Approximated" and "Details" signals decomposed by a discrete Wavelet transform. However, I'm not so sure, to how many layers you might have been decomposed your input vector.

There are two well-known ways of decomposing a signal in Wavelet, one is "packet" (which decomposes both "approximations" and "details" signals, so you would get 2^4=16 sub-signals for decomposing your original signal to 4 layers).

The other decomposition method decomposes the low-frequency part of signals. So, you might need to find out about your level of decomposition that your vector is being decomposed.

Also, if you write your own code, you can decompose it however you wish.

Simple Keys to Understand Wavelet

Shifting (Time) vs Scale (Frequency)

There is one simple thing that if you understand, then Wavelet becomes much easier. First, as you may know, Wavelet is a time-frequency method. However, instead of plotting time vs frequency, you do time vs scale, where scale is "inverse" of frequency.

Children of a Wavelet Function such as DB4

Wavelet transform maps a Wavelet function - such as DB4 - throughout your original signal, and that's how it would compute those numbers that you have printed out, perhaps. One thing to consider is to find a base function, DB4, that would "look like" you original signal. How do you do that?

Basically, you pick a base function, DB4, then Wavelet transform creates multiple forms of that base function (e.g., imagine you name them DB4-0, DB4-2, DB4-3, DB4-4, ...,DB4-15). These children are created based on:

(a) Shifting (in a for loop by incrementing time, sliding a child function, then calculating coefficients), shifting is in relation with time, obviously.

(b) Scaling (means "stretching" a Wavelet function, vertically or horizontally, which would changes the frequency nature of the base function, then sliding it again through time), which is converse relation with frequency, meaning that higher scale, lower frequencies, and vice versa.

Therefore, this depends on how many children functions you may need, based on the decompositions (sub-signals). If you have 16 sub-signals (4 level of decomposition with a packet method), then you will have 16 of those "children" functions mapping throughout your signal. And that's how coefficients vectors being calculated. Then, you may toss those unnecessary sub-signals, and keep focusing on those sub-signals (frequencies) that you might be interested in. The thing is Wavelet reserves (maintain) the time information, as opposed to Fourier.

Normal Decomposition

Also, since you are a good programmer, I'm pretty sure, you can quickly crack the code and I don't think you might be missing anything here. You can just go through their methods and read a few pages of wikipedia, and you would be probably there, if you wish so.

If you might have really exciting details questions, you may try DSP SE. So many signals experts are in there. Sorry about that! Wrote this too fast, also not a good writer/explainer, later hopefully others would edit and provide the right answer. Not really expert.

In short, you are not missing on anything, good method, good luck and best wishes!