Compatibility with Matlab

[maroueneboubakri]

Hi all,

Can you please tell me why I get different results than Matlab ?

Python code:

cA, (cH, cV, cD) = pywt.dwt2(np.ones((6,6)), 'db8')

Matlab code:

[cA,cH,cV,cD] = dwt2(ones(6,6), 'db8')

cH, cV and cD are all different in both matlab and python !

Thanks

[grlee77]

I am not at a computer with Matlab at the moment to check, but my first guess is that the edge mode is different. I think Matlab defaults to periodization (per in Matlab) unless otherwise set by dwtmode. PyWavelets defaults to symmetric. Please set equivalent edge modes in both as in the following table: http://pywavelets.readthedocs.io/en/latest/ref/signal-extension-modes.html#naming-conventions

Also, the shorter mode names such as per used in older versions of PyWavelets have been deprecated in favor of the longer forms such as periodic that are given in the table linked above.

[maroueneboubakri]

I tried all modes, I think there is a bug somewhere in the code, by trying another library (not maintained) it gave me the same results as matlab. This confirms that there is a bug or at least it is not compatible with matlab.

[grlee77]

Hi @maroueneboubakri

I took a closer look in Matlab today. I believe the difference you are seeing is just floating point rounding error. The true values for cH, cV and cD are all zeros for your case of an input that is all ones. Did you try this on some realistic data and see a difference?

I have pasted below identical results in obtained in Matlab R2012a and PyWavelets 0.5.1 when using a 6x6 matrix of non-zero values.

Matlab R2012a output of dwt2:

>>>> r = randn(6, 6)

r =

   0.211719675094344  -0.976407116720665  -0.334855681324846  -0.523694312367449  -0.545694250680012   3.055198864659838
   0.390594421003066  -0.012076616540741  -3.524311861772401  -0.184267189146077  -0.163427296580157   1.174558525849007
  -0.582822801284385   0.586094540238986  -0.957788233357486  -0.706594542625819  -0.768291143390094  -0.750142849593916
   0.223084337845856   1.391196012249572  -2.385548957777934   0.395635844106088  -0.861338231164588  -0.606799666375464
   0.640659641580953   0.635862433121991   1.097195760392661  -1.826871270548409   1.366671753182628   0.521357232510260
   0.914198332219636  -0.913839313340009  -0.352498204266471   0.328620393417846   2.187676486818980  -0.128570343799134

>> save -v6 /tmp/r.mat r   % save for loading in python
>> [cA,cH,cV,cD] = dwt2(r, 'db8', 'mode', 'per')

cA =

  -1.695066438786154  -2.151035980837253   1.425948711981317
   0.410443087074231  -1.229348967072931  -0.357657945353764
  -0.412286470461064   0.328895544887587   2.687350644386791

cH =

  -0.950386599624887   0.678280037233687   0.389100649514621
  -0.549876208358157  -1.088685891863638   1.176274338746321
  -0.631502648694916   0.384405263723653  -0.541964453744911

cV =

   0.981809262184676   1.514683079879554   1.346412881077456
   0.317656806665112   1.887402793317245  -1.518929853295046
   0.054341111864084  -1.582765628489000  -0.548592013923474

cD =

  -1.236825602016070   0.842977981568511  -1.507712780588098
   0.362325143657893   1.774087502304638   0.951467082281602
  -1.040858019200054   1.843153932937395   0.575394939870273

pywt output

import pywt
import scipy.io
r = scipy.io.loadmat('/tmp/r.mat')['r']
cA, (cH, cV, cD) = pywt.dwt2(r, 'db8', mode='periodization')

cA
Out[161]: 
array([[-1.69506644, -2.15103598,  1.42594871],
       [ 0.41044309, -1.22934897, -0.35765795],
       [-0.41228647,  0.32889554,  2.68735064]])

cH
Out[162]: 
array([[-0.9503866 ,  0.67828004,  0.38910065],
       [-0.54987621, -1.08868589,  1.17627434],
       [-0.63150265,  0.38440526, -0.54196445]])

cV
Out[163]: 
array([[ 0.98180926,  1.51468308,  1.34641288],
       [ 0.31765681,  1.88740279, -1.51892985],
       [ 0.05434111, -1.58276563, -0.54859201]])

cD
Out[164]: 
array([[-1.2368256 ,  0.84297798, -1.50771278],
       [ 0.36232514,  1.7740875 ,  0.95146708],
       [-1.04085802,  1.84315393,  0.57539494]])

We have extensive tests on the 1D dwt routines verifying matches to Matlab for all wavelets and edge modes. I don't believe we test against Matlab explicitly for the 2D routines, but I would be surprised if it does not match as identical underlying convolutions routines are used.

I will close this for now as not a bug, but please reopen if you find a case where there is substantial disagreement beyond floating point rounding errors.

[grlee77]

put another way. Both Matlab and PyWavlets are "wrong" for the case of constant input as neither gives a result that is exactly zero. Arguably PyWavelets is doing better as the values it returns are on the order of 1e-16 for cH and cD while Matlab's are around 1e-12.